diff --git a/electrolyzer/components/cell/cell.py b/electrolyzer/components/cell/cell.py index e568084..643b177 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) @@ -20,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") @@ -32,25 +39,29 @@ def setup(self): self.add_input("R_elec", self.config.R_ohmic_elec, units="ohm*(cm**2)") # input profiles - 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)") + # TODO: make either current density or current + 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"kg/({self.dt}*s)" - ) - self.add_output( - "oxygen_produced", val=0.0, copy_shape="current_in", units=f"kg/({self.dt}*s)" - ) - self.add_output( - "water_consumed", val=0.0, copy_shape="current_in", units=f"kg/({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("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_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): """ @@ -60,3 +71,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..347e4fe 100644 --- a/electrolyzer/components/cell/pem_cell.py +++ b/electrolyzer/components/cell/pem_cell.py @@ -1,17 +1,29 @@ +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))) + 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)) + + 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 +31,273 @@ 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): + 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": + # 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": + # 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 / (2 * 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) / (2 * 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.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) + 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.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": + 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 "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"] + ) + 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 "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 "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): + 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]) + + def compute(self, inputs, outputs): + V_cell = self.cell_voltage(inputs) + J_cell = self.get_current_density(inputs) + outputs["V_cell"] = V_cell + outputs["J_out"] = J_cell + + 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_V_cell"] + # outputs["rated_conversion_efficiency"] + # 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 new file mode 100644 index 0000000..8b9f43f --- /dev/null +++ b/electrolyzer/components/cell/test/test_pem_cell.py @@ -0,0 +1,164 @@ +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.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.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.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.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("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.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("O2_out", units="g/h")[0], rel=1e-6, abs=1e-5) + == 297.20412014327 + ) + + 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 + + +# 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.I_in", current_vals, units="A") +# prob.run_model() 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 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]