From 1919d309534c274b5682a80e3b86c9af6123c994 Mon Sep 17 00:00:00 2001 From: elenya-grant <116225007+elenya-grant@users.noreply.github.com> Date: Sat, 6 Jun 2026 14:14:15 -0600 Subject: [PATCH 1/7] added functions to pem cell --- electrolyzer/components/cell/cell.py | 98 +++++++- electrolyzer/components/cell/pem_cell.py | 278 ++++++++++++++++++++++- 2 files changed, 362 insertions(+), 14 deletions(-) diff --git a/electrolyzer/components/cell/cell.py b/electrolyzer/components/cell/cell.py index e568084..f48dbcc 100644 --- a/electrolyzer/components/cell/cell.py +++ b/electrolyzer/components/cell/cell.py @@ -1,7 +1,14 @@ +import numpy as np import openmdao.api as om from attrs import field, define, validators +from openmdao.utils import units +from scipy.constants import physical_constants from electrolyzer.core.utilities import BaseConfig +from electrolyzer.components.constants import H2_MW, O2_MW + + +F, _, _ = physical_constants["Faraday constant"] # Faraday's constant [C/mol] @define(kw_only=True) @@ -32,25 +39,32 @@ def setup(self): self.add_input("R_elec", self.config.R_ohmic_elec, units="ohm*(cm**2)") # input profiles + # TODO: make either current density or current self.add_input("current_in", val=0.0, shape_by_conn=True, units="A") # OR current density? - self.add_input("current_density_in", val=0.0, shape_by_conn=True, units="A/(cm**2)") + # self.add_input("current_density_in", val=0.0, shape_by_conn=True, units="A/(cm**2)") # output profiles self.add_output("cell_voltage", val=0.0, copy_shape="current_in", units="V") self.add_output( - "hydrogen_produced", val=0.0, copy_shape="current_in", units=f"kg/({self.dt}*s)" + "hydrogen_produced", val=0.0, copy_shape="current_in", units=f"g/({self.dt}*s)" ) self.add_output( - "oxygen_produced", val=0.0, copy_shape="current_in", units=f"kg/({self.dt}*s)" + "oxygen_produced", val=0.0, copy_shape="current_in", units=f"g/({self.dt}*s)" ) self.add_output( - "water_consumed", val=0.0, copy_shape="current_in", units=f"kg/({self.dt}*s)" + "water_consumed", val=0.0, copy_shape="current_in", units=f"g/({self.dt}*s)" ) - self.add_output("hydrogen_production_rate", val=0.0, copy_shape="current_in", units="kg/s") - self.add_input("current_density_out", val=0.0, copy_shape="current_in", units="A/(cm**2)") + self.add_output("hydrogen_production_rate", val=0.0, copy_shape="current_in", units="g/s") + self.add_output("oxygen_production_rate", val=0.0, copy_shape="current_in", units="g/s") + + self.add_output("current_density_out", val=0.0, copy_shape="current_in", units="A/(cm**2)") # output design variables self.add_output("rated_cell_voltage", val=0.0, units="V") + self.add_output("rated_conversion_efficiency", val=0.0, units="W*h/kg") + self.add_output("rated_hydrogen_production", val=0.0, units="g/s") + self.add_output("rated_oxygen_production", val=0.0, units="g/s") + self.add_output("rated_cell_power", val=0.0, units="W") def compute(self, inputs, outputs, discrete_inputs, discrete_outputs): """ @@ -60,3 +74,75 @@ def compute(self, inputs, outputs, discrete_inputs, discrete_outputs): """ raise NotImplementedError("This method should be implemented in a subclass.") + + def membrane_conductivity(self, temperature_celsius): + temp_K = units.convert_units(temperature_celsius, "degC", "K") + lambda_water_content = ((-2.89556 + (0.016 * temp_K)) + 1.625) / 0.1875 + + # membrane proton conductivity [S/cm] + sigma = ((0.005139 * lambda_water_content) - 0.00326) * np.exp( + 1268 * ((1 / 303) - (1 / temp_K)) + ) + return sigma + + def water_vaporization_pressure_antoine(self, temperature_celsius): + # coefficient for Antoine formulas + if temperature_celsius <= 100: + A = 8.07131 + B = 1730.63 + C = 233.426 + elif temperature_celsius > 100 and temperature_celsius < 374: + A = 8.14019 + B = 1810.94 + C = 224.485 + + # vapor pressure of water in [mmHg] using Antoine formula + # valid for T<283 Kelvin + p_h2o_sat_mmHg = 10 ** (A - (B / (C + temperature_celsius))) + + # convert mmHg to bar + p_h2O_sat = units.convert_units(p_h2o_sat_mmHg, "torr", "bar") + return p_h2O_sat + + def water_vaporization_pressure_arden_buck(self, temperature_celsius): + # Arden Buck equation T=C, https://www.omnicalculator.com/chemistry/vapour-pressure-of-water#vapor-pressure-formulas # noqa + # Reasonable at temperatures between 0-100C + # different coefficients if temperature < 0 degrees C + p_h2o_sat_Pa = ( + 0.61121 + * np.exp( + (18.678 - (temperature_celsius / 234.5)) + * (temperature_celsius / (257.14 + temperature_celsius)) + ) + ) * 1e3 # (Pa) + + # convert Pa to bar + p_h2o_sat = units.convert_units(p_h2o_sat_Pa, "Pa", "bar") + return p_h2o_sat + + def calculate_faradaic_efficiency(self, J_cell, f1, f2): + J_cell_mA_pr_cm2 = J_cell * 1e3 + + n_F = ((J_cell_mA_pr_cm2**2) / (f1 + (J_cell_mA_pr_cm2**2))) * f2 + + return n_F + + def calculate_h2_production_rate(self, J_cell, I_cell, f1, f2): + n_F = self.calculate_faradaic_efficiency(J_cell, f1, f2) + h2_production_mol_per_s = n_F * (I_cell / (2 * F)) + h2_production_g_per_s = h2_production_mol_per_s * H2_MW + return h2_production_g_per_s + + def calculate_o2_production_rate(self, J_cell, I_cell, f1, f2): + n_F = self.calculate_faradaic_efficiency(J_cell, f1, f2) + o2_production_mol_per_s = n_F * (I_cell / (4 * F)) + o2_production_g_per_s = o2_production_mol_per_s * O2_MW + return o2_production_g_per_s + + def calculate_h2_production(self, J_cell, I_cell, f1, f2): + h2_g_per_s = self.calculate_h2_production_rate(J_cell, I_cell, f1, f2) + return h2_g_per_s * self.dt + + def calculate_o2_production(self, J_cell, I_cell, f1, f2): + o2_g_per_s = self.calculate_o2_production_rate(J_cell, I_cell, f1, f2) + return o2_g_per_s * self.dt diff --git a/electrolyzer/components/cell/pem_cell.py b/electrolyzer/components/cell/pem_cell.py index 250a29b..4642e42 100644 --- a/electrolyzer/components/cell/pem_cell.py +++ b/electrolyzer/components/cell/pem_cell.py @@ -1,17 +1,27 @@ +import numpy as np from attrs import field, define, validators +from openmdao.utils import units +from scipy.constants import R, physical_constants from electrolyzer.tools.validators import contains, range_val from electrolyzer.components.cell.cell import CellBaseClass, CellBaseConfig +from electrolyzer.components.constants import gibbs + + +F, _, _ = physical_constants["Faraday constant"] # Faraday's constant [C/mol] @define(kw_only=True) class PEMCellConfig(CellBaseConfig): f1: float = field(default=250, validator=validators.ge(0)) f2: float = field(default=0.996, validator=range_val(0.5, 1.0)) - i_0a: float = field(default=2.0e-7, validator=range_val(0.0, 1.0)) - i_0c: float = field(default=2.0e-3, validator=range_val(0.0, 1.0)) - alpha_a: float = field(default=2.0, validator=range_val(0.0, 4.0)) - alpha_c: float = field(default=0.5, validator=range_val(0.0, 4.0)) + i_0a: float = field(default=2.0e-7, validator=validators.optional(range_val(0.0, 1.0))) + i_0c: float = field(default=2.0e-3, validator=validators.optional(range_val(0.0, 1.0))) + alpha_a: float = field(default=2.0, validator=range_val(0.25, 4.0)) + alpha_c: float = field(default=0.5, validator=range_val(0.25, 4.0)) + + b_combined: float = field(default=None, validator=validators.optional(validators.ge(0.0))) + i0_combined: float = field(default=None, validator=validators.optional(validators.ge(0.0))) water_vaporization_method: str = field( default="buck", converter=(str.lower, str.strip), validator=contains(["buck", "antoine"]) @@ -19,19 +29,271 @@ class PEMCellConfig(CellBaseConfig): activation_method: str = field( default="arcsinh", validator=contains(["ln", "log10", "arcsinh"]) ) + kinetics_method: str = field( + default="per_electrode", + converter=(str.lower, str.strip), + validator=contains(["per_electrode", "combined"]), + ) + Urev0_calc_method: str = field( default="normal", validator=contains(["normal", "temp_adjusted"]) ) + def __attrs_post_init__(self): + provided_per_electrode_info = all( + getattr(self, param, None) is not None for param in ["i_0a", "i_0c"] + ) + provided_combined_info = all( + getattr(self, param, None) is not None for param in ["b_combined", "i0_combined"] + ) + + extra_inputs_msg = ( + "Extraneous inputs ({extraneous_attrs_msg}) for kinetics_method {kinetics_method}. " + ) + missing_inputs_msg = ( + "Missing inputs ({required_attrs_msg}) for kinetics_method {kinetics_method}. " + ) + + if self.kinetics_method == "per_electrode": + if not provided_per_electrode_info or provided_combined_info: + # check if theres missing info or extraneous info + required_attributes = ["i_0a", "i_0c"] + extreanous_attributes = ["b_combined", "i0_combined"] + extraneous_attrs_msg = ", ".join( + f"`{k}`" for k in extreanous_attributes if getattr(self, k, None) is not None + ) + required_attrs_msg = ", ".join( + f"`{k}`" for k in required_attributes if getattr(self, k, None) is None + ) + txt = "" + if len(extraneous_attrs_msg) > 0: + txt += extra_inputs_msg.format( + extraneous_attrs_msg=extraneous_attrs_msg, + kinetics_method=self.kinetics_method, + ) + + if len(required_attrs_msg) > 0: + txt += missing_inputs_msg.format( + required_attrs_msg=required_attrs_msg, kinetics_method=self.kinetics_method + ) + if len(txt) > 0: + raise AttributeError(txt) + + if self.kinetics_method == "combined": + if not provided_combined_info or provided_per_electrode_info: + # check if theres missing info or extraneous info + required_attributes = ["b_combined", "i0_combined"] + extreanous_attributes = ["i_0a", "i_0c"] + extraneous_attrs_msg = ", ".join( + f"`{k}`" for k in extreanous_attributes if getattr(self, k, None) is not None + ) + required_attrs_msg = ", ".join( + f"`{k}`" for k in required_attributes if getattr(self, k, None) is None + ) + txt = "" + if len(extraneous_attrs_msg) > 0: + txt += extra_inputs_msg.format( + extraneous_attrs_msg=extraneous_attrs_msg, + kinetics_method=self.kinetics_method, + ) + + if len(required_attrs_msg) > 0: + txt += missing_inputs_msg.format( + required_attrs_msg=required_attrs_msg, kinetics_method=self.kinetics_method + ) + + if len(txt) > 0: + raise AttributeError(txt) + class PEMCell(CellBaseClass): def setup(self): + self.n = 2 # number of electrons transferred in reaction + self.config = PEMCellConfig.from_dict(self.options["tech_config"]["cell_parameters"]) + super().setup() + # Design parameters - self.add_input("i_0a", val=self.config.i_0a, shape=1, units="A/(cm**2)") - self.add_input("i_0c", val=self.config.i_0c, shape=1, units="A/(cm**2)") + if self.config.kinetics_method == "per_electrode": + self.add_input("i_0a", val=self.config.i_0a, shape=1, units="A/(cm**2)") + self.add_input("i_0c", val=self.config.i_0c, shape=1, units="A/(cm**2)") + self.add_input("alpha_a", val=self.config.alpha_a, shape=1, units="unitless") + self.add_input("alpha_c", val=self.config.alpha_c, shape=1, units="unitless") + else: + # b_combined is in V/decade + self.add_input("b_combined", val=self.config.b_combined, shape=1, units="V") + self.add_input("i_0combined", val=self.config.i0_combined, shape=1, units="A/(cm**2)") + self.add_input("f1", val=self.config.f1, shape=1, units="(mA**2)/(cm**4)") self.add_input("f2", val=self.config.f2, shape=1, units="unitless") - self.add_input("alpha_a", val=self.config.alpha_a, shape=1, units="unitless") - self.add_input("alpha_c", val=self.config.alpha_c, shape=1, units="unitless") + + def calculate_current_density(self, I_cell, A_cell): + return I_cell / A_cell + + def calculate_current(self, J_cell, A_cell): + return J_cell * A_cell + + def calc_Urev0(self, temperature_celsius): + E_cell0 = gibbs / (self.n * F) # 1.229 # [V] + + # Reversible potential at 25degC - Nerst Equation + if self.config.Urev0_calc_method == "normal": + return E_cell0 + if self.config.Urev0_calc_method == "temp_adjusted": + temp_K = units.convert_units(temperature_celsius, "degC", "K") + Urev0 = E_cell0 - (0.9 * 1e-3 * (temp_K - 298)) + return Urev0 + + def reversible_overpotential(self, temperature_celsius, p_cathode_bar, p_anode_bar): + if self.config.water_vaporization_method == "buck": + p_H2O_sat = self.water_vaporization_pressure_arden_buck(temperature_celsius) + elif self.config.water_vaporization_method == "antoine": + p_H2O_sat = self.water_vaporization_pressure_antoine(temperature_celsius) + + temp_k = units.convert_units(temperature_celsius, "degC", "K") + + # Daltons law of partial pressures + p_H2 = p_cathode_bar - p_H2O_sat + p_O2 = p_anode_bar - p_H2O_sat + + # Nerst Equation + Urev0 = self.calc_Urev0(temperature_celsius) + + E_cell = Urev0 + ((R * temp_k) / (self.n * F)) * ( + np.log((p_H2 * np.sqrt(p_O2)) / p_H2O_sat) + ) + + return E_cell + + def separate_activation_overpotential( + self, J_cell, temperature_celsius, i_0a, i_0c, alpha_a, alpha_c + ): + temp_K = units.convert_units(temperature_celsius, "degC", "K") + + ba = R * temp_K / (alpha_a * F) + bc = R * temp_K / (alpha_c * F) + + anode_ratio = np.clip(J_cell / i_0a, a_min=1e-30) + cathode_ratio = np.clip(J_cell / i_0c, a_min=1e-30) + + if self.config.activation_method == "arcsinh": + V_acta = ba * np.arcsinh(anode_ratio) + V_actc = bc * np.arcsinh(cathode_ratio) + if self.config.activation_method == "ln": + V_acta = ba * np.log(anode_ratio) + V_actc = bc * np.log(cathode_ratio) + if self.config.activation_method == "log10": + V_acta = ba * np.log(10) * np.log10(anode_ratio) + V_actc = bc * np.log(10) * np.log10(cathode_ratio) + return V_acta + V_actc + + def combined_kinetics_activation_overpotential(self, J_cell, b_combined, i_0combined): + ratio = np.clip(J_cell / i_0combined, a_min=1e-30) + if self.config.activation_method == "ln": + V_act = b_combined * np.log(ratio) + if self.config.activation_method == "log10": + V_act = b_combined * np.log10(ratio) + if self.config.activation_method == "arcsinh": + # TODO: put this error raising in the config + msg = "Cannot run combined kinetics with arcsinh activation method" + raise ValueError(msg) + + return V_act + + def ohmic_overpotential(self, J_cell, temperature_celsius, delta_membrane, R_elec): + sigma_membrane = self.membrane_conductivity(temperature_celsius) + + # ionic resistance [ohms*cm^2] + R_membrane = delta_membrane / sigma_membrane + + R_tot = R_membrane + R_elec + + V_ohmic = J_cell * R_tot + + return V_ohmic + + def cell_voltage(self, inputs): + temp_C = inputs["operating_temperature"] + if "current_density_in" in inputs: + J_cell = inputs["current_density_in"] + elif "current_in" in inputs: + J_cell = self.calculate_current_density( + inputs["current_in"], inputs["cell_active_area"] + ) + + Urev = self.reversible_overpotential( + temp_C, inputs["cathode_pressure"], inputs["anode_pressure"] + ) + U_ohmic = self.ohmic_overpotential( + J_cell, temp_C, inputs["membrane_thickness"], inputs["R_elec"] + ) + + if self.config.kinetics_method == "per_electrode": + V_act = self.separate_activation_overpotential( + J_cell, temp_C, inputs["i_0a"], inputs["i_0c"], inputs["alpha_a"], inputs["alpha_c"] + ) + else: + V_act = self.combined_kinetics_activation_overpotential( + J_cell, inputs["b_combined"], inputs["i_0combined"] + ) + + V_cell = Urev + U_ohmic + V_act + return V_cell + + def get_current_density(self, inputs): + if "current_density_in" in inputs: + return inputs["current_density_in"] + if "current_in" in inputs: + J_cell = self.calculate_current_density( + inputs["current_in"], inputs["cell_active_area"] + ) + return J_cell + + def get_current(self, inputs): + if "current_in" in inputs: + return inputs["current_in"] + if "current_density_in" in inputs: + I_cell = self.calculate_current( + inputs["current_density_in"], inputs["cell_active_area"] + ) + return I_cell + + def h2_production_rate(self, inputs): + I_cell = self.get_current(inputs) + J_cell = self.get_current_density(inputs) + h2_grams_per_sec = self.calculate_h2_production_rate( + J_cell, I_cell, inputs["f1"], inputs["f2"] + ) + return h2_grams_per_sec + + def o2_production_rate(self, inputs): + I_cell = self.get_current(inputs) + J_cell = self.get_current_density(inputs) + o2_grams_per_sec = self.calculate_o2_production_rate( + J_cell, I_cell, inputs["f1"], inputs["f2"] + ) + return o2_grams_per_sec + + def h2_production(self, inputs): + h2_grams_per_sec = self.h2_production_rate(inputs) + return h2_grams_per_sec * self.dt + + def o2_production(self, inputs): + o2_grams_per_sec = self.o2_production_rate(inputs) + return o2_grams_per_sec * self.dt + + def power_consumption_rate(self, inputs): + V_cell = self.cell_voltage(inputs) + I_cell = self.get_current(inputs) + return V_cell * I_cell + + def energy_consumption(self, inputs): + V_cell = self.cell_voltage(inputs) + I_cell = self.get_current(inputs) + return V_cell * I_cell * self.dt + + def conversion_efficiency(self, inputs): + h2_grams_per_sec = self.h2_production_rate(inputs) + power_W_per_sec = self.power_consumption_rate(inputs) + return power_W_per_sec / np.max([1e-30, h2_grams_per_sec]) From ab738960b6c537b40479af99d6ad6328d01cf531 Mon Sep 17 00:00:00 2001 From: elenya-grant <116225007+elenya-grant@users.noreply.github.com> Date: Sat, 6 Jun 2026 14:15:50 -0600 Subject: [PATCH 2/7] added test file for pem cell config --- .../cell/test/test_pem_cell_config.py | 118 ++++++++++++++++++ 1 file changed, 118 insertions(+) create mode 100644 electrolyzer/components/cell/test/test_pem_cell_config.py diff --git a/electrolyzer/components/cell/test/test_pem_cell_config.py b/electrolyzer/components/cell/test/test_pem_cell_config.py new file mode 100644 index 0000000..59c4bd8 --- /dev/null +++ b/electrolyzer/components/cell/test/test_pem_cell_config.py @@ -0,0 +1,118 @@ +import pytest +from pytest import fixture + +from electrolyzer.components.cell.pem_cell import PEMCellConfig + + +@fixture +def pem_cell_config(): + config = { + "A_cell": 1000, + "membrane_thickness": 0.02, + "temperature": 80.0, + "P_anode": 1.0, + "P_cathode": 1.0, + "R_ohmic_elec": 50.0e-3, + "f1": 250, + "f2": 0.996, + "water_vaporization_method": "antoine", # "buck", "antoine" + "activation_method": "ln", # "ln", "log10", or "arcsinh" + # "kinetics_method": "per_electrode", # "per_electrode" or "combined" + "Urev0_calc_method": "normal", # "normal" or "temp_adjusted" + } + return config + + +def test_config_combined_kinetics(pem_cell_config, subtests): + correct_inputs = { + "kinetics_method": "combined", + "b_combined": 0.045, + "i0_combined": 2.0e-8, + "i_0a": None, + "i_0c": None, + } + + config = PEMCellConfig.from_dict(pem_cell_config | correct_inputs) + with subtests.test("Correct initialization"): + assert config.b_combined == correct_inputs["b_combined"] + assert config.i0_combined == correct_inputs["i0_combined"] + + with subtests.test("Other args are None"): + assert config.i_0a is None + assert config.i_0c is None + + # Test error if wrong method + # wrong_method = { + # "kinetics_method": "per_electrode", + # "b_combined": 0.045, + # "i0_combined": 2.0e-8, + # } + + # with subtests.test("incorrect parameters"): + # with pytest.raises(AttributeError) as excinfo: + # config = PEMCellConfig.from_dict(pem_cell_config | wrong_method) + # err = str(excinfo.value) + # assert "For kinetics_method of 'per_electrode', the inputs" in err + + # Test error if missing i0_combined + missing_i0 = { + "kinetics_method": "combined", + "b_combined": 0.045, + } + with subtests.test("Missing input i0_combined"): + with pytest.raises(AttributeError) as excinfo: + config = PEMCellConfig.from_dict(pem_cell_config | missing_i0) + err = str(excinfo.value) + assert "Missing inputs (`i0_combined`)" in err + + # Test error if missing b_combined + missing_b = { + "kinetics_method": "combined", + "i0_combined": 2.0e-8, + } + + with subtests.test("Missing input b_combined"): + with pytest.raises(AttributeError) as excinfo: + config = PEMCellConfig.from_dict(pem_cell_config | missing_b) + err = str(excinfo.value) + assert "Missing inputs (`b_combined`)" in err + + # Test error if missing both inputs + missing_both = { + "kinetics_method": "combined", + } + + with subtests.test("Missing both inputs"): + with pytest.raises(AttributeError) as excinfo: + config = PEMCellConfig.from_dict(pem_cell_config | missing_both) + err = str(excinfo.value) + assert "Missing inputs (`b_combined`, `i0_combined`)" in err + + # Test error if one extra arg given + extra_arg = { + "kinetics_method": "combined", + "b_combined": 0.045, + "i0_combined": 2.0e-8, + "i_0a": 2.0e-7, + } + + with subtests.test("Extraneous input"): + with pytest.raises(AttributeError) as excinfo: + config = PEMCellConfig.from_dict(pem_cell_config | extra_arg) + err = str(excinfo.value) + assert "Extraneous inputs (`i_0a`)" in err + + # Test error if multiple extra args given + extra_args = { + "kinetics_method": "combined", + "b_combined": 0.045, + "i0_combined": 2.0e-8, + "i_0a": 2.0e-7, + "i_0c": 2.0e-3, + } + + with subtests.test("Extraneous inputs"): + with pytest.raises(AttributeError) as excinfo: + config = PEMCellConfig.from_dict(pem_cell_config | extra_args) + err = str(excinfo.value) + assert "Extraneous inputs (`i_0a`, `i_0c`)" in err From 6df839ce1d326299d0e0ded1b181a1d398f860a1 Mon Sep 17 00:00:00 2001 From: elenya-grant <116225007+elenya-grant@users.noreply.github.com> Date: Sat, 6 Jun 2026 14:31:39 -0600 Subject: [PATCH 3/7] fixed pem config attrs_post_init method --- electrolyzer/components/cell/pem_cell.py | 97 +++++++++++------------- 1 file changed, 45 insertions(+), 52 deletions(-) diff --git a/electrolyzer/components/cell/pem_cell.py b/electrolyzer/components/cell/pem_cell.py index 4642e42..71f63db 100644 --- a/electrolyzer/components/cell/pem_cell.py +++ b/electrolyzer/components/cell/pem_cell.py @@ -15,8 +15,10 @@ class PEMCellConfig(CellBaseConfig): f1: float = field(default=250, validator=validators.ge(0)) f2: float = field(default=0.996, validator=range_val(0.5, 1.0)) - i_0a: float = field(default=2.0e-7, validator=validators.optional(range_val(0.0, 1.0))) - i_0c: float = field(default=2.0e-3, validator=validators.optional(range_val(0.0, 1.0))) + # i_0a: float = field(default=2.0e-7, validator=validators.optional(range_val(0.0, 1.0))) + # i_0c: float = field(default=2.0e-3, validator=validators.optional(range_val(0.0, 1.0))) + i_0a: float = field(default=None, validator=validators.optional(range_val(0.0, 1.0))) + i_0c: float = field(default=None, validator=validators.optional(range_val(0.0, 1.0))) alpha_a: float = field(default=2.0, validator=range_val(0.25, 4.0)) alpha_c: float = field(default=0.5, validator=range_val(0.25, 4.0)) @@ -40,13 +42,6 @@ class PEMCellConfig(CellBaseConfig): ) def __attrs_post_init__(self): - provided_per_electrode_info = all( - getattr(self, param, None) is not None for param in ["i_0a", "i_0c"] - ) - provided_combined_info = all( - getattr(self, param, None) is not None for param in ["b_combined", "i0_combined"] - ) - extra_inputs_msg = ( "Extraneous inputs ({extraneous_attrs_msg}) for kinetics_method {kinetics_method}. " ) @@ -55,55 +50,53 @@ def __attrs_post_init__(self): ) if self.kinetics_method == "per_electrode": - if not provided_per_electrode_info or provided_combined_info: - # check if theres missing info or extraneous info - required_attributes = ["i_0a", "i_0c"] - extreanous_attributes = ["b_combined", "i0_combined"] - extraneous_attrs_msg = ", ".join( - f"`{k}`" for k in extreanous_attributes if getattr(self, k, None) is not None + # check if theres missing info or extraneous info + required_attributes = ["i_0a", "i_0c"] + extreanous_attributes = ["b_combined", "i0_combined"] + extraneous_attrs_msg = ", ".join( + f"`{k}`" for k in extreanous_attributes if getattr(self, k, None) is not None + ) + required_attrs_msg = ", ".join( + f"`{k}`" for k in required_attributes if getattr(self, k, None) is None + ) + txt = "" + if len(extraneous_attrs_msg) > 0: + txt += extra_inputs_msg.format( + extraneous_attrs_msg=extraneous_attrs_msg, + kinetics_method=self.kinetics_method, ) - required_attrs_msg = ", ".join( - f"`{k}`" for k in required_attributes if getattr(self, k, None) is None + + if len(required_attrs_msg) > 0: + txt += missing_inputs_msg.format( + required_attrs_msg=required_attrs_msg, kinetics_method=self.kinetics_method ) - txt = "" - if len(extraneous_attrs_msg) > 0: - txt += extra_inputs_msg.format( - extraneous_attrs_msg=extraneous_attrs_msg, - kinetics_method=self.kinetics_method, - ) - - if len(required_attrs_msg) > 0: - txt += missing_inputs_msg.format( - required_attrs_msg=required_attrs_msg, kinetics_method=self.kinetics_method - ) - if len(txt) > 0: - raise AttributeError(txt) + if len(txt) > 0: + raise AttributeError(txt) if self.kinetics_method == "combined": - if not provided_combined_info or provided_per_electrode_info: - # check if theres missing info or extraneous info - required_attributes = ["b_combined", "i0_combined"] - extreanous_attributes = ["i_0a", "i_0c"] - extraneous_attrs_msg = ", ".join( - f"`{k}`" for k in extreanous_attributes if getattr(self, k, None) is not None + # check if theres missing info or extraneous info + required_attributes = ["b_combined", "i0_combined"] + extreanous_attributes = ["i_0a", "i_0c"] + extraneous_attrs_msg = ", ".join( + f"`{k}`" for k in extreanous_attributes if getattr(self, k, None) is not None + ) + required_attrs_msg = ", ".join( + f"`{k}`" for k in required_attributes if getattr(self, k, None) is None + ) + txt = "" + if len(extraneous_attrs_msg) > 0: + txt += extra_inputs_msg.format( + extraneous_attrs_msg=extraneous_attrs_msg, + kinetics_method=self.kinetics_method, ) - required_attrs_msg = ", ".join( - f"`{k}`" for k in required_attributes if getattr(self, k, None) is None + + if len(required_attrs_msg) > 0: + txt += missing_inputs_msg.format( + required_attrs_msg=required_attrs_msg, kinetics_method=self.kinetics_method ) - txt = "" - if len(extraneous_attrs_msg) > 0: - txt += extra_inputs_msg.format( - extraneous_attrs_msg=extraneous_attrs_msg, - kinetics_method=self.kinetics_method, - ) - - if len(required_attrs_msg) > 0: - txt += missing_inputs_msg.format( - required_attrs_msg=required_attrs_msg, kinetics_method=self.kinetics_method - ) - - if len(txt) > 0: - raise AttributeError(txt) + + if len(txt) > 0: + raise AttributeError(txt) class PEMCell(CellBaseClass): From 23bf5442158757657c6557b9b0dbaa33d07a10b1 Mon Sep 17 00:00:00 2001 From: elenya-grant <116225007+elenya-grant@users.noreply.github.com> Date: Mon, 6 Jul 2026 10:31:16 -0600 Subject: [PATCH 4/7] minor updates to PEM cell --- electrolyzer/components/cell/pem_cell.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/electrolyzer/components/cell/pem_cell.py b/electrolyzer/components/cell/pem_cell.py index 71f63db..777ff09 100644 --- a/electrolyzer/components/cell/pem_cell.py +++ b/electrolyzer/components/cell/pem_cell.py @@ -101,7 +101,7 @@ def __attrs_post_init__(self): class PEMCell(CellBaseClass): def setup(self): - self.n = 2 # number of electrons transferred in reaction + # self.n = 2 # number of electrons transferred in reaction self.config = PEMCellConfig.from_dict(self.options["tech_config"]["cell_parameters"]) @@ -128,7 +128,7 @@ def calculate_current(self, J_cell, A_cell): return J_cell * A_cell def calc_Urev0(self, temperature_celsius): - E_cell0 = gibbs / (self.n * F) # 1.229 # [V] + E_cell0 = gibbs / (2 * F) # 1.229 # [V] # Reversible potential at 25degC - Nerst Equation if self.config.Urev0_calc_method == "normal": @@ -153,9 +153,7 @@ def reversible_overpotential(self, temperature_celsius, p_cathode_bar, p_anode_b # Nerst Equation Urev0 = self.calc_Urev0(temperature_celsius) - E_cell = Urev0 + ((R * temp_k) / (self.n * F)) * ( - np.log((p_H2 * np.sqrt(p_O2)) / p_H2O_sat) - ) + E_cell = Urev0 + ((R * temp_k) / (2 * F)) * (np.log((p_H2 * np.sqrt(p_O2)) / p_H2O_sat)) return E_cell @@ -290,3 +288,6 @@ def conversion_efficiency(self, inputs): h2_grams_per_sec = self.h2_production_rate(inputs) power_W_per_sec = self.power_consumption_rate(inputs) return power_W_per_sec / np.max([1e-30, h2_grams_per_sec]) + + def compute(self, inputs, outputs, discrete_inputs, discrete_outputs): + self.cell_voltage(inputs) From 92bdc43015332a3949889f0c3f81e09e5cec8af3 Mon Sep 17 00:00:00 2001 From: elenya-grant <116225007+elenya-grant@users.noreply.github.com> Date: Mon, 6 Jul 2026 11:05:34 -0600 Subject: [PATCH 5/7] added tests for pem cell --- electrolyzer/components/cell/pem_cell.py | 25 ++- .../components/cell/test/test_pem_cell.py | 173 ++++++++++++++++++ 2 files changed, 193 insertions(+), 5 deletions(-) create mode 100644 electrolyzer/components/cell/test/test_pem_cell.py diff --git a/electrolyzer/components/cell/pem_cell.py b/electrolyzer/components/cell/pem_cell.py index 777ff09..1b2e07e 100644 --- a/electrolyzer/components/cell/pem_cell.py +++ b/electrolyzer/components/cell/pem_cell.py @@ -165,8 +165,8 @@ def separate_activation_overpotential( ba = R * temp_K / (alpha_a * F) bc = R * temp_K / (alpha_c * F) - anode_ratio = np.clip(J_cell / i_0a, a_min=1e-30) - cathode_ratio = np.clip(J_cell / i_0c, a_min=1e-30) + anode_ratio = np.maximum(J_cell / i_0a, 1e-30) + cathode_ratio = np.maximum(J_cell / i_0c, 1e-30) if self.config.activation_method == "arcsinh": V_acta = ba * np.arcsinh(anode_ratio) @@ -180,7 +180,7 @@ def separate_activation_overpotential( return V_acta + V_actc def combined_kinetics_activation_overpotential(self, J_cell, b_combined, i_0combined): - ratio = np.clip(J_cell / i_0combined, a_min=1e-30) + ratio = np.maximum(J_cell / i_0combined, 1e-30) if self.config.activation_method == "ln": V_act = b_combined * np.log(ratio) if self.config.activation_method == "log10": @@ -289,5 +289,20 @@ def conversion_efficiency(self, inputs): power_W_per_sec = self.power_consumption_rate(inputs) return power_W_per_sec / np.max([1e-30, h2_grams_per_sec]) - def compute(self, inputs, outputs, discrete_inputs, discrete_outputs): - self.cell_voltage(inputs) + def compute(self, inputs, outputs): + V_cell = self.cell_voltage(inputs) + J_cell = self.get_current_density(inputs) + outputs["cell_voltage"] = V_cell + outputs["current_density_out"] = J_cell + + outputs["hydrogen_produced"] = self.h2_production(inputs) + outputs["hydrogen_production_rate"] = self.h2_production_rate(inputs) + outputs["oxygen_produced"] = self.o2_production(inputs) + outputs["oxygen_production_rate"] = self.o2_production_rate(inputs) + # outputs["water_consumed"] + + # outputs["rated_cell_voltage"] + # outputs["rated_conversion_efficiency"] + # outputs["rated_hydrogen_production"] + # outputs["rated_oxygen_production"] + # outputs["rated_cell_power"] diff --git a/electrolyzer/components/cell/test/test_pem_cell.py b/electrolyzer/components/cell/test/test_pem_cell.py new file mode 100644 index 0000000..fd14415 --- /dev/null +++ b/electrolyzer/components/cell/test/test_pem_cell.py @@ -0,0 +1,173 @@ +from types import SimpleNamespace + +import numpy as np +import pytest +import openmdao.api as om +from pytest import fixture + +from electrolyzer.components.cell.pem_cell import PEMCell + + +@fixture +def plant_config(): + return { + "simulation": { + "n_timesteps": 20, # unused by cell at the moment + "dt": 3600, + } + } + + +@fixture +def pem_cell_config(kinetics_method, water_vapor_method, act_method, Urev0_method): + config = { + "A_cell": 1000, + "membrane_thickness": 0.02, + "temperature": 80.0, + "P_anode": 1.0, + "P_cathode": 1.0, + "R_ohmic_elec": 50.0e-3, + "f1": 250, + "f2": 0.996, + "i_0a": 4.0e-7, + "i_0c": 4.0e-3, + "alpha_a": 2.0, + "alpha_c": 0.5, + "water_vaporization_method": water_vapor_method, # "antoine", # "buck", "antoine" + "activation_method": act_method, # "ln", # "ln", "log10", or "arcsinh" + "kinetics_method": kinetics_method, # "per_electrode", # "per_electrode" or "combined" + "Urev0_calc_method": Urev0_method, # "normal", # "normal" or "temp_adjusted" + } + return config + + +@fixture +def cell_namespace_config(kinetics_method, water_vapor_method, act_method, Urev0_method): + return SimpleNamespace( + kinetics_method=kinetics_method, + water_vaporization_method=water_vapor_method, + activation_method=act_method, + Urev0_calc_method=Urev0_method, + ) + + +@pytest.mark.parametrize( + "kinetics_method,water_vapor_method,act_method,Urev0_method", + [("per_electrode", "buck", "arcsinh", "normal")], +) +def test_Urev_buck_normal(cell_namespace_config, subtests): + cell = object.__new__(PEMCell) + cell.config = cell_namespace_config + Urev0 = cell.calc_Urev0(None) + with subtests.test("Urev0 - non-temperature dependent"): + assert pytest.approx(1.229, rel=1e-6, abs=1e-3) == Urev0 + + with subtests.test("Water vaporization pressure (25 C, buck)"): + p_H2O_sat_25C = cell.water_vaporization_pressure_arden_buck(25.0) + assert pytest.approx(0.031685314, rel=1e-6, abs=1e-5) == p_H2O_sat_25C + + with subtests.test("Water vaporization pressure (60 C, buck)"): + p_H2O_sat_60C = cell.water_vaporization_pressure_arden_buck(60.0) + assert pytest.approx(0.19945145, rel=1e-6, abs=1e-5) == p_H2O_sat_60C + + with subtests.test("Water vaporization pressure (80 C, buck)"): + p_H2O_sat_80C = cell.water_vaporization_pressure_arden_buck(80.0) + assert pytest.approx(0.47410267, rel=1e-6, abs=1e-5) == p_H2O_sat_80C + + with subtests.test("Urev - 25 C, 1 bar (buck)"): + Urev_25C = cell.reversible_overpotential(25.0, 1.0, 1.0) + assert pytest.approx(1.273133, rel=1e-6, abs=1e-5) == Urev_25C + + with subtests.test("Urev - 60 C, 1 bar (buck)"): + Urev_60C = cell.reversible_overpotential(60.0, 1.0, 1.0) + assert pytest.approx(1.24776, rel=1e-6, abs=1e-5) == Urev_60C + + with subtests.test("Urev - 80 C, 1 bar (buck)"): + Urev_80C = cell.reversible_overpotential(80.0, 1.0, 1.0) + assert pytest.approx(1.226098, rel=1e-6, abs=1e-5) == Urev_80C + + with subtests.test("Urev - 80 C, 30 bar cathode, 1 bar anode (buck)"): + Urev_80C_30_1 = cell.reversible_overpotential(80.0, 30.0, 1.0) + assert pytest.approx(1.287387, rel=1e-6, abs=1e-5) == Urev_80C_30_1 + + with subtests.test("Urev - 80 C, 1 bar cathode, 30 bar anode (buck)"): + Urev_80C_1_30 = cell.reversible_overpotential(80.0, 1.0, 30.0) + assert pytest.approx(1.2567425, rel=1e-6, abs=1e-5) == Urev_80C_1_30 + + with subtests.test("Urev - 80 C, 10 bar cathode, 30 bar anode (buck)"): + Urev_80C_10_30 = cell.reversible_overpotential(80.0, 10.0, 30.0) + assert pytest.approx(1.300818, rel=1e-6, abs=1e-5) == Urev_80C_10_30 + + +@pytest.mark.parametrize( + "kinetics_method,water_vapor_method,act_method,Urev0_method", + [("per_electrode", "buck", "arcsinh", "normal")], +) +def test_pem_cell_default(plant_config, pem_cell_config, subtests): + prob = om.Problem() + comp = PEMCell(plant_config=plant_config, tech_config={"cell_parameters": pem_cell_config}) + prob.model.add_subsystem("cell", comp, promotes=["*"]) + prob.setup() + prob.set_val("cell.current_in", np.full(20, 1.0 * 1000), units="A") + prob.run_model() + + with subtests.test("Current density is 1 A/cm2"): + J_cell = prob.get_val("cell.current_density_out", units="A/(cm**2)") + assert np.all(J_cell == 1.0) + + with subtests.test("Cell voltage at 1 A/cm2"): + assert ( + pytest.approx(prob.get_val("cell.cell_voltage", units="V")[0], rel=1e-6, abs=1e-5) + == 1.983580501 + ) + + with subtests.test("H2 Production at 1 A/cm2"): + assert ( + pytest.approx( + prob.get_val("cell.hydrogen_produced", units="g/h")[0], rel=1e-6, abs=1e-5 + ) + == 37.45005976 + ) + + with subtests.test("H2 Production Rate at 1 A/cm2"): + assert ( + pytest.approx( + prob.get_val("hydrogen_production_rate", units="g/h")[0], rel=1e-6, abs=1e-5 + ) + == 37.45005976 + ) + + with subtests.test("O2 Production at 1 A/cm2"): + assert ( + pytest.approx(prob.get_val("cell.oxygen_produced", units="g/h")[0], rel=1e-6, abs=1e-5) + == 297.20412014327 + ) + + with subtests.test("O2 Production Rate at 1 A/cm2"): + assert ( + pytest.approx( + prob.get_val("oxygen_production_rate", units="g/h")[0], rel=1e-6, abs=1e-5 + ) + == 297.20412014327 + ) + + power = prob.get_val("cell.cell_voltage", units="V") * prob.get_val( + "cell.current_in", units="A" + ) + eff = (power / 1e3) / (prob.get_val("hydrogen_production_rate", units="kg/h")) + with subtests.test("H2 Conversion efficiency at 1 A/cm2"): + assert pytest.approx(eff[0], rel=1e-6) == 52.96601698232939 + + +# def test_pem_cell_om(pem_cell_config, subtests): + +# prob = om.Problem() + +# comp = PEMCell( +# plant_config=plant_config, +# tech_config=pem_cell_config, +# ) +# prob.model.add_subsystem("cell", comp, promotes=["*"]) +# prob.setup() +# prob.set_val("cell.current_in", current_vals, units="A") +# prob.run_model() From b2bc650f38d0e283fac09a4d0574f782e45803b0 Mon Sep 17 00:00:00 2001 From: elenya-grant <116225007+elenya-grant@users.noreply.github.com> Date: Thu, 9 Jul 2026 16:38:00 -0600 Subject: [PATCH 6/7] renamed inputs/outputs of cell --- electrolyzer/components/cell/cell.py | 43 ++++++++-------- electrolyzer/components/cell/pem_cell.py | 51 +++++++++---------- .../components/cell/test/test_pem_cell.py | 29 ++++------- 3 files changed, 53 insertions(+), 70 deletions(-) diff --git a/electrolyzer/components/cell/cell.py b/electrolyzer/components/cell/cell.py index f48dbcc..643b177 100644 --- a/electrolyzer/components/cell/cell.py +++ b/electrolyzer/components/cell/cell.py @@ -27,11 +27,11 @@ def initialize(self): self.options.declare("tech_config", types=dict) def setup(self): - self.n_timesteps = self.options["plant_config"]["simulation"]["n_timesteps"] + # self.n_timesteps = self.options["plant_config"]["simulation"]["n_timesteps"] self.dt = self.options["plant_config"]["simulation"]["dt"] # design variables - self.add_input("cell_active_area", val=self.config.A_cell, units="A/(cm**2)") + self.add_input("A_cell", val=self.config.A_cell, units="A/(cm**2)") self.add_input("membrane_thickness", val=self.config.membrane_thickness, units="cm") self.add_input("operating_temperature", val=self.config.temperature, units="degC") self.add_input("anode_pressure", self.config.P_anode, units="bar") @@ -40,31 +40,28 @@ def setup(self): # input profiles # TODO: make either current density or current - self.add_input("current_in", val=0.0, shape_by_conn=True, units="A") # OR current density? - # self.add_input("current_density_in", val=0.0, shape_by_conn=True, units="A/(cm**2)") + self.add_input("I_in", val=0.0, shape_by_conn=True, units="A") # OR current density? + # self.add_input("J_in", val=0.0, shape_by_conn=True, units="A/(cm**2)") # output profiles - self.add_output("cell_voltage", val=0.0, copy_shape="current_in", units="V") - self.add_output( - "hydrogen_produced", val=0.0, copy_shape="current_in", units=f"g/({self.dt}*s)" - ) - self.add_output( - "oxygen_produced", val=0.0, copy_shape="current_in", units=f"g/({self.dt}*s)" - ) - self.add_output( - "water_consumed", val=0.0, copy_shape="current_in", units=f"g/({self.dt}*s)" - ) - self.add_output("hydrogen_production_rate", val=0.0, copy_shape="current_in", units="g/s") - self.add_output("oxygen_production_rate", val=0.0, copy_shape="current_in", units="g/s") - - self.add_output("current_density_out", val=0.0, copy_shape="current_in", units="A/(cm**2)") + self.add_output("V_cell", val=0.0, copy_shape="I_in", units="V") + self.add_output("H2_produced", val=0.0, copy_shape="I_in", units=f"g/({self.dt}*s)") + self.add_output("O2_produced", val=0.0, copy_shape="I_in", units=f"g/({self.dt}*s)") + # self.add_output( + # "H2O_consumed", val=0.0, copy_shape="I_in", units=f"g/({self.dt}*s)" + # ) + self.add_output("H2_out", val=0.0, copy_shape="I_in", units="g/s") + self.add_output("O2_out", val=0.0, copy_shape="I_in", units="g/s") + + self.add_output("J_out", val=0.0, copy_shape="I_in", units="A/(cm**2)") + self.add_output("P_cell", val=0.0, copy_shape="I_in", units="W") # output design variables - self.add_output("rated_cell_voltage", val=0.0, units="V") - self.add_output("rated_conversion_efficiency", val=0.0, units="W*h/kg") - self.add_output("rated_hydrogen_production", val=0.0, units="g/s") - self.add_output("rated_oxygen_production", val=0.0, units="g/s") - self.add_output("rated_cell_power", val=0.0, units="W") + # self.add_output("rated_V_cell", val=0.0, units="V") + # self.add_output("rated_conversion_efficiency", val=0.0, units="W*h/kg") + # self.add_output("rated_H2_production", val=0.0, units="g/s") + # self.add_output("rated_O2_production", val=0.0, units="g/s") + # self.add_output("rated_cell_power", val=0.0, units="W") def compute(self, inputs, outputs, discrete_inputs, discrete_outputs): """ diff --git a/electrolyzer/components/cell/pem_cell.py b/electrolyzer/components/cell/pem_cell.py index 1b2e07e..347e4fe 100644 --- a/electrolyzer/components/cell/pem_cell.py +++ b/electrolyzer/components/cell/pem_cell.py @@ -206,12 +206,10 @@ def ohmic_overpotential(self, J_cell, temperature_celsius, delta_membrane, R_ele def cell_voltage(self, inputs): temp_C = inputs["operating_temperature"] - if "current_density_in" in inputs: - J_cell = inputs["current_density_in"] - elif "current_in" in inputs: - J_cell = self.calculate_current_density( - inputs["current_in"], inputs["cell_active_area"] - ) + if "J_in" in inputs: + J_cell = inputs["J_in"] + elif "I_in" in inputs: + J_cell = self.calculate_current_density(inputs["I_in"], inputs["A_cell"]) Urev = self.reversible_overpotential( temp_C, inputs["cathode_pressure"], inputs["anode_pressure"] @@ -233,21 +231,17 @@ def cell_voltage(self, inputs): return V_cell def get_current_density(self, inputs): - if "current_density_in" in inputs: - return inputs["current_density_in"] - if "current_in" in inputs: - J_cell = self.calculate_current_density( - inputs["current_in"], inputs["cell_active_area"] - ) + if "J_in" in inputs: + return inputs["J_in"] + if "I_in" in inputs: + J_cell = self.calculate_current_density(inputs["I_in"], inputs["A_cell"]) return J_cell def get_current(self, inputs): - if "current_in" in inputs: - return inputs["current_in"] - if "current_density_in" in inputs: - I_cell = self.calculate_current( - inputs["current_density_in"], inputs["cell_active_area"] - ) + if "I_in" in inputs: + return inputs["I_in"] + if "J_in" in inputs: + I_cell = self.calculate_current(inputs["J_in"], inputs["A_cell"]) return I_cell def h2_production_rate(self, inputs): @@ -292,17 +286,18 @@ def conversion_efficiency(self, inputs): def compute(self, inputs, outputs): V_cell = self.cell_voltage(inputs) J_cell = self.get_current_density(inputs) - outputs["cell_voltage"] = V_cell - outputs["current_density_out"] = J_cell + outputs["V_cell"] = V_cell + outputs["J_out"] = J_cell - outputs["hydrogen_produced"] = self.h2_production(inputs) - outputs["hydrogen_production_rate"] = self.h2_production_rate(inputs) - outputs["oxygen_produced"] = self.o2_production(inputs) - outputs["oxygen_production_rate"] = self.o2_production_rate(inputs) - # outputs["water_consumed"] + outputs["H2_produced"] = self.h2_production(inputs) + outputs["H2_out"] = self.h2_production_rate(inputs) + outputs["O2_produced"] = self.o2_production(inputs) + outputs["O2_out"] = self.o2_production_rate(inputs) + outputs["P_cell"] = V_cell * inputs["I_in"] + # outputs["H2O_consumed"] - # outputs["rated_cell_voltage"] + # outputs["rated_V_cell"] # outputs["rated_conversion_efficiency"] - # outputs["rated_hydrogen_production"] - # outputs["rated_oxygen_production"] + # outputs["rated_H2_production"] + # outputs["rated_O2_production"] # outputs["rated_cell_power"] diff --git a/electrolyzer/components/cell/test/test_pem_cell.py b/electrolyzer/components/cell/test/test_pem_cell.py index fd14415..8b9f43f 100644 --- a/electrolyzer/components/cell/test/test_pem_cell.py +++ b/electrolyzer/components/cell/test/test_pem_cell.py @@ -108,53 +108,44 @@ def test_pem_cell_default(plant_config, pem_cell_config, subtests): comp = PEMCell(plant_config=plant_config, tech_config={"cell_parameters": pem_cell_config}) prob.model.add_subsystem("cell", comp, promotes=["*"]) prob.setup() - prob.set_val("cell.current_in", np.full(20, 1.0 * 1000), units="A") + prob.set_val("cell.I_in", np.full(20, 1.0 * 1000), units="A") prob.run_model() with subtests.test("Current density is 1 A/cm2"): - J_cell = prob.get_val("cell.current_density_out", units="A/(cm**2)") + J_cell = prob.get_val("cell.J_out", units="A/(cm**2)") assert np.all(J_cell == 1.0) with subtests.test("Cell voltage at 1 A/cm2"): assert ( - pytest.approx(prob.get_val("cell.cell_voltage", units="V")[0], rel=1e-6, abs=1e-5) + pytest.approx(prob.get_val("cell.V_cell", units="V")[0], rel=1e-6, abs=1e-5) == 1.983580501 ) with subtests.test("H2 Production at 1 A/cm2"): assert ( - pytest.approx( - prob.get_val("cell.hydrogen_produced", units="g/h")[0], rel=1e-6, abs=1e-5 - ) + pytest.approx(prob.get_val("cell.H2_produced", units="g/h")[0], rel=1e-6, abs=1e-5) == 37.45005976 ) with subtests.test("H2 Production Rate at 1 A/cm2"): assert ( - pytest.approx( - prob.get_val("hydrogen_production_rate", units="g/h")[0], rel=1e-6, abs=1e-5 - ) - == 37.45005976 + pytest.approx(prob.get_val("H2_out", units="g/h")[0], rel=1e-6, abs=1e-5) == 37.45005976 ) with subtests.test("O2 Production at 1 A/cm2"): assert ( - pytest.approx(prob.get_val("cell.oxygen_produced", units="g/h")[0], rel=1e-6, abs=1e-5) + pytest.approx(prob.get_val("cell.O2_produced", units="g/h")[0], rel=1e-6, abs=1e-5) == 297.20412014327 ) with subtests.test("O2 Production Rate at 1 A/cm2"): assert ( - pytest.approx( - prob.get_val("oxygen_production_rate", units="g/h")[0], rel=1e-6, abs=1e-5 - ) + pytest.approx(prob.get_val("O2_out", units="g/h")[0], rel=1e-6, abs=1e-5) == 297.20412014327 ) - power = prob.get_val("cell.cell_voltage", units="V") * prob.get_val( - "cell.current_in", units="A" - ) - eff = (power / 1e3) / (prob.get_val("hydrogen_production_rate", units="kg/h")) + power = prob.get_val("cell.V_cell", units="V") * prob.get_val("cell.I_in", units="A") + eff = (power / 1e3) / (prob.get_val("H2_out", units="kg/h")) with subtests.test("H2 Conversion efficiency at 1 A/cm2"): assert pytest.approx(eff[0], rel=1e-6) == 52.96601698232939 @@ -169,5 +160,5 @@ def test_pem_cell_default(plant_config, pem_cell_config, subtests): # ) # prob.model.add_subsystem("cell", comp, promotes=["*"]) # prob.setup() -# prob.set_val("cell.current_in", current_vals, units="A") +# prob.set_val("cell.I_in", current_vals, units="A") # prob.run_model() From 611b50faf715930decf1a10ce9596b16eaebad24 Mon Sep 17 00:00:00 2001 From: elenya-grant <116225007+elenya-grant@users.noreply.github.com> Date: Thu, 9 Jul 2026 16:48:35 -0600 Subject: [PATCH 7/7] added constants file --- electrolyzer/components/constants.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) create mode 100644 electrolyzer/components/constants.py diff --git a/electrolyzer/components/constants.py b/electrolyzer/components/constants.py new file mode 100644 index 0000000..101646f --- /dev/null +++ b/electrolyzer/components/constants.py @@ -0,0 +1,12 @@ +from scipy.constants import physical_constants + + +# R is Ideal Gas Constant (J/mol/K) +F, _, _ = physical_constants["Faraday constant"] # Faraday's constant [C/mol] + +H2_MW: float = 2.016 # molecular weight [g/mol] +O2_MW = 31.998 # Molecular weight of Oxygen in g/mol + +gibbs = 237.24e3 # Gibbs Energy of global reaction (J/mol) +H2_LHV_kWh_per_kg: float = 33.33 # lower heating value of H2 [kWh/kg] +H2_HHV_kWh_per_kg: float = 39.41 # higher heating value of H2 [kWh/kg]