diff --git a/docs/alkaline_ref.md b/docs/alkaline_ref.md new file mode 100644 index 0000000..e7904f1 --- /dev/null +++ b/docs/alkaline_ref.md @@ -0,0 +1,73 @@ +1. [Oystein Ulleberg, 2003](https://www.sciencedirect.com/science/article/pii/S0360319902000332?via%3Dihub) + "Modeling of advanced alkaline electrolyzers: a system simulation approach" + + +2. [Gambou, Guilbert,et al 2022](https://www.mdpi.com/1996-1073/15/9/3452) + "A Comprehensive Survey of Alkaline Electrolyzer Modeling: Electrical + Domain and Specific Electrolyte Conductivity" + + +3. [Haug ,Kreitz, 2017](https://www.sciencedirect.com/science/article/pii/S0360319917318633) + "Process modelling of an alkaline water electrolyzer" + + +4. [Henou, Agbossou, 2014](https://www.sciencedirect.com/science/article/pii/S0378775313017527) + - "Simulation tool based on a physics model and an electrical + analogy for an alkaline electrolyser" + - Notes: + - cited by [Gambou, Guilbert,et al 2022] + - HAS ALL THE VALUES FOR VARIABLES USED IN [Gambou, Guilbert,et al 2022] + +5. [Hammoudi,Henao, 2012](https://www.sciencedirect.com/science/article/pii/S036031991201590X) + - "New multi-physics approach for modelling andn design of + alkaline electrolyzers" + - Notes: + - Referenced by [Henou, Agbossou, 2014] for theta calculation + (electrode coverage) + - Eqn 44 for bubble stuff + - j_lim=300 kA/m^2 + - includes other efficiency losses + - cites: https://www.sciencedirect.com/science/article/pii/S0360128509000598 + +6. [Brauns,2021](https://iopscience.iop.org/article/10.1149/1945-7111/abda57/pdf) + - bibtex label: + - "Evaluation of Diaphragms and Membranes as Separators for Alkaline Water Electrolysis" by Jorn Brauns et all 2021. J. Electrochem Soc 168 014510 + - Notes: + - good numbers + - electrolyte flow rate of 350 mL/min + - total electrolyte volume of 10L + - has supplementary material (need to checkout) + - in "material stability" it mentions stuff about DEGRADATION + +7. [NEL Report](https://www.energy.gov/sites/default/files/2022-02/2-Intro-Liquid%20Alkaline%20Workshop.pdf) + - bibtex label: + + +8. [Brauns, Turek 2020](https://www.mdpi.com/2227-9717/8/2/248) + - bibtex label: + - "Alkaline Water Electrolysis Powered by Renewable Energy: A Review" + + +9. [Eigeldinger, Vogt 2000](https://www.sciencedirect.com/science/article/pii/S0013468600005132) + - bibtex label: Eigeldinger_2000 + - "The bubble coverage of gas evolving electrodes in a flowing electrolyte" + - Notes: + - Ref 15 of Henou 2014 for current density and theta stuff + - has current density equation with theta included + +10. [Haug, Koj, Turek 2017](https://www.sciencedirect.com/science/article/pii/S0360319916336588) +"Influence of process conditions on gas purity in alkaline water electrolysis" +by Phillip Haug, Motthias Koj, Thomas Turek [2017] + + +11. [Niroula, Chaudhary, Subedi, Thapa 2003](https://iopscience.iop.org/article/10.1088/1757-899X/1279/1/012005/pdf) + + "Parametric Modelling and Optimization of Alkaline Electrolyzer for the Production + of Green Hydrogen" by S. Niroula, C Chaudhary, A Subedi, and B S Thapa + [2003] doi:10.1088/1757-899X/1279/1/012005 + + +12. [Vogt,Balzer 2005](https://www.sciencedirect.com/science/article/pii/S001346860400948X?via%3Dihub) +"The bubble coverage of gas-evolving electrodes in stagnant electrolytes" +by H. Vogt and R.J. Balzer +Volume 50, Issue 10, 15 March 2005, Pages 2073-2079 diff --git a/electrolyzer/simulation/base.py b/electrolyzer/simulation/base.py new file mode 100644 index 0000000..e40daa8 --- /dev/null +++ b/electrolyzer/simulation/base.py @@ -0,0 +1,34 @@ +from typing import Any, Dict + +import attrs + +from electrolyzer.tools.type_dec import FromDictMixin + + +class BaseClass(FromDictMixin): + """ + BaseClass object class. This class does the logging and MixIn class inheritance. + """ + + @classmethod + def get_model_defaults(cls) -> Dict[str, Any]: + """Produces a dictionary of the keyword arguments and their defaults. + + Returns + ------- + Dict[str, Any] + Dictionary of keyword argument: default. + """ + return {el.name: el.default for el in attrs.fields(cls)} + + def _get_model_dict(self) -> dict: + """Convenience method that wraps the `attrs.asdict` method. Returns the object's + parameters as a dictionary. + + Returns + ------- + dict + The provided or default, if no input provided, + model settings as a dictionary. + """ + return attrs.asdict(self) diff --git a/electrolyzer/simulation/bert.py b/electrolyzer/simulation/bert.py index bc71e34..59264f7 100644 --- a/electrolyzer/simulation/bert.py +++ b/electrolyzer/simulation/bert.py @@ -105,6 +105,55 @@ def _run_electrolyzer_opt(modeling_options, power_signal): return tot_kg, max_curr_density +def run_LCA(elec_sys, plant_life_years: int): + sys_refurb = pd.DataFrame() + sys_aep = pd.DataFrame() + sys_ah2 = pd.DataFrame() + years = list(np.arange(0, plant_life_years, 1)) + for i, stack in enumerate(elec_sys.stacks): + id = i + 1 + refturb_schedule, ahp_kg, aep_kWh = stack.estimate_life_performance_from_year( + plant_life_years + ) + sys_aep = pd.concat( + [ + sys_aep, + pd.DataFrame( + dict(zip(["Stack {}".format(id)], [aep_kWh])), index=years + ), + ] + ) + sys_ah2 = pd.concat( + [ + sys_ah2, + pd.DataFrame(dict(zip(["Stack {}".format(id)], [ahp_kg])), index=years), + ] + ) + sys_refurb = pd.concat( + [ + sys_refurb, + pd.DataFrame( + dict(zip(["Stack {}".format(id)], [refturb_schedule])), index=years + ), + ] + ) + # df_index = [ + # [f"stack {id}"]*plant_life_years, + # list(np.arange(0,plant_life_years,1)) + # ] + # temp_df = pd.DataFrame(dict(zip( + # ["AEP [kWh/year]","Refurb Schedule","AHP [kg/year]"], + # [aep_kWh,ahp_kg,refturb_schedule])), + # index=df_index + # ) + # stack_LCA = pd.concat([stack_LCA,temp_df]) + return { + "Refurb Schedule": sys_refurb, + "AEP [kWh/year]": sys_aep, + "AHP [kg/year]": sys_ah2, + } + + def run_electrolyzer(input_modeling, power_signal, optimize=False): """ Runs an electrolyzer simulation based on a YAML configuration file and power diff --git a/electrolyzer/simulation/cell.py b/electrolyzer/simulation/cell.py deleted file mode 100644 index baeac90..0000000 --- a/electrolyzer/simulation/cell.py +++ /dev/null @@ -1 +0,0 @@ -# This will contain the baseclass for different types of cells diff --git a/electrolyzer/simulation/cell_models/alkaline.py b/electrolyzer/simulation/cell_models/alkaline.py index a16c6ef..1825446 100644 --- a/electrolyzer/simulation/cell_models/alkaline.py +++ b/electrolyzer/simulation/cell_models/alkaline.py @@ -15,87 +15,10 @@ from scipy.constants import R, physical_constants, convert_temperature from electrolyzer.tools.type_dec import FromDictMixin +from electrolyzer.tools.validators import range_val warnings.filterwarnings("ignore") -""" -[Oystein Ulleberg, 2003] - "Modeling of advanced alkaline electrolyzers: a system simulation approach" - https://www.sciencedirect.com/science/article/pii/S0360319902000332?via%3Dihub - -[Gambou, Guilbert,et al 2022] - "A Comprehensive Survey of Alkaline Eelctrolyzer Modeling: Electrical - Domain and Specific Electrolyte Conductivity" - https://www.mdpi.com/1996-1073/15/9/3452 - -[Haug,Kreitz, 2017] - "Process modelling of an alkaline water electrolyzer" - https://www.sciencedirect.com/science/article/pii/S0360319917318633 - -[Henou, Agbossou, 2014] - "Simulation tool based on a physics model and an electrical - analogy for an alkaline electrolyser" - https://www.sciencedirect.com/science/article/pii/S0378775313017527 - ->cited by [Gambou, Guilbert,et al 2022] - -> HAS ALL THE VALUES FOR VARIABLES USED IN [Gambou, Guilbert,et al 2022] - -[Hammoudi,Henao, 2012] - "New multi-physics approach for modelling andn design of - alkaline electrolyzers" - https://www.sciencedirect.com/science/article/pii/S036031991201590X - -> Referenced by [Henou, Agbossou, 2014] for theta calculation - (electrode coverage) - ->Eqn 44 for bubble stuff - ->j_lim=300 kA/m^2 - ->includes other efficiency losses - cites: - https://www.sciencedirect.com/science/article/pii/S0360128509000598 - -[Brauns,2021] -"Evaluation of Diaphragms and Membranes as Separators for Alkaline -Water Electrolysis" - by Jorn Brauns et all 2021. J. Electrochem Soc 168 014510 - https://iopscience.iop.org/article/10.1149/1945-7111/abda57/pdf - ->good numbers - ->electrolyte flow rate of 350 mL/min - ->total electrolyte volume of 10L - -> has supplementary material (need to checkout) - ->in "material stability" it mentions stuff about DEGRADATION - -NEL Report: -https://www.energy.gov/sites/default/files/2022-02/2-Intro-Liquid%20Alkaline%20Workshop.pdf - -[Brauns, Turek 2020] -"Alkaline Water Electrolysis Powered by Renewable Energy: A Review" - https://www.mdpi.com/2227-9717/8/2/248 - - -[Eigeldinger, Vogt 2000] -"The bubble coverage of gas evolving electrodes in a flowing electrolyte" -https://www.sciencedirect.com/science/article/pii/S0013468600005132 - -> Ref 15 of Henou 2014 for current density and theta stuff - -> has current density equation with theta included - -[Haug, Koj, Turek 2017] -"Influence of process conditions on gas purity in alkaline water electrolysis" -by Phillip Haug, Motthias Koj, Thomas Turek [2017] -https://www.sciencedirect.com/science/article/pii/S0360319916336588 - -[Niroula, Chaudhary, Subedi, Thapa 2003] -"Parametric Modelling and Optimization of Alkaline Electrolyzer for the Production -of Green Hydrogen" by S. Niroula, C Chaudhary, A Subedi, and B S Thapa -[2003] doi:10.1088/1757-899X/1279/1/012005 -https://iopscience.iop.org/article/10.1088/1757-899X/1279/1/012005/pdf - -[Vogt,Balzer 2005] -"The bubble coverage of gas-evolving elecrodes in stagnant electrolytes" -by H. Vogt and R.J. Balzer -Volume 50, Issue 10, 15 March 2005, Pages 2073-2079 -https://www.sciencedirect.com/science/article/pii/S001346860400948X?via%3Dihub - - - -""" def ael_electrolyzer_model(X, a, b, c, d, e, f): @@ -129,11 +52,13 @@ class AlkalineCell(FromDictMixin): turndown_ratio: float max_current_density: float - cell_area: float = field(init=False) + f_1: float # faradaic coefficient in mA^2/cm^4 + f_2: float = field(validator=range_val(0, 1)) # faradaic coefficient [unitless] + + cell_area: float = field(init=False) # cell area in cm^2 # Electrode parameters # #################### - A_electrode: float = field(init=False) # [cm^2] e_a: float = field(init=False) # [cm] anode thickness e_c: float = field(init=False) # [cm] cathode thickness d_am: float = field(init=False) # [cm] Anode-membrane gap @@ -151,7 +76,6 @@ class AlkalineCell(FromDictMixin): # Membrane parameters # #################### - A_membrane: float = field(init=False) # [cm^2] e_m: float = field(init=False) # [cm] membrane thickness # THIS ONE IS PRIMARLY BASED ON @@ -166,6 +90,7 @@ class AlkalineCell(FromDictMixin): M_K: float = 39.0983 # molecular weight of Potassium [g/mol] lhv: float = 33.33 # lower heating value of H2 [kWh/kg] hhv: float = 39.41 # higher heating value of H2 [kWh/kg] + gibbs: float = 237.24e3 # Gibbs Energy of global reaction (J/mol) def __attrs_post_init__(self) -> None: # Cell parameters # @@ -174,12 +99,12 @@ def __attrs_post_init__(self) -> None: # Electrode parameters # ######################## - self.A_electrode = self.electrode["A_electrode"] - self.e_a = self.electrode["e_a"] - self.e_c = self.electrode["e_c"] - self.d_am = self.electrode["d_am"] - self.d_cm = self.electrode["d_cm"] - self.d_ac = self.electrode["d_ac"] + self.e_a = self.electrode["e_e"] # anode thickness + self.e_c = self.electrode["e_e"] # cathode thickness + + self.d_am = self.electrode["d_em"] # anode-membrane distance + self.d_cm = self.electrode["d_em"] # cathode-membrane distance + self.d_ac = self.electrode["d_ac"] # distance between electrodes # Electrolyte parameters # ########################## @@ -187,7 +112,6 @@ def __attrs_post_init__(self) -> None: # Membrane parameters # ####################### - self.A_membrane = self.membrane["A_membrane"] self.e_m = self.membrane["e_m"] # calcluate molarity and molality of KOH solution @@ -214,7 +138,7 @@ def calculate_bubble_rate_coverage(self, T_C, I): T_k = convert_temperature([T_C], "C", "K")[0] J_lim = 30 # [A/cm^2] [Vogt,Balzer 2005] T_amb = T_k = convert_temperature([25], "C", "K")[0] - j = I / self.A_electrode # [A/cm^2] "nominal current density" + j = I / self.cell_area # [A/cm^2] "nominal current density" # Eqn 19 of [Gambou, Guilbert,et al 2022] Pv_H20 = np.exp( @@ -247,7 +171,8 @@ def calc_current_density(self, T_C, I): # actual current density reflecting impact of bubble rate coverage theta_epsilon = self.calculate_bubble_rate_coverage(T_C, I) theta = theta_epsilon[0] - A_electrode_eff = self.A_electrode * (1 - theta) # [cm^2] + # A_electrode_eff = self.A_electrode * (1 - theta) # [cm^2] + A_electrode_eff = self.cell_area * (1 - theta) # [cm^2] current_dens = I / A_electrode_eff # [A/cm^2] return current_dens @@ -369,7 +294,8 @@ def calc_Urev(self, T_C, P): # Eqn 16 Pv_KOH = np.exp((2.302 * a) + (b * np.log(Pv_H20))) - Urev0 = self.calc_open_circuit_voltage(T_C) + Urev0 = self.calc_reversible_voltage() + # Urev0 = self.calc_open_circuit_voltage(T_C) # Eqn 14 U_rev = Urev0 + ((R * T_K) / (self.z * F)) * np.log( @@ -452,13 +378,13 @@ def calc_electrode_resistance(self, T_C): # Eqn 20 - resistivity of anode Ra = ( rho_nickle_eff - * (self.e_a / self.A_electrode) + * (self.e_a / self.cell_area) * (1 + (temp_coeff * (T_C - tref))) ) # Eqn 20 - resistivity of cathode Rc = ( rho_nickle_eff - * (self.e_c / self.A_electrode) + * (self.e_c / self.cell_area) * (1 + (temp_coeff * (T_C - tref))) ) @@ -504,9 +430,8 @@ def calc_electrolyte_resistance(self, T_C, I): # R_ele_bf: Bubble-free electrolyte resistance # Eqn 32 of [Gambou, Guilbert,et al 2022] and Eqn 19 of [Henou, Agbossou, 2014] - R_ele_bf = (100 / sigma_bf) * ( - (self.d_am / self.A_electrode) + (self.d_cm / self.A_electrode) + (self.d_am / self.cell_area) + (self.d_cm / self.cell_area) ) # Eqn 2 of [Brauns,2021] says R_eg=(1/(sigma_koh))*((dcs+das)/A_el) # where A_el is the metal area, not the entire area (which has holes) @@ -514,7 +439,6 @@ def calc_electrolyte_resistance(self, T_C, I): # Resistance due to bubbles theta_epsilon = self.calculate_bubble_rate_coverage(T_C, I) epsilon = theta_epsilon[1] - # R_ele_b=R_ele_bf*((1/(1-epsilon)**(3/2))-1) # R_ele_b: Bubble resistance R_ele_b = R_ele_bf * ((1 / ((1 - epsilon) ** (3 / 2))) - 1) # ^Bruggman equation @@ -531,7 +455,7 @@ def calc_membrane_resistance(self, T_C): Reference: [Gambou, Guilbert,et al 2022]: Eqn 36 [Henou, Agbossou, 2014]: Eqn 21 - [NEL]: TODO add slide + [NEL]: TODO add reference to slide """ # NOTE: THIS HAS BEEN VERIFIED @@ -540,8 +464,9 @@ def calc_membrane_resistance(self, T_C): # [Gambou, Guilbert,et al 2022] # S_mem=54.48 # membrane surface area in cm^2 + # Equation 36 - Ohms Rmem = (0.06 + 80 * np.exp(T_C / 50)) / ( - 10000 * self.A_membrane + 10000 * self.cell_area ) # Equation 36 - Ohms # ^ Equation 21 of [Henou, Agbossou, 2014] # output: Rmem=0.23 ohm*cm^2 @@ -595,29 +520,11 @@ def calc_ohmic_overpotential(self, T_C, I): V_ohm = I * R_tot # [V/cell] return V_ohm - def calc_open_circuit_voltage(self, T_C): + def calc_reversible_voltage(self): """ - I [A]: current - T_C [C]: temperature - return :: E_rev0 [V/cell]: open-circuit voltage - - TODO: Are we correcting for temperature twice? U_rev0 should be just 1.229 and - never change (possibly?) - - Reference: [Gambou, Guilbert,et al 2022]: Eqn 14 + Calculates reversible cell potential at standard state. """ - # General Nerst Equation - # Eqn 14 of [Gambou, Guilbert,et al 2022] - T_K = convert_temperature([T_C], "C", "K")[0] - E_rev0 = ( - 1.5184 - - (1.5421 * (10 ** (-3)) * T_K) - + (9.523 * (10 ** (-5)) * T_K * np.log(T_K)) - + (9.84 * (10 ** (-8)) * (T_K**2)) - ) - # OR should this just be 1.229? - # E_rev_fake = 1.229 - return E_rev0 + return self.gibbs / (self.z * F) def calc_faradaic_efficiency(self, T_C, I): """ @@ -627,8 +534,8 @@ def calc_faradaic_efficiency(self, T_C, I): Reference: [Oystein Ulleberg, 2003] Eqn 9, Table 3 """ # f1 and f2 values from Table 3 - f1 = 250 # [mA^2/cm^4] - f2 = 0.96 # [-] + # f1 = 250 # [mA^2/cm^4] + # f2 = 0.96 # [-] j = self.calc_current_density(T_C, I) # [A/cm^2] j *= 1000 # [mA/cm^2] @@ -640,7 +547,7 @@ def calc_faradaic_efficiency(self, T_C, I): # f_2=[0.99,0.985,0.98,0.96] #[0-1] # Eqn 9 from [Oystein Ulleberg, 2003] - eta_F = f2 * (j**2) / (f1 + j**2) + eta_F = self.f_2 * (j**2) / (self.f_1 + j**2) return eta_F def calc_mass_flow_rate(self, T_C, I): @@ -653,10 +560,9 @@ def calc_mass_flow_rate(self, T_C, I): eta_F = self.calc_faradaic_efficiency(T_C, I) # Eqn 10 [mol/sec] - h2_prod_mol = eta_F * I / (self.z * F) + h2_prod_mol = eta_F * I / (self.z * F) # mol/s mfr = self.M_H * self.z * h2_prod_mol # [g/sec] # z is valency number of electrons transferred per ion # for oxygen, z=4 mfr = mfr / 1e3 # [kg/sec] return mfr - # h2_prod is in mol/s diff --git a/electrolyzer/simulation/cell_models/cell.py b/electrolyzer/simulation/cell_models/cell.py new file mode 100644 index 0000000..9c5f7ad --- /dev/null +++ b/electrolyzer/simulation/cell_models/cell.py @@ -0,0 +1 @@ +# this file will contain a cell base-class diff --git a/electrolyzer/simulation/cell_models/pem.py b/electrolyzer/simulation/cell_models/pem.py index 8ea8407..e0c892c 100644 --- a/electrolyzer/simulation/cell_models/pem.py +++ b/electrolyzer/simulation/cell_models/pem.py @@ -42,7 +42,17 @@ class PEMCell(FromDictMixin): cell_area: float turndown_ratio: float max_current_density: float - + # everything below is new + p_anode: float # bar + p_cathode: float # bar + alpha_a: float + alpha_c: float + i_0_a: float + i_0_c: float + e_m: float + R_ohmic_elec: float + f_1: float # faradaic coefficient in mA^2/cm^4 + f_2: float # faradaic coefficient [unitless] # If we rework this class to be even more generic, we can have these be specified # as configuration params @@ -64,8 +74,8 @@ def calc_open_circuit_voltage(self, temperature): """Calculates open circuit voltage using the Nernst equation.""" T_K = convert_temperature([temperature], "C", "K")[0] E_rev_0 = self.calc_reversible_voltage() - p_anode = P_ATMO # (Pa) assumed atmo - p_cathode = 30 * P_STD # (Pa) 30 bars + p_anode = self.p_anode * P_STD # (Pa) + p_cathode = self.p_cathode * P_STD # (Pa) # noqa: E501 # Arden Buck equation T=C, https://www.omnicalculator.com/chemistry/vapour-pressure-of-water#vapor-pressure-formulas # noqa @@ -83,6 +93,9 @@ def calc_open_circuit_voltage(self, temperature): p_o2 = p_anode - p_h2O_sat # General Nernst equation, 10.1016/j.ijhydene.2017.03.046 + # E_cell = E_rev_0 + ((R * T_K) / (self.n * F)) * ( + # np.log((p_h2 / p_h2O_sat) * np.sqrt(p_o2 / p_h2O_sat)) + # ) E_cell = E_rev_0 + ((R * T_K) / (self.n * F)) * ( np.log((p_h2 / P_ATMO) * np.sqrt(p_o2 / P_ATMO)) ) @@ -105,20 +118,24 @@ def calc_activation_overpotential(self, i, temperature): T_cathode = T_K # TODO: updated with realistic anode temperature? # anode charge transfer coefficient TODO: is this a realistic value? - alpha_a = 2 + # alpha_a = 2 # cathode charge transfer coefficient TODO: is this a realistic value? - alpha_c = 0.5 + # alpha_c = 0.5 # anode exchange current density TODO: update to be f(T)? - i_0_a = 2e-7 + # i_0_a = 2e-7 # cathode exchange current density TODO: update to be f(T)? - i_0_c = 2e-3 + # i_0_c = 2e-3 # derived from Butler-Volmer eqs - V_act_a = ((R * T_anode) / (alpha_a * F)) * np.arcsinh(i / (2 * i_0_a)) - V_act_c = ((R * T_cathode) / (alpha_c * F)) * np.arcsinh(i / (2 * i_0_c)) + V_act_a = ((R * T_anode) / (self.alpha_a * F)) * np.arcsinh( + i / (2 * self.i_0_a) + ) + V_act_c = ((R * T_cathode) / (self.alpha_c * F)) * np.arcsinh( + i / (2 * self.i_0_c) + ) # alternate equations for Activation overpotential # Option 2: Dakota: I believe this may be more accurate, found more @@ -144,16 +161,16 @@ def calc_ohmic_overpotential(self, i, temperature): # pulled from https://www.sciencedirect.com/science/article/pii/S0360319917309278?via%3Dihub # noqa # TODO: pulled from empirical data, is there a better eq? lambda_nafion = ((-2.89556 + (0.016 * T_K)) + 1.625) / 0.1875 - t_nafion = 0.02 # (cm) confirmed that membrane thickness is <0.02. + # t_nafion = 0.02 # (cm) confirmed that membrane thickness is <0.02. # TODO: confirm with Nel, is there a better eq? sigma_nafion = ((0.005139 * lambda_nafion) - 0.00326) * np.exp( 1268 * ((1 / 303) - (1 / T_K)) ) - R_ohmic_ionic = t_nafion / sigma_nafion + R_ohmic_ionic = self.e_m / sigma_nafion # TODO: confirm realistic value with Nel https://www.sciencedirect.com/science/article/pii/S0378775315001901 # noqa - R_ohmic_elec = 50e-3 + # R_ohmic_elec = 50e-3 # Alternate R_ohmic_elec from https://www.sciencedirect.com/science/article/pii/S0360319918309017 # noqa # rho = (ohm*m) material resistivity @@ -161,7 +178,7 @@ def calc_ohmic_overpotential(self, i, temperature): # A_path = (m2) cross-sectional area of conductor path # R_ohmic_elec = ((rho*l_path)/A_path) - V_ohmic = i * (R_ohmic_elec + R_ohmic_ionic) + V_ohmic = i * (self.R_ohmic_elec + R_ohmic_ionic) return V_ohmic @@ -225,25 +242,23 @@ def calc_faradaic_efficiency(self, T_C, I): return :: eta_F [-]: Faraday's efficiency Reference: https://res.mdpi.com/d_attachment/energies/energies-13-04792/article_deploy/energies-13-04792-v2.pdf """ # noqa - f_1 = 250 # (mA2/cm4) - f_2 = 0.996 + I *= 1000 eta_F = ( - ((I / self.cell_area) ** 2) / (f_1 + ((I / self.cell_area) ** 2)) - ) * f_2 + ((I / self.cell_area) ** 2) / (self.f_1 + ((I / self.cell_area) ** 2)) + ) * self.f_2 return eta_F - def calc_mass_flow_rate(self, T_C, Idc, dryer_loss=6.5): + def calc_mass_flow_rate(self, T_C, Idc): """ Idc [A]: stack current - dryer_loss [%]: loss of drying H2 T_C [C]: cell temperature (currently unused) return :: mfr [kg/s]: mass flow rate """ eta_F = self.calc_faradaic_efficiency(T_C, Idc) - mfr = eta_F * Idc * self.M / (self.n * F) * (1 - dryer_loss / 100.0) # [g/s] + mfr = eta_F * Idc * self.M / (self.n * F) # [g/s] # mfr = mfr / 1000. * 3600. # [kg/h] mfr = mfr / 1e3 # [kg/s] return mfr diff --git a/electrolyzer/simulation/cluster.py b/electrolyzer/simulation/cluster.py index e69de29..fcfd460 100644 --- a/electrolyzer/simulation/cluster.py +++ b/electrolyzer/simulation/cluster.py @@ -0,0 +1 @@ +# this file will contain the cluster class diff --git a/electrolyzer/simulation/stack.py b/electrolyzer/simulation/stack.py index 9087569..1283756 100644 --- a/electrolyzer/simulation/stack.py +++ b/electrolyzer/simulation/stack.py @@ -8,8 +8,10 @@ import rainflow from attrs import field, define from scipy.signal import tf2ss, cont2discrete +from scipy.constants import physical_constants from electrolyzer.tools.type_dec import NDArrayFloat, FromDictMixin, array_converter +from electrolyzer.tools.validators import contains from electrolyzer.simulation.cell_models.pem import PEMCell, PEM_electrolyzer_model from electrolyzer.simulation.cell_models.alkaline import ( AlkalineCell, @@ -17,12 +19,15 @@ ) +F, _, _ = physical_constants["Faraday constant"] # Faraday's constant [C/mol] + + @define class Stack(FromDictMixin): # Stack parameters # #################### dt: float = field() - cell_type: str = field() + cell_type: str = field(validator=contains(["PEM", "alkaline"])) temperature: float = field() n_cells: int = field() @@ -31,6 +36,12 @@ class Stack(FromDictMixin): stack_rating_kW: float = field(default=None) include_degradation_penalty: bool = field(default=True) + # If degradation results in hydrogen losses, hydrogen_degradation_penalty=True + # If degradation results in more power consumed, hydrogen_degradation_penalty=False + hydrogen_degradation_penalty: bool = field(default=True) + eol_eff_percent_loss: float = field(init=False) + d_eol: float = field(init=False) + max_current: float = field(default=1000) # TODO this is a bad default, fix later min_power: float = field(default=None) @@ -74,6 +85,12 @@ class Stack(FromDictMixin): voltage_history: NDArrayFloat = field( init=False, default=[], converter=array_converter ) + degradation_history: NDArrayFloat = field( + init=False, default=[], converter=array_converter + ) + power_input_history: NDArrayFloat = field( + init=False, default=[], converter=array_converter + ) # [V] degradation from fluctuating power only d_f: float = field(init=False, default=0) @@ -116,7 +133,9 @@ class Stack(FromDictMixin): # [s] total time of simulation time: float = field(init=False, default=0) - # [s] time constant https://www.sciencedirect.com/science/article/pii/S0360319911021380 section 3.4 # noqa + # [s] time constant + # https://www.sciencedirect.com/science/article/pii/S0360319911021380 + # section 3.4 # noqa tau: float = 5 stack_state: float = field(init=False, default=0) @@ -130,7 +149,7 @@ class Stack(FromDictMixin): def __attrs_post_init__(self) -> None: # Stack parameters # #################### - + self.eol_eff_percent_loss = self.degradation["eol_eff_percent_loss"] if self.cell_type == "PEM": # initialize electrolzyer cell model self.cell = PEMCell.from_dict(self.cell_params["PEM_params"]) @@ -194,34 +213,107 @@ def __attrs_post_init__(self) -> None: ] ) - # self.h2_pres_out = 31 # H2 outlet pressure (bar) - self.DTSS = self.calc_state_space() + self.d_eol = self.calc_end_of_life_voltage() - def run(self, P_in): + def run_power_deg_penalty(self, P_in): + """Run the stack with a degradation penalty applied to the power consumption. + Power consumed may be greater than power input. + + Args: + P_in (float): stack power input in Wdc + + Returns: + 2-element tuple containing + + - **I_stack** (float): stack current in Amps + - **V_cell** (float): cell voltage in Volts + """ + + I_stack = self.electrolyzer_model( + (P_in / 1e3, self.temperature), *self.fit_params + ) + V_cell = self.cell.calc_cell_voltage(I_stack, self.temperature) + + return I_stack, V_cell + + def run_h2_deg_penalty(self, P_in): + """Run the stack with a degradation penalty applied to the hydrogen production. + Power consumed should be nearly equal to power input. + + Args: + P_in (float): stack power input in Wdc + + Returns: + 2-element tuple containing + + - **I_stack** (float): stack current in Amps + - **V_cell** (float): cell voltage in Volts """ - P_in [Wdc]: stack power input - return :: H2_mfr [kg/s]: hydrogen mass flow rate - return :: H2_mass_out [kg]: hydrogen mass - return :: power_left [W]: difference in P_in and power consumed + + I_nom = self.electrolyzer_model( + (P_in / 1e3, self.temperature), *self.fit_params + ) + V_cell = self.cell.calc_cell_voltage(I_nom, self.temperature) + eff_mult = np.nan_to_num((V_cell + self.V_degradation) / V_cell) + I_stack = np.nan_to_num(I_nom / eff_mult) + + return I_stack, V_cell + + def run(self, P_in): + """Run the stack for smoe input power. + + Args: + P_in (float): stack power input in Wdc + + Returns: + 3-element tuple containing + + - **H2_mfr** (float): hydrogen mass flow rate in kg/sec + - **H2_mass_out** (float): hydrogen mass produced in kg/dt + - **power_left** (float): remaining power, P_in and power consumed + """ self.update_status() + if self.hydrogen_degradation_penalty: + I_stack, V_cell = self.run_h2_deg_penalty(P_in) + else: + I_stack, V_cell = self.run_power_deg_penalty(P_in) + + H2_mfr, H2_mass_out, power_left = self.run_stack(I_stack, V_cell, P_in) + return H2_mfr, H2_mass_out, power_left - I = self.electrolyzer_model((P_in / 1e3, self.temperature), *self.fit_params) - V = self.cell.calc_cell_voltage(I, self.temperature) + def run_stack(self, I_stack, V_cell, P_in): + """Run the stack for a given stack current, cell voltage, and power input. + Updates the cell degradation, updates dynamics, calculates hydrogen mass + flow rate and production, and remaining power. + + Args: + I_stack (float): stack current input in Amps + V_cell (float): cell voltage in Volts + P_in (float): stack power input in Wdc + Returns: + 3-element tuple containing + + - **H2_mfr** (float): hydrogen mass flow rate in kg/sec + - **H2_mass_out** (float): hydrogen mass produced in kg/dt + - **power_left** (float): difference in P_in and power consumed in Watts + """ if self.stack_on: power_left = P_in - self.I = I + self.I = I_stack if self.include_degradation_penalty: - V += self.V_degradation + V_cell += self.V_degradation - self.update_temperature(I, V) + self.update_temperature(I_stack, V_cell) self.update_degradation() - power_left -= self.calc_stack_power(I, V) * 1e3 - H2_mfr = self.cell.calc_mass_flow_rate(self.temperature, I) * self.n_cells + power_left -= self.calc_stack_power(I_stack, V_cell) * 1e3 + H2_mfr = ( + self.cell.calc_mass_flow_rate(self.temperature, I_stack) * self.n_cells + ) self.stack_state, H2_mfr = self.update_dynamics(H2_mfr, self.stack_state) H2_mass_out = H2_mfr * self.dt @@ -234,17 +326,20 @@ def run(self, P_in): if self.stack_waiting: self.uptime += self.dt - self.I = I - self.update_temperature(I, V) + self.I = I_stack + self.update_temperature(I_stack, V_cell) self.update_degradation() power_left = 0 else: power_left = P_in - V = 0 - - self.cell_voltage = V - self.voltage_history = np.append(self.voltage_history, [V]) + V_cell = 0 + self.cell_voltage = V_cell + self.voltage_history = np.append(self.voltage_history, [V_cell]) + self.degradation_history = np.append( + self.degradation_history, [self.V_degradation] + ) + self.power_input_history = np.append(self.power_input_history, [P_in]) # check if it is an hour to decide whether to calculate fatigue hourly_temp = self.hourly_counter self.time += self.dt @@ -290,25 +385,37 @@ def create_polarization(self): return fitobj def convert_power_to_current(self, Pdc, T): - """ - Pdc [kW]: stack power - T [degC]: stack temperature - return :: Idc [A]: stack current + """Estimate stack current for a given power operating point and temperature. + + Args: + Pdc (float): stack power in kWdc + T (float): stack temperature in Celsius + + Returns: + float: ``Idc`` stack current in Amps """ Idc = self.electrolyzer_model((Pdc, T), *self.fit_params) return Idc def curtail_power(self, P_in): - """ - P_in [kWdc]: input power - Curtail power if over electrolyzer rating: + """Curtail power if power exceeds stack rating + + Args: + P_in (float): input power in kWdc + + Returns: + float | array: input power in kWdc saturated at stack rated power. """ return np.where(P_in > self.stack_rating_kW, self.stack_rating_kW, P_in) def calc_fatigue_degradation(self, voltage_signal): """ - voltage_signal: the voltage signal from the last 3600 seconds - return:: voltage_penalty: the degradation penalty + Args: + voltage_signal (float | array): the voltage signal from the last + 3600 seconds + Returns: + float | array: ``voltage_penalty`` the degradation penalty from + variable operation in Volts """ # based off degradation due to square waves of different frequencies # from results in https://iopscience.iop.org/article/10.1149/2.0231915jes @@ -372,13 +479,18 @@ def update_temperature(self, I, V): return self.temperature def update_dynamics(self, H2_mfr_ss, stack_state): - """ - H2_mfr_ss: steady state mass flow rate - stack_state: previous mfr state - return :: next_state: next mfr state - return :: H2_mfr_actual: actual mfr according to the filter + """This is really just a filter on the steady state mfr from + time step to time step. + + Args: + H2_mfr_ss (float): steady state mass flow rate + stack_state (float): previous mfr state - This is really just a filter on the steady state mfr from time step to time step + Returns: + 2-element tuple containing + + - **next_state** (float): next mfr state + - **H2_mfr_actual** (float): actual mfr according to the filter """ if self.ignore_dynamics: @@ -406,6 +518,9 @@ def calc_state_space(self): return [ss_d[0], ss_d[1], ss_d[2], ss_d[3]] def update_status(self): + """Update the stack status if the stack is waiting. If the stack is waiting + and has waited long enough to be on, this method updates the stack status to on. + """ # Change the stack to be truly on if it has waited long enough if self.stack_on: return @@ -416,13 +531,15 @@ def update_status(self): self.stack_on = True def turn_stack_off(self): + """Turn the stack off if the stack is on or watiting. + Updates the cycle count, waiting period, turn off time, and stack status. + """ if self.stack_on or self.stack_waiting: # record turn off time to adjust waiting period self.turn_off_time = self.time self.stack_on = False self.stack_waiting = False self.cycle_count += 1 - # self.stack_state = 0 # adjust waiting period self.wait_time = np.max( @@ -430,6 +547,9 @@ def turn_stack_off(self): ) def turn_stack_on(self): + """Turn the stack on if the stack is off or watiting. + Updates the waiting period, turn on time, and stack status. + """ if self.stack_on: return @@ -450,9 +570,11 @@ def turn_stack_on(self): def calc_stack_power(self, Idc, V=None): """ Args: - Idc [A]: stack current - V (optional): stack voltage - return :: Pdc [kW]: stack power + Idc (float): stack current in Amps + V (float, optional): stack voltage + + Returns: + float: ``Pdc`` [kW] stack power """ V = V or (self.cell.calc_cell_voltage(Idc, self.temperature)) Pdc = Idc * V * self.n_cells @@ -461,13 +583,142 @@ def calc_stack_power(self, Idc, V=None): return Pdc def calc_electrolysis_efficiency(self, Pdc, mfr): - """ - Pdc [kW]: stack power - mfr [kg/h]: mass flow rate - return :: eta_kWh_per_kg, eta_hhv_percent, and eta_lhv_percent: efficiencies + """Calculate the efficiency of the stack in kWh/kg, %-HHV and %-LHV. + + Args: + Pdc (float): stack power in kW + mfr (float): mass flow rate of hydrogen in kg/hr + + Returns: + 3-element tuple containing + + - **eta_kWh_per_kg**: efficiency in kWh/kg + - **eta_hhv_percent**: efficiency as %-HHV + - **eta_lhv_percent**: efficiency as %-LHV """ eta_kWh_per_kg = Pdc / mfr eta_hhv_percent = self.cell.hhv / eta_kWh_per_kg * 100.0 eta_lhv_percent = self.cell.lhv / eta_kWh_per_kg * 100.0 return (eta_kWh_per_kg, eta_hhv_percent, eta_lhv_percent) + + def calc_end_of_life_voltage(self): + """Calculate the end-of-life cell degradation voltage based on the + ``eol_eff_percent_loss`` parameter. + + Returns: + d_eol (float): cell degradation in Volts that indicates end-of-life. + """ + + # efficiency drop that indicates end-of-life as a percentage + eol_eff_mult = (100 + self.eol_eff_percent_loss) / 100 + V_cell_bol = self.cell.calc_cell_voltage(self.max_current, self.temperature) + H2_mfr_bol = ( + self.cell.calc_mass_flow_rate(self.temperature, self.max_current) + * self.n_cells + ) + + H2_mfr_eol = H2_mfr_bol / eol_eff_mult + i_eol_no_faradaic_loss = (H2_mfr_eol * 1e3 * self.cell.n * F) / ( + 1 * self.n_cells * self.cell.M * self.dt + ) + n_f = self.cell.calc_faradaic_efficiency( + self.temperature, i_eol_no_faradaic_loss + ) + i_eol = (H2_mfr_eol * 1e3 * self.cell.n * F) / ( + n_f * self.n_cells * self.cell.M * self.dt + ) + + d_eol = ((self.max_current * V_cell_bol) / i_eol) - V_cell_bol + return d_eol + + def estimate_time_until_replacement(self): + """Estimate the time until replacement based on fraction of life used, + which is the ratio of the stack degradation to the end of life degradation. + + Returns: + float: Number of hours until stack should be replaced with respect to + simulation duration (alternatively, with respect to the number of + hours the stack has existed in the plant.) + """ + + frac_of_life_used = self.V_degradation / self.d_eol + # time between replacement [hrs] based on time its existed (whether on or off) + time_between_replacement = (1 / frac_of_life_used) * (self.time / 3600) # [hrs] + return time_between_replacement + + def estimate_stack_life(self): + """Estimate the stack life based on fraction of life used, which is the ratio + of the stack degradation to the end of life degradation. + + Returns: + float: Stack life in hours with respect to number of hours the stack has + been operational. + """ + # stack life [hrs] based on time its been operational + frac_of_life_used = self.V_degradation / self.d_eol + stack_life = (1 / frac_of_life_used) * (self.uptime / 3600) # [hrs] + return stack_life + + def estimate_life_performance_from_year(self, plant_life_years: int): + """Estimate future performance of the stack assuming the + same operation from the simulation for the duration of the plant life. + + Note: + This function is not tested and may only work for simulations + with an hourly timestep and a simulation length of 8760 hours. + + Args: + plant_life_years (int): number of years in the plant life. + + Returns: + 3-element tuple containing + + - **refurb_schedule** (list): refurbishment schedule of the stack, + a value of 1 represents a year that the stack has to be replaced. + - **ahp_kg** (list): annual hydrogen production of the stack per year + of the ``plant_life_years``. Each element has units of kg-H2/year + - **aep_kWh** (list): annual energy consumption of the stack per year + of the ``plant_life_years``. Each element has units of kWh/year. + """ + refurb_schedule = np.zeros(plant_life_years) + ahp_kg = np.zeros(plant_life_years) + aep_kWh = np.zeros(plant_life_years) + + V_deg = np.array(self.degradation_history) + V_cell_bol = np.array(self.voltage_history) - V_deg + sim_length = len(V_cell_bol) + + Vdeg0 = 0 + for i in range(plant_life_years): + V_deg_pr_sim = Vdeg0 + V_deg + if np.max(V_deg_pr_sim) > self.d_eol: + idx_dead = np.argwhere(V_deg_pr_sim > self.d_eol)[0][0] + V_deg_pr_sim = np.concatenate( + [V_deg_pr_sim[0:idx_dead], V_deg[idx_dead:sim_length]] + ) + refurb_schedule[i] = 1 + if self.hydrogen_degradation_penalty: + I_nom = self.electrolyzer_model( + (self.power_input_history / 1e3, self.temperature), *self.fit_params + ) + V_cell = self.cell.calc_cell_voltage(I_nom, self.temperature) + eff_mult = (V_cell + self.V_degradation) / V_cell # 1 + efficiency drop + I_stack = I_nom / eff_mult + + else: + I_stack = self.electrolyzer_model( + (self.power_input_history / 1e3, self.temperature), *self.fit_params + ) + V_cell = self.cell.calc_cell_voltage(I_stack, self.temperature) + + H2_mass_out = ( + self.cell.calc_mass_flow_rate(self.temperature, I_stack) + * self.n_cells + * self.dt + ) + power_usage_kW = self.calc_stack_power(I_stack, V_cell + V_deg_pr_sim) + ahp_kg[i] = np.sum(H2_mass_out) + aep_kWh[i] = np.sum(power_usage_kW) + Vdeg0 = V_deg_pr_sim[sim_length - 1] + return refurb_schedule, ahp_kg, aep_kWh diff --git a/electrolyzer/tools/modeling_schema.yaml b/electrolyzer/tools/modeling_schema.yaml index 2617121..8e5e16c 100644 --- a/electrolyzer/tools/modeling_schema.yaml +++ b/electrolyzer/tools/modeling_schema.yaml @@ -131,10 +131,15 @@ properties: type: boolean default: True description: Toggle whether degradation is applied to the Stack operation + hydrogen_degradation_penalty: + type: boolean + default: True + description: Toggle whether degradation is applied to hydrogen (True) or power (False) degradation: type: object default: + eol_eff_percent_loss: 10 PEM_params: rate_steady: 1.41737929e-10 rate_fatigue: 3.33330244e-07 @@ -145,6 +150,12 @@ properties: rate_onoff: 3.0726072607260716e-04 description: Degradation rates for PEM or alkaline electrolyzer properties: + eol_eff_percent_loss: + type: number + description: End of life efficiency drop + minimum: 5 + default: 10 + maximum: 100 PEM_params: type: object default: @@ -189,12 +200,12 @@ properties: cell_params: type: object - default: - cell_type: PEM - max_current_density: 2 - PEM_params: - cell_area: 1000 - turndown_ratio: 0.1 + default: {} + # cell_type: PEM + # max_current_density: 2 + # PEM_params: + # cell_area: 1000 + # turndown_ratio: 0.1 description: Cell parameters for PEM or alkaline cells properties: cell_type: @@ -203,14 +214,16 @@ properties: description: Cell electrochemistry (PEM, alkaline, or other) PEM_params: type: object - default: - cell_area: 1001 - turndown_ratio: 0.1 - max_current_density: 2 + default: {} + # cell_area: 1001 + # turndown_ratio: 0.1 + # max_current_density: 2 description: PEM cell parameters properties: cell_area: type: number + minimum: 1 + maximum: 10000 default: 1000 description: Area of individual Cells in the Stack (cm^2) unit: cm^2 @@ -222,9 +235,82 @@ properties: description: Minimum operating power as a fraction of rated power max_current_density: type: number + minimum: 1.5 + maximum: 6 default: 2 units: A/cm^2 description: Maximum current density + p_anode: + type: number + minimum: 1 + maximum: 30 + default: 1.01325 + units: bar + description: anode pressure + p_cathode: + type: number + minimum: 1 + maximum: 30 + default: 30 + units: bar + description: cathode pressure + alpha_a: + type: number + minimum: 0 + maximum: 4 + default: 2 + units: none + description: anode charge transfer coefficient + alpha_c: + type: number + minimum: 0 + maximum: 4 + default: 0.5 + units: none + description: cathode charge transfer coefficient + i_0_a: + type: number + minimum: 0 + maximum: 1 + default: 2.0e-7 + units: A/cm^2 + description: anode exchange current density + i_0_c: + type: number + minimum: 0 + maximum: 1 + default: 2.0e-3 + units: A/cm^2 + description: cathode exchange current density + e_m: + type: number + minimum: 0.002 + maximum: 0.08 + default: 0.02 + units: cm + description: membrane thickness + R_ohmic_elec: + type: number + minimum: 0 + maximum: 1 + default: 50.0e-3 + units: A*cm^2 + description: electrolyte resistance + f_1: + type: number + minimum: 150 + maximum: 400 + default: 250 + units: mA^2/cm^4 + description: faradaic coefficient + f_2: + type: number + minimum: 0.5 + maximum: 1.0 + default: 0.996 + units: none + description: faradaic coefficient + ALK_params: type: object default: {} @@ -249,35 +335,41 @@ properties: type: number units: A/cm^2 description: Maximum current density + f_1: + type: number + minimum: 150 + maximum: 400 + default: 250 + units: mA^2/cm^4 + description: faradaic coefficient + f_2: + type: number + minimum: 0.5 + maximum: 1.0 + default: 0.96 + units: none + description: faradaic coefficient + cell_area: + type: number + # minimum: 1 + # maximum: 10000 + default: 300 + description: Area of individual Cells in the Stack (cm^2) + unit: cm^2 electrode: type: object default: {} description: Alkaline electrode parameters properties: - A_electrode: - type: number - default: 300 - description: electrode total area - units: cm^2 - e_a: + e_e: type: number default: 0.2 - description: anode thickness + description: electrode thickness units: cm - e_c: - type: number - default: 0.2 - description: cathode thickness - units: cm - d_am: + d_em: type: number default: 0.125 - description: anode-membrane distance - units: cm - d_cm: - type: number - default: 0.125 - description: cathode-membrane distance + description: electrode-membrane distance units: cm d_ac: type: number @@ -301,11 +393,6 @@ properties: default: {} description: Alkaline membrane parameters properties: - A_membrane: - type: number - default: 300 - units: cm^2 - description: membrane surface area e_m: type: number default: 0.05 diff --git a/electrolyzer/tools/validators.py b/electrolyzer/tools/validators.py new file mode 100644 index 0000000..543bdd1 --- /dev/null +++ b/electrolyzer/tools/validators.py @@ -0,0 +1,30 @@ +""" +This module contains validator functions for use with `attrs` class definitions. +""" + + +def gt_zero(instance, attribute, value): + """Validates that an attribute's value is greater than zero.""" + if value <= 0: + raise ValueError(f"{attribute} must be greater than zero") + + +def range_val(min, max): + """Validates that an attribute's value is between two values, + inclusive ([min, max]).""" + + def validator(instance, attribute, value): + if value < min or value > max: + raise ValueError(f"{attribute} must be in range [{min}, {max}]") + + return validator + + +def contains(items): + """Validates that an item is part of a given list.""" + + def validator(instance, attribute, value): + if value not in items: + raise ValueError(f"Item {value} not found in list for {attribute}: {items}") + + return validator diff --git a/tests/glue_code/test_run_electrolyzer.py b/tests/glue_code/test_run_electrolyzer.py index 18743d9..3555304 100644 --- a/tests/glue_code/test_run_electrolyzer.py +++ b/tests/glue_code/test_run_electrolyzer.py @@ -135,7 +135,7 @@ def test_regression(result): _, df = result # Test total kg H2 produced - assert_almost_equal(df["kg_rate"].sum(), 228.0749041980651, decimal=1) + assert_almost_equal(df["kg_rate"].sum(), 243.93038801007688, decimal=1) # Test degradation state of stacks degradation = df[[col for col in df if "deg" in col]] diff --git a/tests/glue_code/test_run_lcoh.py b/tests/glue_code/test_run_lcoh.py index c4970af..b453e24 100644 --- a/tests/glue_code/test_run_lcoh.py +++ b/tests/glue_code/test_run_lcoh.py @@ -11,18 +11,23 @@ lcoh_breakdown = pd.DataFrame( { - "Life Totals [$]": [5.388657e06, 1.079412e06, 1.1979687e07, 1.225039e06], + "Life Totals [$]": [ + 5388657.433992826, + 1079412.3726892117, + 11981942.92917099, + 1225039.714867523, + ], "Life Totals [$/kg-H2]": [ - 1.3285182681325287, - 0.2661180588173578, - 2.9534690901554037, - 0.3020209876624967, + 1.242164580703914, + 0.24882038499422945, + 2.7620135992953014, + 0.2823896234644343, ], }, index=["CapEx", "OM", "Feedstock", "Stack Rep"], ) -RESULT = (lcoh_breakdown, 4.850126404767788) +RESULT = (lcoh_breakdown, 4.535388188457879) ROOT = Path(__file__).parent.parent.parent @@ -37,11 +42,11 @@ def test_run_lcoh(): test_signal_angle = np.linspace(0, 8 * np.pi, 3600 * 8 + 10) base_value = (turbine_rating / 2) + 0.2 variation_value = turbine_rating - base_value - power_test_signal = (base_value + variation_value * np.cos(test_signal_angle)) * 1e6 + power_signal = (base_value + variation_value * np.cos(test_signal_angle)) * 1e6 lcoe = 44.18 * (1 / 1000) - calc_result = run_lcoh(fname_input_modeling, power_test_signal, lcoe) + calc_result = run_lcoh(fname_input_modeling, power_signal, lcoe) assert_frame_equal( calc_result[0]["LCOH Breakdown"], @@ -65,12 +70,12 @@ def test_run_lcoh_opt(): test_signal_angle = np.linspace(0, 8 * np.pi, 3600 * 8 + 10) base_value = (turbine_rating / 2) + 0.2 variation_value = turbine_rating - base_value - power_test_signal = (base_value + variation_value * np.cos(test_signal_angle)) * 1e6 + power_signal = (base_value + variation_value * np.cos(test_signal_angle)) * 1e6 lcoe = 44.18 * (1 / 1000) - h2_result = run_electrolyzer(fname_input_modeling, power_test_signal, optimize=True) - lcoh_result = run_lcoh(fname_input_modeling, power_test_signal, lcoe, optimize=True) + h2_result = run_electrolyzer(fname_input_modeling, power_signal, optimize=True) + lcoh_result = run_lcoh(fname_input_modeling, power_signal, lcoe, optimize=True) # h2 prod, max curr density, LCOH assert len(lcoh_result) == 3 @@ -78,4 +83,4 @@ def test_run_lcoh_opt(): # results from regular optimize run should match assert_array_almost_equal(h2_result, lcoh_result[:2]) - assert np.isclose(lcoh_result[2], 4.7541633856347625) + assert np.isclose(lcoh_result[2], 4.44566272289819) diff --git a/tests/test_PEM_cell.py b/tests/test_PEM_cell.py index 53491bc..e114492 100644 --- a/tests/test_PEM_cell.py +++ b/tests/test_PEM_cell.py @@ -1,13 +1,210 @@ -"""This module provides unit tests for `PEMCell`""" +"""This module provides unit tests for `Cell`.""" import pytest +from numpy.testing import assert_almost_equal from electrolyzer.simulation.cell_models.pem import PEMCell as Cell -# from numpy.testing import assert_almost_equal - - @pytest.fixture def cell(): - return PEMCell.from_dict() + return Cell.from_dict( + { + "cell_area": 1000, + "turndown_ratio": 0.1, + "max_current_density": 2, + "p_anode": 1.01325, + "p_cathode": 30, + "alpha_a": 2, + "alpha_c": 0.5, + "i_0_a": 2.0e-7, + "i_0_c": 2.0e-3, + "e_m": 0.02, + "R_ohmic_elec": 50.0e-3, + "f_1": 250, + "f_2": 0.996, + } + ) + + +def test_init(): + """`Cell` should initialize properly from a Dictionary.""" + cell = Cell.from_dict( + { + "cell_area": 1000, + "turndown_ratio": 0.1, + "max_current_density": 2, + "p_anode": 1.01325, + "p_cathode": 30, + "alpha_a": 2, + "alpha_c": 0.5, + "i_0_a": 2.0e-7, + "i_0_c": 2.0e-3, + "e_m": 0.02, + "R_ohmic_elec": 50.0e-3, + "f_1": 250, + "f_2": 0.996, + } + ) + + assert cell.cell_area == 1000 + assert cell.n == 2 + assert cell.gibbs == 237.24e3 + assert cell.M == 2.016 # molecular weight [g/mol] + assert cell.lhv == 33.33 # lower heating value of H2 [kWh/kg] + assert cell.hhv == 39.41 # higher heating value of H2 [kWh/kg] + # assert cell.p_anode ==1 + assert cell.p_anode == 1.01325 + assert cell.p_cathode == 30 + assert cell.alpha_a == 2 + assert cell.alpha_c == 0.5 + assert cell.i_0_a == 2.0e-7 + assert cell.i_0_c == 2.0e-3 + assert cell.e_m == 0.02 + assert cell.R_ohmic_elec == 50.0e-3 + assert cell.f_1 == 250 + assert cell.f_2 == 0.996 + + +def test_calc_reversible_voltage(cell: Cell): + """Reversible cell potential should match literature.""" + E_rev = cell.calc_reversible_voltage() + + assert_almost_equal(E_rev, 1.229, decimal=3) + + +def test_calc_open_circuit_voltage(cell: Cell): + """Open circuit voltage should follow Nernst equation.""" + T = 60 # C + E_rev = cell.calc_reversible_voltage() + + # should be greater than reversible cell voltage + E_cell = cell.calc_open_circuit_voltage(T) + assert E_cell > E_rev + + # should approach E_rev at near 100C (valid temperature range) + E_cell_25 = cell.calc_open_circuit_voltage(99.9725) + assert_almost_equal(E_cell_25, E_rev, decimal=3) + + +def test_calc_activation_overpotential(cell: Cell): + """Activation overpotential should follow Butler-Volmer equations.""" + T = 60 # C + I = 2000 # current + i = I / cell.cell_area # current density + V_act = cell.calc_activation_overpotential(i, T) + V_act_a, V_act_c = V_act + + # should be gte 0 + assert sum(V_act) >= 0 + + # cathode should have a higher overpotential, under assumption that reaction + # kinetics are faster at the anode (oxidation) than cathode (reduction) + assert V_act_c > V_act_a + + # should increase with temperature + T2 = 80 + V_act_T2 = cell.calc_activation_overpotential(i, T2) + assert sum(V_act_T2) > sum(V_act) + + # should increase with current density + cell.cell_area /= 2 + i2 = I / cell.cell_area # current density increases + V_act_i2 = cell.calc_activation_overpotential(i2, T) + assert sum(V_act_i2) > sum(V_act) + + +def test_calc_ohmic_overpotential(cell: Cell): + """Ohmic overpotential should reflect standard Ohm's law.""" + T = 60 # C + I = 2000 # current + i = I / cell.cell_area # current density + + V_ohm = cell.calc_ohmic_overpotential(i, T) + + # should be gte 0 + assert V_ohm >= 0 + + # should decrease with temperature due to increased membrane conductivity + T2 = 80 + V_ohm_T2 = cell.calc_ohmic_overpotential(i, T2) + assert V_ohm_T2 < V_ohm + + # should increase with current density (V = IR) + cell.cell_area /= 2 + i2 = I / cell.cell_area + V_ohm_i2 = cell.calc_ohmic_overpotential(i2, T) + assert V_ohm_i2 > V_ohm + + +def test_calc_concentration_overpotential(cell: Cell): + """TODO: This method is not yet implemented on the class.""" + assert cell.calc_concentration_overpotential() == 0 + + +def test_calc_overpotentials(cell: Cell): + """Should calculate and return all overpotentials.""" + T = 60 # C + I = 2000 # current + i = I / cell.cell_area + overpotentials = cell.calc_overpotentials(i, T) + + # Activation (cathode/anode), Ohmic, Concentration + assert len(overpotentials) == 4 + + # Should all be gte zero + assert all([v >= 0 for v in overpotentials]) + + +def test_calc_cell_voltage(cell: Cell): + """Should calculate cell voltage""" + T = 60 # C + I = 2000 # current + E_rev = cell.calc_reversible_voltage() + + V_cell = cell.calc_cell_voltage(I, T) + + # should be higher than reversible cell voltage due to overpotentials + assert V_cell > E_rev + + # should increase with current + I2 = 2500 + V_cell_I2 = cell.calc_cell_voltage(I2, T) + assert V_cell_I2 > V_cell + + # should decrease with temperature + T2 = 100 + V_cell_T2 = cell.calc_cell_voltage(I, T2) + assert V_cell_T2 < V_cell + + +def test_calc_faradaic_efficiency(cell: Cell): + """Should calculate Faraday's efficiency.""" + I = 500 + + # should increase with current + eta = cell.calc_faradaic_efficiency(60, I) + + I2 = 2000 + eta2 = cell.calc_faradaic_efficiency(60, I2) + + assert eta2 > eta + + # should approach 1 as current increases + I3 = 5000 + eta3 = cell.calc_faradaic_efficiency(60, I3) + + assert_almost_equal(eta3, 0.996, decimal=3) + + +def test_calc_mass_flow_rate(cell: Cell): + """Should calculate the mass flow rate of H2.""" + I = 1000 + + mfr = cell.calc_mass_flow_rate(60, I) + + # should increase with current + I2 = 2000 + mfr2 = cell.calc_mass_flow_rate(60, I2) + + assert mfr2 > mfr diff --git a/tests/test_cell.py b/tests/test_cell.py deleted file mode 100644 index e504631..0000000 --- a/tests/test_cell.py +++ /dev/null @@ -1,171 +0,0 @@ -"""This module provides unit tests for `Cell`.""" - -import pytest -from numpy.testing import assert_almost_equal - -from electrolyzer.simulation.cell_models.pem import PEMCell as Cell - - -@pytest.fixture -def cell(): - return Cell.from_dict( - {"cell_area": 1000, "turndown_ratio": 0.1, "max_current_density": 2} - ) - - -def test_init(): - """`Cell` should initialize properly from a Dictionary.""" - cell = Cell.from_dict( - {"cell_area": 1000, "turndown_ratio": 0.1, "max_current_density": 2} - ) - - assert cell.cell_area == 1000 - assert cell.n == 2 - assert cell.gibbs == 237.24e3 - assert cell.M == 2.016 # molecular weight [g/mol] - assert cell.lhv == 33.33 # lower heating value of H2 [kWh/kg] - assert cell.hhv == 39.41 # higher heating value of H2 [kWh/kg] - - -def test_calc_reversible_voltage(cell: Cell): - """Reversible cell potential should match literature.""" - E_rev = cell.calc_reversible_voltage() - - assert_almost_equal(E_rev, 1.229, decimal=3) - - -def test_calc_open_circuit_voltage(cell: Cell): - """Open circuit voltage should follow Nernst equation.""" - T = 60 # C - E_rev = cell.calc_reversible_voltage() - - # should be greater than reversible cell voltage - E_cell = cell.calc_open_circuit_voltage(T) - assert E_cell > E_rev - - # should approach E_rev at near 100C (valid temperature range) - E_cell_25 = cell.calc_open_circuit_voltage(99.9725) - assert_almost_equal(E_cell_25, E_rev, decimal=3) - - -def test_calc_activation_overpotential(cell: Cell): - """Activation overpotential should follow Butler-Volmer equations.""" - T = 60 # C - I = 2000 # current - i = I / cell.cell_area # current density - V_act = cell.calc_activation_overpotential(i, T) - V_act_a, V_act_c = V_act - - # should be gte 0 - assert sum(V_act) >= 0 - - # cathode should have a higher overpotential, under assumption that reaction - # kinetics are faster at the anode (oxidation) than cathode (reduction) - assert V_act_c > V_act_a - - # should increase with temperature - T2 = 80 - V_act_T2 = cell.calc_activation_overpotential(i, T2) - assert sum(V_act_T2) > sum(V_act) - - # should increase with current density - cell.cell_area /= 2 - i2 = I / cell.cell_area # current density increases - V_act_i2 = cell.calc_activation_overpotential(i2, T) - assert sum(V_act_i2) > sum(V_act) - - -def test_calc_ohmic_overpotential(cell: Cell): - """Ohmic overpotential should reflect standard Ohm's law.""" - T = 60 # C - I = 2000 # current - i = I / cell.cell_area # current density - - V_ohm = cell.calc_ohmic_overpotential(i, T) - - # should be gte 0 - assert V_ohm >= 0 - - # should decrease with temperature due to increased membrane conductivity - T2 = 80 - V_ohm_T2 = cell.calc_ohmic_overpotential(i, T2) - assert V_ohm_T2 < V_ohm - - # should increase with current density (V = IR) - cell.cell_area /= 2 - i2 = I / cell.cell_area - V_ohm_i2 = cell.calc_ohmic_overpotential(i2, T) - assert V_ohm_i2 > V_ohm - - -def test_calc_concentration_overpotential(cell: Cell): - """TODO: This method is not yet implemented on the class.""" - assert cell.calc_concentration_overpotential() == 0 - - -def test_calc_overpotentials(cell: Cell): - """Should calculate and return all overpotentials.""" - T = 60 # C - I = 2000 # current - i = I / cell.cell_area - overpotentials = cell.calc_overpotentials(i, T) - - # Activation (cathode/anode), Ohmic, Concentration - assert len(overpotentials) == 4 - - # Should all be gte zero - assert all([v >= 0 for v in overpotentials]) - - -def test_calc_cell_voltage(cell: Cell): - """Should calculate cell voltage""" - T = 60 # C - I = 2000 # current - E_rev = cell.calc_reversible_voltage() - - V_cell = cell.calc_cell_voltage(I, T) - - # should be higher than reversible cell voltage due to overpotentials - assert V_cell > E_rev - - # should increase with current - I2 = 2500 - V_cell_I2 = cell.calc_cell_voltage(I2, T) - assert V_cell_I2 > V_cell - - # should decrease with temperature - T2 = 100 - V_cell_T2 = cell.calc_cell_voltage(I, T2) - assert V_cell_T2 < V_cell - - -def test_calc_faradaic_efficiency(cell: Cell): - """Should calculate Faraday's efficiency.""" - I = 500 - - # should increase with current - eta = cell.calc_faradaic_efficiency(60, I) - - I2 = 2000 - eta2 = cell.calc_faradaic_efficiency(60, I2) - - assert eta2 > eta - - # should approach 1 as current increases - I3 = 5000 - eta3 = cell.calc_faradaic_efficiency(60, I3) - - assert_almost_equal(eta3, 0.996, decimal=3) - - -def test_calc_mass_flow_rate(cell: Cell): - """Should calculate the mass flow rate of H2.""" - I = 1000 - - mfr = cell.calc_mass_flow_rate(60, I) - - # should increase with current - I2 = 2000 - mfr2 = cell.calc_mass_flow_rate(60, I2) - - assert mfr2 > mfr diff --git a/tests/test_stack.py b/tests/test_stack.py index 80e40ef..a45a809 100644 --- a/tests/test_stack.py +++ b/tests/test_stack.py @@ -17,11 +17,12 @@ def create_stack(): "n_cells": 100, # "stack_rating_kW": 750, "degradation": { + "eol_eff_percent_loss": 10, "PEM_params": { "rate_steady": 1.41737929e-10, "rate_fatigue": 3.33330244e-07, "rate_onoff": 1.47821515e-04, - } + }, }, "cell_params": { "cell_type": "PEM", @@ -29,6 +30,16 @@ def create_stack(): "cell_area": 1000, "turndown_ratio": 0.1, "max_current_density": 2, + "p_anode": 1.01325, + "p_cathode": 30, + "alpha_a": 2, + "alpha_c": 0.5, + "i_0_a": 2.0e-7, + "i_0_c": 2.0e-3, + "e_m": 0.02, + "R_ohmic_elec": 50.0e-3, + "f_1": 250, + "f_2": 0.996, }, }, } @@ -57,11 +68,12 @@ def test_init(mocker): "n_cells": 100, "stack_rating_kW": 750, "degradation": { + "eol_eff_percent_loss": 10, "PEM_params": { "rate_steady": 1.41737929e-10, "rate_fatigue": 3.33330244e-07, "rate_onoff": 1.47821515e-04, - } + }, }, "cell_params": { "cell_type": "PEM", @@ -69,8 +81,19 @@ def test_init(mocker): "cell_area": 1000, "turndown_ratio": 0.1, "max_current_density": 2, + "p_anode": 1.01325, + "p_cathode": 30, + "alpha_a": 2, + "alpha_c": 0.5, + "i_0_a": 2.0e-7, + "i_0_c": 2.0e-3, + "e_m": 0.02, + "R_ohmic_elec": 50.0e-3, + "f_1": 250, + "f_2": 0.996, }, }, + "hydrogen_degradation_penalty": False, } stack = Stack.from_dict(stack_dict) @@ -86,6 +109,9 @@ def test_init(mocker): assert stack.min_power == 0.1 * stack.stack_rating assert stack.include_degradation_penalty is True + assert ( + stack.hydrogen_degradation_penalty == stack_dict["hydrogen_degradation_penalty"] + ) assert stack.rf_track == 0.0 assert stack.V_degradation == 0.0 assert stack.uptime == 0.0 @@ -164,6 +190,7 @@ def test_run(mocker): def test_create_polarization(stack: Stack): + # TODO remake this """ Should create a polarization curve based on fit for the specified model over a range of temperatures. @@ -411,6 +438,7 @@ def test_calc_stack_power(stack: Stack): def test_calc_electrolysis_efficiency(stack: Stack): + # TODO: remake this test - something is weird about it """ Should calculate values of electrolysis efficiency for given DC Power input and MFR. """ @@ -425,7 +453,7 @@ def test_calc_electrolysis_efficiency(stack: Stack): assert len(eta_values) == 3 # efficiency should decrease as we approach max current due to overpotentials - assert eta_values[0] > 80 # highest efficiency around 80% capacity + assert eta_values[0] > 70 # highest efficiency around 80% capacity H2_mfr2 = ( stack.cell.calc_mass_flow_rate(stack.temperature, stack.max_current) * stack.n_cells @@ -445,11 +473,12 @@ def test_dt_behavior(): "n_cells": 100, # "stack_rating_kW": 750, "degradation": { + "eol_eff_percent_loss": 10, "PEM_params": { "rate_steady": 1.41737929e-10, "rate_fatigue": 3.33330244e-07, "rate_onoff": 1.47821515e-04, - } + }, }, "cell_params": { "cell_type": "PEM", @@ -457,6 +486,16 @@ def test_dt_behavior(): "cell_area": 1000, "turndown_ratio": 0.1, "max_current_density": 2, + "p_anode": 1.01325, + "p_cathode": 30, + "alpha_a": 2, + "alpha_c": 0.5, + "i_0_a": 2.0e-7, + "i_0_c": 2.0e-3, + "e_m": 0.02, + "R_ohmic_elec": 50.0e-3, + "f_1": 250, + "f_2": 0.996, }, }, }