From 4551f87b7c0578f65915f573875e5027eb3db99e Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Wed, 20 May 2026 14:10:05 -0400 Subject: [PATCH 01/21] Initial fuel cell model with psuedo code --- .../converters/hydrogen/PEM_h2_fuel_cell.py | 296 ++++++++++++++++++ 1 file changed, 296 insertions(+) create mode 100644 h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py diff --git a/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py b/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py new file mode 100644 index 000000000..9a5a204df --- /dev/null +++ b/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py @@ -0,0 +1,296 @@ +from attrs import field, define + +from h2integrate.core.utilities import BaseConfig, merge_shared_inputs +from h2integrate.core.validators import gte_zero, range_val +from h2integrate.core.model_baseclasses import ( + CostModelBaseClass, + CostModelBaseConfig, + PerformanceModelBaseClass, +) + + +@define(kw_only=True) +class PEMH2FuelCellPerformanceConfig(BaseConfig): + """Configuration class for the hydrogen fuel cell performance model. + + Attributes: + system_capacity_kw (float): The capacity of the fuel cell system in kilowatts (kW). + fuel_cell_efficiency_hhv (float): The higher heating value efficiency of the + fuel cell (0 <= efficiency <= 1). + """ + + # TODO: how to size the fuel cell? N_cells + N_stacks? + # How does N_cells translate to electricity rating? + + system_capacity_kw: float = field(validator=gte_zero) + stack_temperature_K: float + fuel_cell_efficiency_hhv: float = field(validator=range_val(0, 1)) + + +class PEMH2FuelCellPerformanceModel(PerformanceModelBaseClass): + """ + Performance model for a hydrogen fuel cell. + + The model implements the relationship: + electricity_out = hydrogen_in * fuel_cell_efficiency_hhv * HHV_hydrogen + + where: + - hydrogen_in is the mass flow rate of hydrogen in kg/hr + - fuel_cell_efficiency is the efficiency of the fuel cell (0 <= efficiency <= 1) + - HHV_hydrogen is the higher heating value of hydrogen (approximately 142 MJ/kg) + """ + + _time_step_bounds = ( + 3600, + 3600, + ) # (min, max) time step lengths (in seconds) compatible with this model + + def initialize(self): + super().initialize() + self.commodity = "electricity" + self.commodity_rate_units = "kW" + self.commodity_amount_units = "kW*h" + + def setup(self): + super().setup() + + self.config = PEMH2FuelCellPerformanceConfig.from_dict( + merge_shared_inputs(self.options["tech_config"]["model_inputs"], "performance"), + additional_cls_name=self.__class__.__name__, + ) + + # Add natural gas input, default to 0 --> set using feedstock component + # or upstream hydrogen converter component + self.add_input( + "hydrogen_in", + val=0.0, + shape=self.n_timesteps, + units="kg/h", + ) + + self.add_input( + "oxygen_in", + val=0.0, + shape=self.n_timesteps, + units="kg/h", + ) + + self.add_input( + "stack_temperature", + val=self.config.stack_temperature_K, + units="K", + desc="Operating temperature of the stack", + ) + + # self.add_input( + # "fuel_cell_efficiency", + # val=self.config.fuel_cell_efficiency_hhv, + # units=None, + # desc="HHV efficiency of the fuel cell (0 <= efficiency <= 1)", + # ) + + # Add rated capacity as an input with config value as default + self.add_input( + "system_capacity", + val=self.config.system_capacity_kw, + units="kW", + desc="Capacity of the h2 fuel cell system", + ) + + self.add_output( + "hydrogen_consumed", + val=0.0, + shape=self.n_timesteps, + units="kg/h", + desc="Mass flow rate of hydrogen consumed by the fuel cell", + ) + + self.add_output( + "oxygen_consumed", + val=0.0, + shape=self.n_timesteps, + units="kg/h", + desc="Mass flow rate of oxygen consumed by the fuel cell", + ) + + self.add_output( + "water_out", + val=0.0, + shape=self.n_timesteps, + units="kg/h", + desc="Mass flow rate of water consumed by the fuel cell", + ) + + # Default the electricity set point input as the rated capacity + self.add_input( + f"{self.commodity}_set_point", + val=self.config.system_capacity_kw, + shape=self.n_timesteps, + units=self.commodity_rate_units, + desc="Electricity set point for natural gas plant", + ) + + def compute(self, inputs, outputs): + """ + Compute electricity output from the fuel cell based on hydrogen input + and fuel cell HHV efficiency. + + Args: + inputs: OpenMDAO inputs object containing hydrogen_in, fuel cell + HHV efficiency, electricity_set_point, and system_capacity. + outputs: OpenMDAO outputs object for electricity_out, + hydrogen_consumed. + """ + + # calculate max input and output + inputs["system_capacity"] # plant capacity in kW + inputs["hydrogen_in"] # kg/h + inputs["oxygen_in"] # kg/h + inputs["stack_temperature"] + # fuel_cell_efficiency = inputs["fuel_cell_efficiency"] + + # Set calculation constants: + self.f_c = 96485.33 # Faraday's constant in Coulombs/mol + self.M_H2 = 0.002016 # Molar mass of H2 in kg/mol + self.M_O2 = 0.32 # Molar mass of O2 in kg/mol + self.Tref = 298.15 # Standard room temperature in K [25 deg Celsius] + self.cp_H2 = 14300 # Specific heat of H2 in J/(kg*K) + self.cp_air = 1005 # Specific heat of air in J/(kg*K) + self.cp_water = 4184 # Specific heat of water in J/(kg*K) + self.cp_N2 = 1040 # Specific heat of nitrogen in J/(kg*K) + self.hhv_h2 = 141.8 * 1e6 # Higher heating value of hydrogen in J/kg + self.hhv_air = 0 # No higher heating value of air + self.hhv_water = 0 # ???? Need to check this + + # PSUEDO CODE: + """ + 1. Receive power setpoint into fuel cell + 2. Find current with I-V curve - NEED THIS + 3. Calculate power out with current - NEED TO FIND CELL ELECTRICITY CALC + 4. Calculate H2 consumed, O2 consumed, water out + 5. See if H2 in and O2 in can provide this + 6. Repeat step 4 if H2 or O2 limit power out + + """ + + # # conversion factor: kW electricity to kg/h hydrogen, units: (kg/h)/kW + # kw_to_kgh_h2 = (3600.0 * 0.001) / (fuel_cell_efficiency * HHV_H2_MJ_PER_KG) + + # # TODO: Calculate max H2 and O2 consumption of the stacks + # max_h2_consumption = system_capacity * kw_to_kgh_h2 + + # # electrical set point, saturated at maximum rated system capacity + # electricity_set_point = np.where( + # inputs["electricity_set_point"] > system_capacity, + # system_capacity, + # inputs["electricity_set_point"], + # ) + + # h2_demand = electricity_set_point * kw_to_kgh_h2 + + # # available feedstock, saturated at maximum system feedstock consumption + # h2_available = np.where( + # inputs["hydrogen_in"] > max_h2_consumption, + # max_h2_consumption, + # inputs["hydrogen_in"], + # ) + + # # h2 consumed is minimum between available feedstock and output demand + # hydrogen_in = np.minimum(h2_available, h2_demand) + + # # make any negative hydrogen input zero + # hydrogen_in = np.maximum(hydrogen_in, 0.0) + + # # calculate electricity output in kW + # electricity_out_kw = hydrogen_in / kw_to_kgh_h2 + + # # clip the electricity output to the system capacity + # outputs["electricity_out"] = np.minimum(electricity_out_kw, system_capacity) + # outputs["total_electricity_produced"] = np.sum(outputs["electricity_out"]) * ( + # self.dt / 3600 + # ) + # outputs["rated_electricity_production"] = system_capacity + # outputs["annual_electricity_produced"] = outputs["total_electricity_produced"] * ( + # 1 / self.fraction_of_year_simulated + # ) + # outputs["capacity_factor"] = outputs["total_electricity_produced"] / ( + # system_capacity * self.n_timesteps * (self.dt / 3600) + # ) + # outputs["hydrogen_consumed"] = outputs["electricity_out"] * kw_to_kgh_h2 + + # def enthalpy_flow(self, m, cp, T, Tref, h0): + # """Mass-specific enthalpy flow: Hdot = m * (cp*(T - Tref) + h0)""" + # return m * (cp * (T - Tref) + h0) + + +@define(kw_only=True) +class H2FuelCellCostConfig(CostModelBaseConfig): + """Configuration class for the hydrogen fuel cell cost model. + + Fields include `system_capacity_kw`, `capex_per_kw`, and `fixed_opex_per_kw_per_year`. + The `cost_year` field is inherited from `CostModelBaseConfig`. + """ + + system_capacity_kw: float = field(validator=gte_zero) + capex_per_kw: float = field(validator=gte_zero) + fixed_opex_per_kw_per_year: float = field(validator=gte_zero) + + +class H2FuelCellCostModel(CostModelBaseClass): + """ + Cost model for a hydrogen fuel cell system. + + The model calculates capital and fixed operating costs based on system capacity and + specified cost parameters. + """ + + _time_step_bounds = ( + 3600, + 3600, + ) # (min, max) time step lengths (in seconds) compatible with this model + + def setup(self): + self.config = H2FuelCellCostConfig.from_dict( + merge_shared_inputs(self.options["tech_config"]["model_inputs"], "cost"), + additional_cls_name=self.__class__.__name__, + ) + + super().setup() + + self.add_input( + "system_capacity", + val=self.config.system_capacity_kw, + units="kW", + desc="Capacity of the h2 fuel cell system", + ) + + self.add_input( + "unit_capex", + val=self.config.capex_per_kw, + units="USD/kW", + desc="Capital cost per unit capacity", + ) + + self.add_input( + "fixed_opex_per_year", + val=self.config.fixed_opex_per_kw_per_year, + units="(USD/kW)/year", + desc="Fixed operating expenses per unit capacity per year", + ) + + def compute(self, inputs, outputs, discrete_inputs, discrete_outputs): + """ + Compute capital and fixed operating costs for the fuel cell system. + + Args: + inputs: OpenMDAO inputs object containing system_capacity. + outputs: OpenMDAO outputs object for capital_cost and fixed_operating_cost_per_year. + """ + + system_capacity_kw = inputs["system_capacity"] + + # Calculate capital cost + outputs["CapEx"] = system_capacity_kw * inputs["unit_capex"] + + # Calculate fixed operating cost per year + outputs["OpEx"] = system_capacity_kw * inputs["fixed_opex_per_year"] From 9115bbec8bc1ad5c70658927938f1bf4b2eeb088 Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Wed, 27 May 2026 16:14:59 -0400 Subject: [PATCH 02/21] Initial fuel cell model --- .../converters/hydrogen/PEM_h2_fuel_cell.py | 126 +++++++++++++++++- 1 file changed, 125 insertions(+), 1 deletion(-) diff --git a/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py b/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py index 9a5a204df..07938bc47 100644 --- a/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py +++ b/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py @@ -1,4 +1,6 @@ +import numpy as np from attrs import field, define +from scipy.interpolate import make_interp_spline from h2integrate.core.utilities import BaseConfig, merge_shared_inputs from h2integrate.core.validators import gte_zero, range_val @@ -23,10 +25,78 @@ class PEMH2FuelCellPerformanceConfig(BaseConfig): # How does N_cells translate to electricity rating? system_capacity_kw: float = field(validator=gte_zero) + n_stacks = int stack_temperature_K: float + min_system_power_fraction_kw: float fuel_cell_efficiency_hhv: float = field(validator=range_val(0, 1)) +def calc_current(power_ref, cell_area, stack_number): + # Calculates the current and voltage from IV curve based on power reference + current_curve = [ + 0.0356, + 0.05413333, + 0.0796, + 0.11366667, + 0.244, + 0.454, + 0.70366667, + 0.96933333, + 1.24, + 1.52666667, + 1.80333333, + 2.07, + 2.32, + 2.54333333, + 2.73666667, + 2.9, + ] # in A + voltage_curve = [ + 0.987, + 0.936, + 0.884, + 0.838, + 0.786, + 0.736, + 0.686, + 0.636, + 0.586, + 0.53566667, + 0.486, + 0.436, + 0.386, + 0.33533333, + 0.286, + 0.236, + ] # in V + power_curve = [ + 35.16666667, + 50.53333333, + 70.33333333, + 95.46666667, + 191.66666667, + 334.33333333, + 482.66666667, + 616.66666667, + 729.0, + 817.0, + 875.33333333, + 902.0, + 895.0, + 854.33333333, + 782.66666667, + 684.33333333, + ] / 1e6 # change from mW to kW + + power_I_curve = make_interp_spline(power_curve, current_curve, k=3) + V_I_curve = make_interp_spline(current_curve, voltage_curve, k=5) + power_density = power_ref / cell_area / stack_number + + I_cell = power_I_curve(power_density) + V_cell = V_I_curve(I_cell) + return I_cell, V_cell + + class PEMH2FuelCellPerformanceModel(PerformanceModelBaseClass): """ Performance model for a hydrogen fuel cell. @@ -127,7 +197,7 @@ def setup(self): val=self.config.system_capacity_kw, shape=self.n_timesteps, units=self.commodity_rate_units, - desc="Electricity set point for natural gas plant", + desc="Electricity set point for PEM fuel cell", ) def compute(self, inputs, outputs): @@ -162,6 +232,11 @@ def compute(self, inputs, outputs): self.hhv_air = 0 # No higher heating value of air self.hhv_water = 0 # ???? Need to check this + self.n_cells = 100 # number of cells per stack - how to size this?? + # is n_cells = N_series? + self.stack_size = self.config.system_capacity_kw / self.config.n_stacks + self.cell_active_area = 1250 # [cm^2] + # PSUEDO CODE: """ 1. Receive power setpoint into fuel cell @@ -173,6 +248,55 @@ def compute(self, inputs, outputs): """ + h2_consumed = np.zeros(self.n_timesteps) + o2_consumed = np.zeros(self.n_timesteps) + commodity_out = np.zeros(self.n_timesteps) + + for i in range(self.n_timesteps): + power_reference = inputs[f"{self.commodity}_set_point"][i] + H2in = inputs["hydrogen_in"][i] + O2in = inputs["oxygen_in"][i] + + # Find current and voltage from IV curve with power setpoint + I_cell, V_cell = calc_current(power_reference) + + # Calculate hydrogen and oxygen consumed + H2_consumed_rate = ( + I_cell * self.n_cells / (2.0 * self.f_c) * self.M_H2 * self.dt + ) # kg/time step + O2_consumed_rate = ( + I_cell * self.n_cells / (4.0 * self.f_c) * self.M_O2 * self.dt + ) # kg/time step + + if H2_consumed_rate > H2in or O2_consumed_rate > O2in: + print("Not enough H2 or O2 for this power point") + # implement an adjustment based on H2 & O2 available + + # Compute electricity from the system + electricity_produced = V_cell * I_cell * self.n_cells * self.config.n_stacks + + # Compute H2O out (is this needed?) + + h2_consumed[i] = H2_consumed_rate + o2_consumed[i] = O2_consumed_rate + commodity_out[i] = electricity_produced + + # Set Outputs + # clip the electricity output to the system capacity + outputs["electricity_out"] = np.minimum(commodity_out, self.config.system_capacity_kw) + outputs["total_electricity_produced"] = np.sum(outputs["electricity_out"]) * ( + self.dt / 3600 + ) + outputs["rated_electricity_production"] = self.config.system_capacity_kw + outputs["annual_electricity_produced"] = outputs["total_electricity_produced"] * ( + 1 / self.fraction_of_year_simulated + ) + outputs["capacity_factor"] = outputs["total_electricity_produced"] / ( + self.config.system_capacity_kw * self.n_timesteps * (self.dt / 3600) + ) + outputs["hydrogen_consumed"] = h2_consumed + outputs["oxygen_consumed"] = o2_consumed + # # conversion factor: kW electricity to kg/h hydrogen, units: (kg/h)/kW # kw_to_kgh_h2 = (3600.0 * 0.001) / (fuel_cell_efficiency * HHV_H2_MJ_PER_KG) From 9529cf6a2ef86e9592b75aad6fe83f67f744c32b Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Thu, 4 Jun 2026 09:37:45 -0400 Subject: [PATCH 03/21] fuel cell example - draft --- examples/37_pem_fc/37_pem_fc.yaml | 5 ++ examples/37_pem_fc/driver_config.yaml | 4 ++ examples/37_pem_fc/hb_inputs.csv | 5 ++ examples/37_pem_fc/plant_config.yaml | 57 ++++++++++++++++++ examples/37_pem_fc/run_pem_fuel_cell.py | 18 ++++++ examples/37_pem_fc/tech_config.yaml | 58 ++++++++++++++++++ .../converters/hydrogen/PEM_h2_fuel_cell.py | 60 ++++++++++++------- h2integrate/core/supported_models.py | 2 + 8 files changed, 186 insertions(+), 23 deletions(-) create mode 100644 examples/37_pem_fc/37_pem_fc.yaml create mode 100644 examples/37_pem_fc/driver_config.yaml create mode 100644 examples/37_pem_fc/hb_inputs.csv create mode 100644 examples/37_pem_fc/plant_config.yaml create mode 100644 examples/37_pem_fc/run_pem_fuel_cell.py create mode 100644 examples/37_pem_fc/tech_config.yaml diff --git a/examples/37_pem_fc/37_pem_fc.yaml b/examples/37_pem_fc/37_pem_fc.yaml new file mode 100644 index 000000000..a2e06f587 --- /dev/null +++ b/examples/37_pem_fc/37_pem_fc.yaml @@ -0,0 +1,5 @@ +name: H2Integrate_config +system_summary: This reference contains a hydrogen fuel cell meeting an electricity demand +driver_config: driver_config.yaml +technology_config: tech_config.yaml +plant_config: plant_config.yaml diff --git a/examples/37_pem_fc/driver_config.yaml b/examples/37_pem_fc/driver_config.yaml new file mode 100644 index 000000000..e6b823fec --- /dev/null +++ b/examples/37_pem_fc/driver_config.yaml @@ -0,0 +1,4 @@ +name: driver_config +description: This analysis runs a hybrid plant to match the first example in H2Integrate +general: + folder_output: outputs diff --git a/examples/37_pem_fc/hb_inputs.csv b/examples/37_pem_fc/hb_inputs.csv new file mode 100644 index 000000000..64d8ca180 --- /dev/null +++ b/examples/37_pem_fc/hb_inputs.csv @@ -0,0 +1,5 @@ +Index 0,Index 1,Index 2,Index 3,Index 4,Index 5,Type,Haber Bosch Big,Haber Bosch Small +technologies,ammonia,model_inputs,shared_parameters,production_capacity,,float,100000,10000 +technologies,electrolyzer,model_inputs,performance_parameters,n_clusters,,int,16,2 +technologies,electrolyzer,model_inputs,performance_parameters,include_degradation_penalty,,bool,TRUE,FALSE +technologies,electrolyzer,model_inputs,financial_parameters,capital_items,replacement_cost_percent,float,0.15,0.25 diff --git a/examples/37_pem_fc/plant_config.yaml b/examples/37_pem_fc/plant_config.yaml new file mode 100644 index 000000000..b4cd89bc6 --- /dev/null +++ b/examples/37_pem_fc/plant_config.yaml @@ -0,0 +1,57 @@ +name: plant_config +description: This plant is located in MN, USA... +sites: + site: + latitude: 32.31714 + longitude: -100.18 +# array of arrays containing left-to-right technology +# interconnections; can support bidirectional connections +# with the reverse definition. +# this will naturally grow as we mature the interconnected tech +technology_interconnections: + # connect feedstocks to fuel cell + - [h2_feedstock, h2_fuel_cell, hydrogen, pipe] + - [o2_feedstock, h2_fuel_cell, oxygen, pipe] + # connect fuel cell to demand component + - [h2_fuel_cell, electricity_load_demand, electricity, cable] +plant: + plant_life: 1 +finance_parameters: + finance_groups: + finance_model: ProFastLCO + model_inputs: + params: + analysis_start_year: 2032 + installation_time: 36 # months + inflation_rate: 0.0 # 0 for real analysis + discount_rate: 0.06 # nominal return based on 2024 ATB baseline workbook for land-based wind + debt_equity_ratio: 0.724 # 2024 ATB uses 72.4% debt for land-based wind + property_tax_and_insurance: 0.025 # percent of CAPEX estimated based on https://www.nlr.gov/docs/fy25osti/91775.pdf https://www.house.mn.gov/hrd/issinfo/clsrates.aspx + total_income_tax_rate: 0.2574 # 0.257 tax rate in 2024 atb baseline workbook, value here is based on federal (21%) and state in MN (9.8) + capital_gains_tax_rate: 0.15 # H2FAST default + sales_tax_rate: 0.0 # average combined state and local sales tax https://taxfoundation.org/location/texas/ + debt_interest_rate: 0.07 # based on 2024 ATB nominal interest rate for land-based wind + debt_type: Revolving debt # can be "Revolving debt" or "One time loan". Revolving debt is H2FAST default and leads to much lower LCOH + loan_period_if_used: 0 # H2FAST default, not used for revolving debt + cash_onhand_months: 1 # H2FAST default + admin_expense: 0.00 # percent of sales H2FAST default + capital_items: + depr_type: MACRS # can be "MACRS" or "Straight line" + depr_period: 7 # 5 years - for clean energy facilities as specified by the IRS MACRS schedule https://www.irs.gov/publications/p946#en_US_2020_publink1000107507 + refurb: [0.] + cost_adjustment_parameters: + cost_year_adjustment_inflation: 0.025 # used to adjust modeled costs to target_dollar_year + target_dollar_year: 2022 + finance_subgroups: + h2: + commodity: hydrogen + commodity_stream: h2_feedstock + technologies: [h2_feedstock] + o2: + commodity: oxygen + commodity_stream: o2_feedstock + technologies: [o2_feedstock] + electricity: + commodity: electricity + technologies: + - h2_fuel_cell diff --git a/examples/37_pem_fc/run_pem_fuel_cell.py b/examples/37_pem_fc/run_pem_fuel_cell.py new file mode 100644 index 000000000..4b86610c8 --- /dev/null +++ b/examples/37_pem_fc/run_pem_fuel_cell.py @@ -0,0 +1,18 @@ +import numpy as np + +from h2integrate import H2IntegrateModel + + +# Create a H2Integrate model +model = H2IntegrateModel("37_pem_fc.yaml") + +# Setup the model +model.setup() + +# Set fuel cell demand profile +demand_profile = np.ones(8760) * 20000 +model.prob.set_val("h2_fuel_cell.electricity_set_point", demand_profile, units="kW") + +# Run model +model.run() +model.post_process() diff --git a/examples/37_pem_fc/tech_config.yaml b/examples/37_pem_fc/tech_config.yaml new file mode 100644 index 000000000..d72a12287 --- /dev/null +++ b/examples/37_pem_fc/tech_config.yaml @@ -0,0 +1,58 @@ +name: technology_config +description: This hybrid plant produces ammonia +technologies: + h2_fuel_cell: + performance_model: + model: PEMH2FuelCellPerformanceModel + cost_model: + model: H2FuelCellCostModel + model_inputs: + shared_parameters: + system_capacity_kw: 30000 + performance_parameters: + n_stacks: 30 + stack_temperature_K: 278 # Not used yet + cost_parameters: + capex_per_kw: 150 + fixed_opex_per_kw_per_year: 35 + cost_year: 2022 + electricity_load_demand: + performance_model: + model: GenericDemandComponent + model_inputs: + performance_parameters: + commodity: electricity + commodity_rate_units: kW + demand_profile: 20000 + h2_feedstock: + performance_model: + model: FeedstockPerformanceModel + cost_model: + model: FeedstockCostModel + model_inputs: + shared_parameters: + commodity: hydrogen + commodity_rate_units: kg/h + performance_parameters: + rated_capacity: 50.0 # kg/h of hydrogen + cost_parameters: + cost_year: 2022 + price: 5.0 + annual_cost: 0. + start_up_cost: 0.0 + o2_feedstock: + performance_model: + model: FeedstockPerformanceModel + cost_model: + model: FeedstockCostModel + model_inputs: + shared_parameters: + commodity: oxygen + commodity_rate_units: kg/h + performance_parameters: + rated_capacity: 50.0 # kg/h of oxygen + cost_parameters: + cost_year: 2022 + price: 0.5 # $5/kg of oxygen. price is extremely variable + annual_cost: 0. + start_up_cost: 0.0 diff --git a/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py b/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py index 07938bc47..3890f3f6e 100644 --- a/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py +++ b/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py @@ -3,7 +3,7 @@ from scipy.interpolate import make_interp_spline from h2integrate.core.utilities import BaseConfig, merge_shared_inputs -from h2integrate.core.validators import gte_zero, range_val +from h2integrate.core.validators import gte_zero from h2integrate.core.model_baseclasses import ( CostModelBaseClass, CostModelBaseConfig, @@ -25,13 +25,13 @@ class PEMH2FuelCellPerformanceConfig(BaseConfig): # How does N_cells translate to electricity rating? system_capacity_kw: float = field(validator=gte_zero) - n_stacks = int + n_stacks: int stack_temperature_K: float - min_system_power_fraction_kw: float - fuel_cell_efficiency_hhv: float = field(validator=range_val(0, 1)) + # min_system_power_fraction_kw: float + # fuel_cell_efficiency_hhv: float = field(validator=range_val(0, 1)) -def calc_current(power_ref, cell_area, stack_number): +def calc_current(power_ref, cell_area, n_cells, stack_number): # Calculates the current and voltage from IV curve based on power reference current_curve = [ 0.0356, @@ -46,10 +46,10 @@ def calc_current(power_ref, cell_area, stack_number): 1.52666667, 1.80333333, 2.07, - 2.32, - 2.54333333, - 2.73666667, - 2.9, + # 2.32, + # 2.54333333, + # 2.73666667, + # 2.9, ] # in A voltage_curve = [ 0.987, @@ -64,10 +64,10 @@ def calc_current(power_ref, cell_area, stack_number): 0.53566667, 0.486, 0.436, - 0.386, - 0.33533333, - 0.286, - 0.236, + # 0.386, + # 0.33533333, + # 0.286, + # 0.236, ] # in V power_curve = [ 35.16666667, @@ -82,15 +82,20 @@ def calc_current(power_ref, cell_area, stack_number): 817.0, 875.33333333, 902.0, - 895.0, - 854.33333333, - 782.66666667, - 684.33333333, - ] / 1e6 # change from mW to kW + # 895.0, + # 854.33333333, + # 782.66666667, + # 684.33333333, + ] + + # Change power from mW to kW + power_curve = [x / 1e6 for x in power_curve] + print(power_curve) power_I_curve = make_interp_spline(power_curve, current_curve, k=3) V_I_curve = make_interp_spline(current_curve, voltage_curve, k=5) - power_density = power_ref / cell_area / stack_number + power_density = power_ref / cell_area / stack_number / n_cells + print("Power density", power_density) I_cell = power_I_curve(power_density) V_cell = V_I_curve(I_cell) @@ -220,7 +225,7 @@ def compute(self, inputs, outputs): # fuel_cell_efficiency = inputs["fuel_cell_efficiency"] # Set calculation constants: - self.f_c = 96485.33 # Faraday's constant in Coulombs/mol + self.f_c = 96485.33 # Faraday's constant in A/mol self.M_H2 = 0.002016 # Molar mass of H2 in kg/mol self.M_O2 = 0.32 # Molar mass of O2 in kg/mol self.Tref = 298.15 # Standard room temperature in K [25 deg Celsius] @@ -234,6 +239,7 @@ def compute(self, inputs, outputs): self.n_cells = 100 # number of cells per stack - how to size this?? # is n_cells = N_series? + self.N_series = 1 self.stack_size = self.config.system_capacity_kw / self.config.n_stacks self.cell_active_area = 1250 # [cm^2] @@ -258,16 +264,24 @@ def compute(self, inputs, outputs): O2in = inputs["oxygen_in"][i] # Find current and voltage from IV curve with power setpoint - I_cell, V_cell = calc_current(power_reference) + I_cell, V_cell = calc_current( + power_reference, self.cell_active_area, self.n_cells, self.config.n_stacks + ) # Calculate hydrogen and oxygen consumed + print("I_cell", I_cell) + print("f_c", self.f_c) + print("M_H2", self.M_H2) + print("self.dt", self.dt) H2_consumed_rate = ( - I_cell * self.n_cells / (2.0 * self.f_c) * self.M_H2 * self.dt + (I_cell * self.N_series) / (2.0 * self.f_c) * self.M_H2 * self.dt ) # kg/time step O2_consumed_rate = ( - I_cell * self.n_cells / (4.0 * self.f_c) * self.M_O2 * self.dt + (I_cell * self.N_series) / (4.0 * self.f_c) * self.M_O2 * self.dt ) # kg/time step + print("H2 and O2 consumed per second", H2_consumed_rate, O2_consumed_rate) + if H2_consumed_rate > H2in or O2_consumed_rate > O2in: print("Not enough H2 or O2 for this power point") # implement an adjustment based on H2 & O2 available diff --git a/h2integrate/core/supported_models.py b/h2integrate/core/supported_models.py index b2abecd78..b858ec2d0 100644 --- a/h2integrate/core/supported_models.py +++ b/h2integrate/core/supported_models.py @@ -76,6 +76,7 @@ from h2integrate.resource.wind.nlr_developer_wtk_api import WTKNLRDeveloperAPIWindResource from h2integrate.converters.hydrogen.basic_cost_model import BasicElectrolyzerCostModel from h2integrate.converters.hydrogen.pem_electrolyzer import ECOElectrolyzerPerformanceModel +from h2integrate.converters.hydrogen.PEM_h2_fuel_cell import PEMH2FuelCellPerformanceModel from h2integrate.converters.solar.atb_res_com_pv_cost import ATBResComPVCostModel from h2integrate.converters.solar.atb_utility_pv_cost import ATBUtilityPVCostModel from h2integrate.converters.iron.martin_mine_cost_model import MartinIronMineCostComponent @@ -224,6 +225,7 @@ "WOMBATElectrolyzerModel": WOMBATElectrolyzerModel, "LinearH2FuelCellPerformanceModel": LinearH2FuelCellPerformanceModel, "H2FuelCellCostModel": H2FuelCellCostModel, + "PEMH2FuelCellPerformanceModel": PEMH2FuelCellPerformanceModel, "SteamMethaneReformerPerformanceModel": SteamMethaneReformerPerformanceModel, "SteamMethaneReformerCostModel": SteamMethaneReformerCostModel, "SimpleASUCostModel": SimpleASUCostModel, From 040cf3bcd2b8c01fd82a141e3c793b4281d7b280 Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Thu, 4 Jun 2026 18:26:48 -0400 Subject: [PATCH 04/21] Update example --- examples/37_pem_fc/plant_config.yaml | 2 +- examples/37_pem_fc/tech_config.yaml | 6 +- .../converters/hydrogen/PEM_h2_fuel_cell.py | 79 +++++++++++-------- 3 files changed, 51 insertions(+), 36 deletions(-) diff --git a/examples/37_pem_fc/plant_config.yaml b/examples/37_pem_fc/plant_config.yaml index b4cd89bc6..b3e0eb53c 100644 --- a/examples/37_pem_fc/plant_config.yaml +++ b/examples/37_pem_fc/plant_config.yaml @@ -15,7 +15,7 @@ technology_interconnections: # connect fuel cell to demand component - [h2_fuel_cell, electricity_load_demand, electricity, cable] plant: - plant_life: 1 + plant_life: 2 finance_parameters: finance_groups: finance_model: ProFastLCO diff --git a/examples/37_pem_fc/tech_config.yaml b/examples/37_pem_fc/tech_config.yaml index d72a12287..66e3bae32 100644 --- a/examples/37_pem_fc/tech_config.yaml +++ b/examples/37_pem_fc/tech_config.yaml @@ -10,7 +10,7 @@ technologies: shared_parameters: system_capacity_kw: 30000 performance_parameters: - n_stacks: 30 + n_stacks: 1200 stack_temperature_K: 278 # Not used yet cost_parameters: capex_per_kw: 150 @@ -34,7 +34,7 @@ technologies: commodity: hydrogen commodity_rate_units: kg/h performance_parameters: - rated_capacity: 50.0 # kg/h of hydrogen + rated_capacity: 10.0 # kg/h of hydrogen cost_parameters: cost_year: 2022 price: 5.0 @@ -50,7 +50,7 @@ technologies: commodity: oxygen commodity_rate_units: kg/h performance_parameters: - rated_capacity: 50.0 # kg/h of oxygen + rated_capacity: 100.0 # kg/h of oxygen cost_parameters: cost_year: 2022 price: 0.5 # $5/kg of oxygen. price is extremely variable diff --git a/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py b/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py index 3890f3f6e..80fc58fe3 100644 --- a/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py +++ b/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py @@ -40,12 +40,12 @@ def calc_current(power_ref, cell_area, n_cells, stack_number): 0.11366667, 0.244, 0.454, - 0.70366667, - 0.96933333, - 1.24, - 1.52666667, - 1.80333333, - 2.07, + # 0.70366667, + # 0.96933333, + # 1.24, + # 1.52666667, + # 1.80333333, + # 2.07, # 2.32, # 2.54333333, # 2.73666667, @@ -58,12 +58,12 @@ def calc_current(power_ref, cell_area, n_cells, stack_number): 0.838, 0.786, 0.736, - 0.686, - 0.636, - 0.586, - 0.53566667, - 0.486, - 0.436, + # 0.686, + # 0.636, + # 0.586, + # 0.53566667, + # 0.486, + # 0.436, # 0.386, # 0.33533333, # 0.286, @@ -76,12 +76,12 @@ def calc_current(power_ref, cell_area, n_cells, stack_number): 95.46666667, 191.66666667, 334.33333333, - 482.66666667, - 616.66666667, - 729.0, - 817.0, - 875.33333333, - 902.0, + # 482.66666667, + # 616.66666667, + # 729.0, + # 817.0, + # 875.33333333, + # 902.0, # 895.0, # 854.33333333, # 782.66666667, @@ -89,16 +89,19 @@ def calc_current(power_ref, cell_area, n_cells, stack_number): ] # Change power from mW to kW - power_curve = [x / 1e6 for x in power_curve] - print(power_curve) + power_curve = [x / 1e3 for x in power_curve] power_I_curve = make_interp_spline(power_curve, current_curve, k=3) V_I_curve = make_interp_spline(current_curve, voltage_curve, k=5) + + # convert power_ref to Watts + power_ref = power_ref * 1e3 power_density = power_ref / cell_area / stack_number / n_cells print("Power density", power_density) I_cell = power_I_curve(power_density) V_cell = V_I_curve(I_cell) + I_cell = I_cell * cell_area return I_cell, V_cell @@ -237,11 +240,14 @@ def compute(self, inputs, outputs): self.hhv_air = 0 # No higher heating value of air self.hhv_water = 0 # ???? Need to check this - self.n_cells = 100 # number of cells per stack - how to size this?? + self.max_cell_power_density = 0.000334 # is n_cells = N_series? self.N_series = 1 self.stack_size = self.config.system_capacity_kw / self.config.n_stacks - self.cell_active_area = 1250 # [cm^2] + self.cell_active_area = 400 # [cm^2] from Battelle (https://www.energy.gov/sites/prod/files/2018/02/f49/fcto_battelle_mfg_cost_analysis_1%20_to_25kw_pp_chp_fc_systems_jan2017_0.pdf) + self.n_cells = round( + self.stack_size / (self.cell_active_area * self.max_cell_power_density) + ) # PSUEDO CODE: """ @@ -269,25 +275,34 @@ def compute(self, inputs, outputs): ) # Calculate hydrogen and oxygen consumed - print("I_cell", I_cell) - print("f_c", self.f_c) - print("M_H2", self.M_H2) - print("self.dt", self.dt) - H2_consumed_rate = ( - (I_cell * self.N_series) / (2.0 * self.f_c) * self.M_H2 * self.dt + # print("I_cell", I_cell) + # print("f_c", self.f_c) + # print("M_H2", self.M_H2) + H2_consumed_rate = ((I_cell * self.N_series * self.M_H2) / (2.0 * self.f_c)) * ( + self.dt * self.config.n_stacks ) # kg/time step - O2_consumed_rate = ( - (I_cell * self.N_series) / (4.0 * self.f_c) * self.M_O2 * self.dt + O2_consumed_rate = ((I_cell * self.N_series * self.M_O2) / (4.0 * self.f_c)) * ( + self.dt * self.config.n_stacks ) # kg/time step - print("H2 and O2 consumed per second", H2_consumed_rate, O2_consumed_rate) + print("H2 and O2 consumed per hour", H2_consumed_rate, O2_consumed_rate) if H2_consumed_rate > H2in or O2_consumed_rate > O2in: print("Not enough H2 or O2 for this power point") # implement an adjustment based on H2 & O2 available # Compute electricity from the system - electricity_produced = V_cell * I_cell * self.n_cells * self.config.n_stacks + electricity_produced = ( + V_cell * I_cell * self.n_cells * self.config.n_stacks / 1e3 + ) # Calculated in watts, convert to kW + print( + "electricity produced:", + V_cell, + I_cell, + self.n_cells, + self.config.n_stacks, + electricity_produced, + ) # Compute H2O out (is this needed?) From 7e2c69364b47a14a18f1fe55fbb0ec1f8c960458 Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Thu, 4 Jun 2026 18:37:19 -0400 Subject: [PATCH 05/21] updates --- examples/37_pem_fc/run_pem_fuel_cell.py | 2 +- examples/37_pem_fc/tech_config.yaml | 6 +++--- h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py | 8 +++++--- 3 files changed, 9 insertions(+), 7 deletions(-) diff --git a/examples/37_pem_fc/run_pem_fuel_cell.py b/examples/37_pem_fc/run_pem_fuel_cell.py index 4b86610c8..6cbc09505 100644 --- a/examples/37_pem_fc/run_pem_fuel_cell.py +++ b/examples/37_pem_fc/run_pem_fuel_cell.py @@ -10,7 +10,7 @@ model.setup() # Set fuel cell demand profile -demand_profile = np.ones(8760) * 20000 +demand_profile = np.ones(8760) * 1000 model.prob.set_val("h2_fuel_cell.electricity_set_point", demand_profile, units="kW") # Run model diff --git a/examples/37_pem_fc/tech_config.yaml b/examples/37_pem_fc/tech_config.yaml index 66e3bae32..86054b95a 100644 --- a/examples/37_pem_fc/tech_config.yaml +++ b/examples/37_pem_fc/tech_config.yaml @@ -8,9 +8,9 @@ technologies: model: H2FuelCellCostModel model_inputs: shared_parameters: - system_capacity_kw: 30000 + system_capacity_kw: 1500 performance_parameters: - n_stacks: 1200 + n_stacks: 60 stack_temperature_K: 278 # Not used yet cost_parameters: capex_per_kw: 150 @@ -23,7 +23,7 @@ technologies: performance_parameters: commodity: electricity commodity_rate_units: kW - demand_profile: 20000 + demand_profile: 1000 h2_feedstock: performance_model: model: FeedstockPerformanceModel diff --git a/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py b/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py index 80fc58fe3..876c0f04a 100644 --- a/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py +++ b/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py @@ -88,7 +88,7 @@ def calc_current(power_ref, cell_area, n_cells, stack_number): # 684.33333333, ] - # Change power from mW to kW + # Change power from mW to W power_curve = [x / 1e3 for x in power_curve] power_I_curve = make_interp_spline(power_curve, current_curve, k=3) @@ -96,6 +96,7 @@ def calc_current(power_ref, cell_area, n_cells, stack_number): # convert power_ref to Watts power_ref = power_ref * 1e3 + print(power_ref, cell_area, stack_number, n_cells) power_density = power_ref / cell_area / stack_number / n_cells print("Power density", power_density) @@ -279,13 +280,14 @@ def compute(self, inputs, outputs): # print("f_c", self.f_c) # print("M_H2", self.M_H2) H2_consumed_rate = ((I_cell * self.N_series * self.M_H2) / (2.0 * self.f_c)) * ( - self.dt * self.config.n_stacks + self.dt * self.config.n_stacks * self.n_cells ) # kg/time step O2_consumed_rate = ((I_cell * self.N_series * self.M_O2) / (4.0 * self.f_c)) * ( - self.dt * self.config.n_stacks + self.dt * self.config.n_stacks * self.n_cells ) # kg/time step print("H2 and O2 consumed per hour", H2_consumed_rate, O2_consumed_rate) + print(self.stack_size, self.n_cells) if H2_consumed_rate > H2in or O2_consumed_rate > O2in: print("Not enough H2 or O2 for this power point") From 03529107b0f3e299c41ec383121f6caa434703ac Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Fri, 5 Jun 2026 13:25:17 -0400 Subject: [PATCH 06/21] Remove print statements --- examples/37_pem_fc/plant_config.yaml | 13 +++---- .../converters/hydrogen/PEM_h2_fuel_cell.py | 34 +++++++++---------- 2 files changed, 24 insertions(+), 23 deletions(-) diff --git a/examples/37_pem_fc/plant_config.yaml b/examples/37_pem_fc/plant_config.yaml index b3e0eb53c..689a9d0f2 100644 --- a/examples/37_pem_fc/plant_config.yaml +++ b/examples/37_pem_fc/plant_config.yaml @@ -8,12 +8,13 @@ sites: # interconnections; can support bidirectional connections # with the reverse definition. # this will naturally grow as we mature the interconnected tech -technology_interconnections: - # connect feedstocks to fuel cell - - [h2_feedstock, h2_fuel_cell, hydrogen, pipe] - - [o2_feedstock, h2_fuel_cell, oxygen, pipe] - # connect fuel cell to demand component - - [h2_fuel_cell, electricity_load_demand, electricity, cable] +# technology_interconnections: [] + # - [] + # # connect feedstocks to fuel cell + # - [h2_feedstock, h2_fuel_cell, hydrogen, pipe] + # - [o2_feedstock, h2_fuel_cell, oxygen, pipe] + # # connect fuel cell to demand component + # - [h2_fuel_cell, electricity_load_demand, electricity, cable] plant: plant_life: 2 finance_parameters: diff --git a/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py b/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py index 876c0f04a..cf4e326e5 100644 --- a/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py +++ b/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py @@ -96,9 +96,9 @@ def calc_current(power_ref, cell_area, n_cells, stack_number): # convert power_ref to Watts power_ref = power_ref * 1e3 - print(power_ref, cell_area, stack_number, n_cells) + # print(power_ref, cell_area, stack_number, n_cells) power_density = power_ref / cell_area / stack_number / n_cells - print("Power density", power_density) + # print("Power density", power_density) I_cell = power_I_curve(power_density) V_cell = V_I_curve(I_cell) @@ -267,8 +267,8 @@ def compute(self, inputs, outputs): for i in range(self.n_timesteps): power_reference = inputs[f"{self.commodity}_set_point"][i] - H2in = inputs["hydrogen_in"][i] - O2in = inputs["oxygen_in"][i] + inputs["hydrogen_in"][i] + inputs["oxygen_in"][i] # Find current and voltage from IV curve with power setpoint I_cell, V_cell = calc_current( @@ -286,25 +286,25 @@ def compute(self, inputs, outputs): self.dt * self.config.n_stacks * self.n_cells ) # kg/time step - print("H2 and O2 consumed per hour", H2_consumed_rate, O2_consumed_rate) - print(self.stack_size, self.n_cells) + # print("H2 and O2 consumed per hour", H2_consumed_rate, O2_consumed_rate) + # print(self.stack_size, self.n_cells) - if H2_consumed_rate > H2in or O2_consumed_rate > O2in: - print("Not enough H2 or O2 for this power point") - # implement an adjustment based on H2 & O2 available + # if H2_consumed_rate > H2in or O2_consumed_rate > O2in: + # print("Not enough H2 or O2 for this power point") + # implement an adjustment based on H2 & O2 available # Compute electricity from the system electricity_produced = ( V_cell * I_cell * self.n_cells * self.config.n_stacks / 1e3 ) # Calculated in watts, convert to kW - print( - "electricity produced:", - V_cell, - I_cell, - self.n_cells, - self.config.n_stacks, - electricity_produced, - ) + # print( + # "electricity produced:", + # V_cell, + # I_cell, + # self.n_cells, + # self.config.n_stacks, + # electricity_produced, + # ) # Compute H2O out (is this needed?) From 6bb575124088ff26c4f8448f4db0e34be2969083 Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Fri, 5 Jun 2026 13:32:53 -0400 Subject: [PATCH 07/21] add new model to hydrogen init file --- examples/37_pem_fc/plant_config.yaml | 13 ++++++------- h2integrate/converters/hydrogen/__init__.py | 3 +++ 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/examples/37_pem_fc/plant_config.yaml b/examples/37_pem_fc/plant_config.yaml index 689a9d0f2..7b0bf2776 100644 --- a/examples/37_pem_fc/plant_config.yaml +++ b/examples/37_pem_fc/plant_config.yaml @@ -8,13 +8,12 @@ sites: # interconnections; can support bidirectional connections # with the reverse definition. # this will naturally grow as we mature the interconnected tech -# technology_interconnections: [] - # - [] - # # connect feedstocks to fuel cell - # - [h2_feedstock, h2_fuel_cell, hydrogen, pipe] - # - [o2_feedstock, h2_fuel_cell, oxygen, pipe] - # # connect fuel cell to demand component - # - [h2_fuel_cell, electricity_load_demand, electricity, cable] +technology_interconnections: + # connect feedstocks to fuel cell + - [h2_feedstock, h2_fuel_cell, hydrogen, pipe] + - [o2_feedstock, h2_fuel_cell, oxygen, pipe] + # connect fuel cell to demand component + - [h2_fuel_cell, electricity_load_demand, electricity, cable] plant: plant_life: 2 finance_parameters: diff --git a/h2integrate/converters/hydrogen/__init__.py b/h2integrate/converters/hydrogen/__init__.py index b39fb2556..5bf3e68c2 100644 --- a/h2integrate/converters/hydrogen/__init__.py +++ b/h2integrate/converters/hydrogen/__init__.py @@ -9,6 +9,9 @@ LinearH2FuelCellPerformanceModel, H2FuelCellCostModel, ) +from h2integrate.converters.hydrogen.PEM_h2_fuel_cell import ( + PEMH2FuelCellPerformanceModel, +) from h2integrate.converters.hydrogen.steam_methane_reformer import ( SteamMethaneReformerPerformanceModel, SteamMethaneReformerCostModel, From bee84780e334a3c8a9aaf06445d62223fba8833f Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Thu, 18 Jun 2026 16:07:54 -0400 Subject: [PATCH 08/21] Updates to model --- .../converters/hydrogen/PEM_h2_fuel_cell.py | 78 +++++++++++++++---- 1 file changed, 65 insertions(+), 13 deletions(-) diff --git a/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py b/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py index cf4e326e5..9e88ecd53 100644 --- a/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py +++ b/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py @@ -231,16 +231,18 @@ def compute(self, inputs, outputs): # Set calculation constants: self.f_c = 96485.33 # Faraday's constant in A/mol self.M_H2 = 0.002016 # Molar mass of H2 in kg/mol - self.M_O2 = 0.32 # Molar mass of O2 in kg/mol + self.M_O2 = 0.032 # Molar mass of O2 in kg/mol self.Tref = 298.15 # Standard room temperature in K [25 deg Celsius] self.cp_H2 = 14300 # Specific heat of H2 in J/(kg*K) self.cp_air = 1005 # Specific heat of air in J/(kg*K) - self.cp_water = 4184 # Specific heat of water in J/(kg*K) + self.cp_H2O = 4184 # Specific heat of water in J/(kg*K) self.cp_N2 = 1040 # Specific heat of nitrogen in J/(kg*K) + self.cp_O2 = 918 # Specific heat of oxygen in J/(kg*K) self.hhv_h2 = 141.8 * 1e6 # Higher heating value of hydrogen in J/kg self.hhv_air = 0 # No higher heating value of air - self.hhv_water = 0 # ???? Need to check this + self.hhv_H2O = 2260 # Higher heating value of water in J/kg + # Sizing the cells self.max_cell_power_density = 0.000334 # is n_cells = N_series? self.N_series = 1 @@ -286,8 +288,8 @@ def compute(self, inputs, outputs): self.dt * self.config.n_stacks * self.n_cells ) # kg/time step - # print("H2 and O2 consumed per hour", H2_consumed_rate, O2_consumed_rate) - # print(self.stack_size, self.n_cells) + print("H2 and O2 consumed per hour", H2_consumed_rate, O2_consumed_rate) + print(self.stack_size, self.n_cells) # if H2_consumed_rate > H2in or O2_consumed_rate > O2in: # print("Not enough H2 or O2 for this power point") @@ -297,14 +299,9 @@ def compute(self, inputs, outputs): electricity_produced = ( V_cell * I_cell * self.n_cells * self.config.n_stacks / 1e3 ) # Calculated in watts, convert to kW - # print( - # "electricity produced:", - # V_cell, - # I_cell, - # self.n_cells, - # self.config.n_stacks, - # electricity_produced, - # ) + + # Need to implement this function + # self.calculate_water_production(H2_consumed_rate, O2_consumed_rate) # Compute H2O out (is this needed?) @@ -373,10 +370,65 @@ def compute(self, inputs, outputs): # ) # outputs["hydrogen_consumed"] = outputs["electricity_out"] * kw_to_kgh_h2 + # ############################################################################## + # # Helper functions for energy balance and water production calculations + # def enthalpy_flow(self, m, cp, T, Tref, h0): # """Mass-specific enthalpy flow: Hdot = m * (cp*(T - Tref) + h0)""" # return m * (cp * (T - Tref) + h0) + # def calculate_water_production(self, h2_consumed_rate, o2_consumed_rate): + # # Calculate water production based on stoichiometry of the reaction + # # 2H2 + O2 --> 2H2O + # # For every 2 moles of H2 consumed, 2 moles of H2O are produced + # # For every 1 mole of O2 consumed, 2 moles of H2O are produced + + # # Calculuate the energy balance of the reaction to find the water production + # # Enthalpy change of the reaction: + # ΔH = (m_H2 * hhv_h2) + (m_O2 * hhv_o2) - (m_H2O * hhv_water) + # # Assuming hhv_o2 and hhv_water are 0, we can simplify to: + # # ΔH = (m_H2 * hhv_h2) - (m_H2O * hhv_water) + # # Since hhv_water is 0, we can further simplify to: + # # ΔH = m_H2 * hhv_h2 + # # The energy released by the reaction is equal to the energy produced by the fuel cell + # which is the power output (electricity produced) plus the heat produced (which we can + # assume is a certain percentage of the energy released by the reaction, say 50% for a PEM + # fuel cell). + # # Therefore, we can calculate the water production based on the energy balance of the + # reaction and the power output of the fuel cell. + # H_H2_in = self.enthalpy_flow(mH2_in, self.cp_H2, T_H2_in, self.Tref, self.hhv_h2) + # H_H2_out = self.enthalpy_flow(mH2_out, self.cp_H2, T_H2_out, self.Tref, self.hhv_h2) + + # H_air_in = self.enthalpy_flow(mO2_in, self.cp_air, T_air_in, self.Tref, 0.0) + # H_air_out = ( + # self.enthalpy_flow(mO2_out, self.cp_O2, T_air_out, self.Tref, 0.0) + + # self.enthalpy_flow(mN2_out, self.cp_N2, T_air_out, self.Tref, 0.0) + # ) # consider just oxygen for a first pass + + # H_H2O_in = self.enthalpy_flow(mH2O_in, + # self.cp_H2O, T_H2O_in, self.Tref, self.hhv_H2O) + # H_H2O_out = self.enthalpy_flow(mH2O_out, + # self.cp_H2O, T_H2O_out, self.Tref, self.hhv_H2O) + + # return (H_H2_in + H_air_in + H_H2O_in + # - H_H2_out - H_air_out - H_H2O_out + # - Wel - Q) + + # ############################################################################## + + # # Convert mass flow rates to molar flow rates + # h2_molar_flow = h2_consumed_rate / self.M_H2 # mol/time step + # o2_molar_flow = o2_consumed_rate / self.M_O2 # mol/time step + + # # Calculate water production based on limiting reactant + # h2o_from_h2 = h2_molar_flow * (self.M_H2O / 1) # kg/time step + # h2o_from_o2 = o2_molar_flow * (self.M_H2O / 0.5) # kg/time step + + # # The actual water produced is the minimum of the two calculations + # water_produced = np.minimum(h2o_from_h2, h2o_from_o2) + + # return water_produced + @define(kw_only=True) class H2FuelCellCostConfig(CostModelBaseConfig): From 30205b3f4d4b42812a9cfe2afec96f754b9ed6df Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Fri, 19 Jun 2026 19:01:18 -0400 Subject: [PATCH 09/21] Add more detailed FC cost model, update for passthrough controller --- examples/37_pem_fc/tech_config.yaml | 15 +- .../system_level/cost_minimization_control.py | 2 + .../system_level/system_level_control_base.py | 3 + .../converters/hydrogen/PEM_h2_fuel_cell.py | 138 ++++++------------ h2integrate/converters/hydrogen/__init__.py | 1 + h2integrate/core/supported_models.py | 1 + 6 files changed, 61 insertions(+), 99 deletions(-) diff --git a/examples/37_pem_fc/tech_config.yaml b/examples/37_pem_fc/tech_config.yaml index 86054b95a..0edd34f1e 100644 --- a/examples/37_pem_fc/tech_config.yaml +++ b/examples/37_pem_fc/tech_config.yaml @@ -5,7 +5,7 @@ technologies: performance_model: model: PEMH2FuelCellPerformanceModel cost_model: - model: H2FuelCellCostModel + model: PEMH2FuelCellCostModel model_inputs: shared_parameters: system_capacity_kw: 1500 @@ -13,9 +13,16 @@ technologies: n_stacks: 60 stack_temperature_K: 278 # Not used yet cost_parameters: - capex_per_kw: 150 - fixed_opex_per_kw_per_year: 35 - cost_year: 2022 + capex_stack_per_kw: 800 + capex_hydrogen_supply_per_kw: 140.67 + capex_air_supply_per_kw: 138.79 + capex_cooling_per_kw: 64.14 + capex_controls_instrumentation_per_kw: 106.52 + capex_electrical_per_kw: 447.83 + capex_assembly_per_kw: 63.51 + capex_additional_labor_per_kw: 137.96 + fixed_opex_per_kw_per_year: 31 + cost_year: 2026 electricity_load_demand: performance_model: model: GenericDemandComponent diff --git a/h2integrate/control/control_strategies/system_level/cost_minimization_control.py b/h2integrate/control/control_strategies/system_level/cost_minimization_control.py index 869976d50..43be86f23 100644 --- a/h2integrate/control/control_strategies/system_level/cost_minimization_control.py +++ b/h2integrate/control/control_strategies/system_level/cost_minimization_control.py @@ -74,11 +74,13 @@ def compute(self, inputs, outputs): remaining = np.maximum(demand, 0.0) marginal_costs = self._compute_marginal_costs(inputs) + print(f"Marginal costs: {marginal_costs}") # TODO: remove debug print # Merit order: sort by mean marginal cost (cheapest first) mean_costs = np.array([mc.mean() for mc in marginal_costs]) dispatch_order = np.argsort(mean_costs) + # jkjk # Initialize all dispatchable set-point outputs to zero for set_point_name in self.dispatchable_set_point_names: outputs[set_point_name] = np.zeros(self.n_timesteps) diff --git a/h2integrate/control/control_strategies/system_level/system_level_control_base.py b/h2integrate/control/control_strategies/system_level/system_level_control_base.py index d1b8c0348..89f974fe1 100644 --- a/h2integrate/control/control_strategies/system_level/system_level_control_base.py +++ b/h2integrate/control/control_strategies/system_level/system_level_control_base.py @@ -641,6 +641,9 @@ def _compute_marginal_costs(self, inputs): marginal_cost = np.zeros(self.n_timesteps) marginal_costs.append(marginal_cost) + print( + f"Marginal cost for {marginal_cost_type}: {marginal_cost}" + ) # TODO: remove debug print return marginal_costs diff --git a/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py b/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py index 9e88ecd53..f0e2a19dc 100644 --- a/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py +++ b/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py @@ -96,11 +96,10 @@ def calc_current(power_ref, cell_area, n_cells, stack_number): # convert power_ref to Watts power_ref = power_ref * 1e3 - # print(power_ref, cell_area, stack_number, n_cells) power_density = power_ref / cell_area / stack_number / n_cells # print("Power density", power_density) - I_cell = power_I_curve(power_density) + I_cell = max(power_I_curve(power_density), 0) V_cell = V_I_curve(I_cell) I_cell = I_cell * cell_area return I_cell, V_cell @@ -123,6 +122,7 @@ class PEMH2FuelCellPerformanceModel(PerformanceModelBaseClass): 3600, 3600, ) # (min, max) time step lengths (in seconds) compatible with this model + _control_classifier = "dispatchable" def initialize(self): super().initialize() @@ -200,13 +200,13 @@ def setup(self): desc="Mass flow rate of water consumed by the fuel cell", ) - # Default the electricity set point input as the rated capacity + # Default the electricity command value input as the rated capacity self.add_input( - f"{self.commodity}_set_point", + f"{self.commodity}_command_value", val=self.config.system_capacity_kw, shape=self.n_timesteps, units=self.commodity_rate_units, - desc="Electricity set point for PEM fuel cell", + desc="Electricity command value for natural gas plant", ) def compute(self, inputs, outputs): @@ -216,7 +216,7 @@ def compute(self, inputs, outputs): Args: inputs: OpenMDAO inputs object containing hydrogen_in, fuel cell - HHV efficiency, electricity_set_point, and system_capacity. + HHV efficiency, electricity_command_value, and system_capacity. outputs: OpenMDAO outputs object for electricity_out, hydrogen_consumed. """ @@ -232,6 +232,7 @@ def compute(self, inputs, outputs): self.f_c = 96485.33 # Faraday's constant in A/mol self.M_H2 = 0.002016 # Molar mass of H2 in kg/mol self.M_O2 = 0.032 # Molar mass of O2 in kg/mol + self.M_H2O = 0.018 # Molar mass of H2O in kg/mol self.Tref = 298.15 # Standard room temperature in K [25 deg Celsius] self.cp_H2 = 14300 # Specific heat of H2 in J/(kg*K) self.cp_air = 1005 # Specific heat of air in J/(kg*K) @@ -265,10 +266,11 @@ def compute(self, inputs, outputs): h2_consumed = np.zeros(self.n_timesteps) o2_consumed = np.zeros(self.n_timesteps) + h2o_generated = np.zeros(self.n_timesteps) commodity_out = np.zeros(self.n_timesteps) for i in range(self.n_timesteps): - power_reference = inputs[f"{self.commodity}_set_point"][i] + power_reference = inputs[f"{self.commodity}_command_value"][i] inputs["hydrogen_in"][i] inputs["oxygen_in"][i] @@ -278,9 +280,6 @@ def compute(self, inputs, outputs): ) # Calculate hydrogen and oxygen consumed - # print("I_cell", I_cell) - # print("f_c", self.f_c) - # print("M_H2", self.M_H2) H2_consumed_rate = ((I_cell * self.N_series * self.M_H2) / (2.0 * self.f_c)) * ( self.dt * self.config.n_stacks * self.n_cells ) # kg/time step @@ -288,10 +287,10 @@ def compute(self, inputs, outputs): self.dt * self.config.n_stacks * self.n_cells ) # kg/time step - print("H2 and O2 consumed per hour", H2_consumed_rate, O2_consumed_rate) - print(self.stack_size, self.n_cells) + # print("H2 and O2 consumed per hour", H2_consumed_rate, O2_consumed_rate) + # print(self.stack_size, self.n_cells) - # if H2_consumed_rate > H2in or O2_consumed_rate > O2in: + # TODO: if H2_consumed_rate > H2in or O2_consumed_rate > O2in: # print("Not enough H2 or O2 for this power point") # implement an adjustment based on H2 & O2 available @@ -300,13 +299,16 @@ def compute(self, inputs, outputs): V_cell * I_cell * self.n_cells * self.config.n_stacks / 1e3 ) # Calculated in watts, convert to kW - # Need to implement this function - # self.calculate_water_production(H2_consumed_rate, O2_consumed_rate) - - # Compute H2O out (is this needed?) + # Compute H2O out + H2O_generated = ( + (I_cell * self.N_series / (2 * self.f_c)) + * self.M_H2O + * (self.dt * self.config.n_stacks * self.n_cells) + ) # in kg/time step h2_consumed[i] = H2_consumed_rate o2_consumed[i] = O2_consumed_rate + h2o_generated[i] = H2O_generated commodity_out[i] = electricity_produced # Set Outputs @@ -324,22 +326,11 @@ def compute(self, inputs, outputs): ) outputs["hydrogen_consumed"] = h2_consumed outputs["oxygen_consumed"] = o2_consumed + outputs["water_out"] = h2o_generated # # conversion factor: kW electricity to kg/h hydrogen, units: (kg/h)/kW # kw_to_kgh_h2 = (3600.0 * 0.001) / (fuel_cell_efficiency * HHV_H2_MJ_PER_KG) - # # TODO: Calculate max H2 and O2 consumption of the stacks - # max_h2_consumption = system_capacity * kw_to_kgh_h2 - - # # electrical set point, saturated at maximum rated system capacity - # electricity_set_point = np.where( - # inputs["electricity_set_point"] > system_capacity, - # system_capacity, - # inputs["electricity_set_point"], - # ) - - # h2_demand = electricity_set_point * kw_to_kgh_h2 - # # available feedstock, saturated at maximum system feedstock consumption # h2_available = np.where( # inputs["hydrogen_in"] > max_h2_consumption, @@ -370,80 +361,30 @@ def compute(self, inputs, outputs): # ) # outputs["hydrogen_consumed"] = outputs["electricity_out"] * kw_to_kgh_h2 - # ############################################################################## - # # Helper functions for energy balance and water production calculations - - # def enthalpy_flow(self, m, cp, T, Tref, h0): - # """Mass-specific enthalpy flow: Hdot = m * (cp*(T - Tref) + h0)""" - # return m * (cp * (T - Tref) + h0) - - # def calculate_water_production(self, h2_consumed_rate, o2_consumed_rate): - # # Calculate water production based on stoichiometry of the reaction - # # 2H2 + O2 --> 2H2O - # # For every 2 moles of H2 consumed, 2 moles of H2O are produced - # # For every 1 mole of O2 consumed, 2 moles of H2O are produced - - # # Calculuate the energy balance of the reaction to find the water production - # # Enthalpy change of the reaction: - # ΔH = (m_H2 * hhv_h2) + (m_O2 * hhv_o2) - (m_H2O * hhv_water) - # # Assuming hhv_o2 and hhv_water are 0, we can simplify to: - # # ΔH = (m_H2 * hhv_h2) - (m_H2O * hhv_water) - # # Since hhv_water is 0, we can further simplify to: - # # ΔH = m_H2 * hhv_h2 - # # The energy released by the reaction is equal to the energy produced by the fuel cell - # which is the power output (electricity produced) plus the heat produced (which we can - # assume is a certain percentage of the energy released by the reaction, say 50% for a PEM - # fuel cell). - # # Therefore, we can calculate the water production based on the energy balance of the - # reaction and the power output of the fuel cell. - # H_H2_in = self.enthalpy_flow(mH2_in, self.cp_H2, T_H2_in, self.Tref, self.hhv_h2) - # H_H2_out = self.enthalpy_flow(mH2_out, self.cp_H2, T_H2_out, self.Tref, self.hhv_h2) - - # H_air_in = self.enthalpy_flow(mO2_in, self.cp_air, T_air_in, self.Tref, 0.0) - # H_air_out = ( - # self.enthalpy_flow(mO2_out, self.cp_O2, T_air_out, self.Tref, 0.0) + - # self.enthalpy_flow(mN2_out, self.cp_N2, T_air_out, self.Tref, 0.0) - # ) # consider just oxygen for a first pass - - # H_H2O_in = self.enthalpy_flow(mH2O_in, - # self.cp_H2O, T_H2O_in, self.Tref, self.hhv_H2O) - # H_H2O_out = self.enthalpy_flow(mH2O_out, - # self.cp_H2O, T_H2O_out, self.Tref, self.hhv_H2O) - - # return (H_H2_in + H_air_in + H_H2O_in - # - H_H2_out - H_air_out - H_H2O_out - # - Wel - Q) - - # ############################################################################## - - # # Convert mass flow rates to molar flow rates - # h2_molar_flow = h2_consumed_rate / self.M_H2 # mol/time step - # o2_molar_flow = o2_consumed_rate / self.M_O2 # mol/time step - - # # Calculate water production based on limiting reactant - # h2o_from_h2 = h2_molar_flow * (self.M_H2O / 1) # kg/time step - # h2o_from_o2 = o2_molar_flow * (self.M_H2O / 0.5) # kg/time step - - # # The actual water produced is the minimum of the two calculations - # water_produced = np.minimum(h2o_from_h2, h2o_from_o2) - - # return water_produced - @define(kw_only=True) -class H2FuelCellCostConfig(CostModelBaseConfig): +class PEMH2FuelCellCostConfig(CostModelBaseConfig): """Configuration class for the hydrogen fuel cell cost model. - Fields include `system_capacity_kw`, `capex_per_kw`, and `fixed_opex_per_kw_per_year`. - The `cost_year` field is inherited from `CostModelBaseConfig`. + Fields include `system_capacity_kw`, `capex_stack_per_kw`, `capex_hydrogen_supply_per_kw`, + `capex_air_supply_per_kw`, `capex_cooling_per_kw`, `capex_controls_instrumentation_per_kw`, + `capex_electrical_per_kw`, `capex_assembly_per_kw`, `capex_additional_labor_per_kw`, + and `fixed_opex_per_kw_per_year`. The `cost_year` field is inherited from `CostModelBaseConfig`. """ system_capacity_kw: float = field(validator=gte_zero) - capex_per_kw: float = field(validator=gte_zero) + capex_stack_per_kw: float = field(validator=gte_zero) + capex_hydrogen_supply_per_kw: float = field(validator=gte_zero) + capex_air_supply_per_kw: float = field(validator=gte_zero) + capex_cooling_per_kw: float = field(validator=gte_zero) + capex_controls_instrumentation_per_kw: float = field(validator=gte_zero) + capex_electrical_per_kw: float = field(validator=gte_zero) + capex_assembly_per_kw: float = field(validator=gte_zero) + capex_additional_labor_per_kw: float = field(validator=gte_zero) fixed_opex_per_kw_per_year: float = field(validator=gte_zero) -class H2FuelCellCostModel(CostModelBaseClass): +class PEMH2FuelCellCostModel(CostModelBaseClass): """ Cost model for a hydrogen fuel cell system. @@ -457,7 +398,7 @@ class H2FuelCellCostModel(CostModelBaseClass): ) # (min, max) time step lengths (in seconds) compatible with this model def setup(self): - self.config = H2FuelCellCostConfig.from_dict( + self.config = PEMH2FuelCellCostConfig.from_dict( merge_shared_inputs(self.options["tech_config"]["model_inputs"], "cost"), additional_cls_name=self.__class__.__name__, ) @@ -473,7 +414,14 @@ def setup(self): self.add_input( "unit_capex", - val=self.config.capex_per_kw, + val=self.config.capex_stack_per_kw + + self.config.capex_hydrogen_supply_per_kw + + self.config.capex_air_supply_per_kw + + self.config.capex_cooling_per_kw + + self.config.capex_controls_instrumentation_per_kw + + self.config.capex_electrical_per_kw + + self.config.capex_assembly_per_kw + + self.config.capex_additional_labor_per_kw, units="USD/kW", desc="Capital cost per unit capacity", ) diff --git a/h2integrate/converters/hydrogen/__init__.py b/h2integrate/converters/hydrogen/__init__.py index 5bf3e68c2..a5be25669 100644 --- a/h2integrate/converters/hydrogen/__init__.py +++ b/h2integrate/converters/hydrogen/__init__.py @@ -11,6 +11,7 @@ ) from h2integrate.converters.hydrogen.PEM_h2_fuel_cell import ( PEMH2FuelCellPerformanceModel, + PEMH2FuelCellCostModel, ) from h2integrate.converters.hydrogen.steam_methane_reformer import ( SteamMethaneReformerPerformanceModel, diff --git a/h2integrate/core/supported_models.py b/h2integrate/core/supported_models.py index fb469c3f2..8f9a78360 100644 --- a/h2integrate/core/supported_models.py +++ b/h2integrate/core/supported_models.py @@ -81,6 +81,7 @@ def copy(self): "WOMBATElectrolyzerModel": "converters.hydrogen:WOMBATElectrolyzerModel", "LinearH2FuelCellPerformanceModel": "converters.hydrogen:LinearH2FuelCellPerformanceModel", "PEMH2FuelCellPerformanceModel": "converters.hydrogen:PEMH2FuelCellPerformanceModel", + "PEMH2FuelCellCostModel": "converters.hydrogen:PEMH2FuelCellCostModel", "H2FuelCellCostModel": "converters.hydrogen:H2FuelCellCostModel", "SteamMethaneReformerPerformanceModel": "converters.hydrogen:SteamMethaneReformerPerformanceModel", "SteamMethaneReformerCostModel": "converters.hydrogen:SteamMethaneReformerCostModel", From 8b47b456382205533f23e7df0f3a6e575e6dac01 Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Wed, 24 Jun 2026 14:36:37 -0400 Subject: [PATCH 10/21] Update interpolation model and add NG FC model --- .../system_level/cost_minimization_control.py | 2 +- .../system_level/system_level_control_base.py | 6 +- .../converters/hydrogen/PEM_h2_fuel_cell.py | 7 +- .../converters/natural_gas/SO_NG_fuel_cell.py | 472 ++++++++++++++++++ .../converters/natural_gas/__init__.py | 4 + h2integrate/core/supported_models.py | 2 + 6 files changed, 486 insertions(+), 7 deletions(-) create mode 100644 h2integrate/converters/natural_gas/SO_NG_fuel_cell.py diff --git a/h2integrate/control/control_strategies/system_level/cost_minimization_control.py b/h2integrate/control/control_strategies/system_level/cost_minimization_control.py index 43be86f23..126e8bd4c 100644 --- a/h2integrate/control/control_strategies/system_level/cost_minimization_control.py +++ b/h2integrate/control/control_strategies/system_level/cost_minimization_control.py @@ -74,7 +74,7 @@ def compute(self, inputs, outputs): remaining = np.maximum(demand, 0.0) marginal_costs = self._compute_marginal_costs(inputs) - print(f"Marginal costs: {marginal_costs}") # TODO: remove debug print + # print(f"Marginal costs: {marginal_costs}") # TODO: remove debug print # Merit order: sort by mean marginal cost (cheapest first) mean_costs = np.array([mc.mean() for mc in marginal_costs]) diff --git a/h2integrate/control/control_strategies/system_level/system_level_control_base.py b/h2integrate/control/control_strategies/system_level/system_level_control_base.py index 89f974fe1..493e50ce1 100644 --- a/h2integrate/control/control_strategies/system_level/system_level_control_base.py +++ b/h2integrate/control/control_strategies/system_level/system_level_control_base.py @@ -637,13 +637,13 @@ def _compute_marginal_costs(self, inputs): marginal_cost = self._varopex_marginal_cost(inputs, marginal_cost_data) elif marginal_cost_type == "feedstock": marginal_cost = self._feedstock_marginal_cost(inputs, marginal_cost_data) + # print( + # f"Marginal cost for {marginal_cost_type}: {marginal_cost}" + # ) # TODO: remove debug print else: marginal_cost = np.zeros(self.n_timesteps) marginal_costs.append(marginal_cost) - print( - f"Marginal cost for {marginal_cost_type}: {marginal_cost}" - ) # TODO: remove debug print return marginal_costs diff --git a/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py b/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py index f0e2a19dc..e5fc269c2 100644 --- a/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py +++ b/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py @@ -1,6 +1,5 @@ import numpy as np from attrs import field, define -from scipy.interpolate import make_interp_spline from h2integrate.core.utilities import BaseConfig, merge_shared_inputs from h2integrate.core.validators import gte_zero @@ -91,8 +90,10 @@ def calc_current(power_ref, cell_area, n_cells, stack_number): # Change power from mW to W power_curve = [x / 1e3 for x in power_curve] - power_I_curve = make_interp_spline(power_curve, current_curve, k=3) - V_I_curve = make_interp_spline(current_curve, voltage_curve, k=5) + power_coefs = np.polyfit(power_curve, current_curve, 5) + power_I_curve = np.poly1d(power_coefs) + V_coefs = np.polyfit(current_curve, voltage_curve, 5) + V_I_curve = np.poly1d(V_coefs) # convert power_ref to Watts power_ref = power_ref * 1e3 diff --git a/h2integrate/converters/natural_gas/SO_NG_fuel_cell.py b/h2integrate/converters/natural_gas/SO_NG_fuel_cell.py new file mode 100644 index 000000000..9bef1819f --- /dev/null +++ b/h2integrate/converters/natural_gas/SO_NG_fuel_cell.py @@ -0,0 +1,472 @@ +import numpy as np +from attrs import field, define + +from h2integrate.core.utilities import BaseConfig, merge_shared_inputs +from h2integrate.core.validators import gte_zero +from h2integrate.core.model_baseclasses import ( + CostModelBaseClass, + CostModelBaseConfig, + PerformanceModelBaseClass, +) + + +@define(kw_only=True) +class SO_NG_FuelCellPerformanceConfig(BaseConfig): + """Configuration class for the hydrogen fuel cell performance model. + + Attributes: + system_capacity_kw (float): The capacity of the fuel cell system in kilowatts (kW). + fuel_cell_efficiency_hhv (float): The higher heating value efficiency of the + fuel cell (0 <= efficiency <= 1). + """ + + # TODO: how to size the fuel cell? N_cells + N_stacks? + # How does N_cells translate to electricity rating? + + system_capacity_kw: float = field(validator=gte_zero) + n_stacks: int + stack_temperature_K: float + # min_system_power_fraction_kw: float + # fuel_cell_efficiency_hhv: float = field(validator=range_val(0, 1)) + + +def calc_current(power_ref, cell_area, n_cells, stack_number): + # Calculates the current and voltage from IV curve based on power reference + current_curve = [ + 0.0356, + 0.05413333, + 0.0796, + 0.11366667, + 0.244, + 0.454, + # 0.70366667, + # 0.96933333, + # 1.24, + # 1.52666667, + # 1.80333333, + # 2.07, + # 2.32, + # 2.54333333, + # 2.73666667, + # 2.9, + ] # in A + voltage_curve = [ + 0.987, + 0.936, + 0.884, + 0.838, + 0.786, + 0.736, + # 0.686, + # 0.636, + # 0.586, + # 0.53566667, + # 0.486, + # 0.436, + # 0.386, + # 0.33533333, + # 0.286, + # 0.236, + ] # in V + power_curve = [ + 35.16666667, + 50.53333333, + 70.33333333, + 95.46666667, + 191.66666667, + 334.33333333, + # 482.66666667, + # 616.66666667, + # 729.0, + # 817.0, + # 875.33333333, + # 902.0, + # 895.0, + # 854.33333333, + # 782.66666667, + # 684.33333333, + ] + + # Change power from mW to W + power_curve = [x / 1e3 for x in power_curve] + + power_coefs = np.polyfit(power_curve, current_curve, 5) + power_I_curve = np.poly1d(power_coefs) + V_coefs = np.polyfit(current_curve, voltage_curve, 5) + V_I_curve = np.poly1d(V_coefs) + + # convert power_ref to Watts + power_ref = power_ref * 1e3 + power_density = power_ref / cell_area / stack_number / n_cells + # print("Power density", power_density) + + I_cell = max(power_I_curve(power_density), 0) + V_cell = V_I_curve(I_cell) + I_cell = I_cell * cell_area + return I_cell, V_cell + + +class SO_NG_FuelCellPerformanceModel(PerformanceModelBaseClass): + """ + Performance model for a hydrogen fuel cell. + + The model implements the relationship: + electricity_out = hydrogen_in * fuel_cell_efficiency_hhv * HHV_hydrogen + + where: + - hydrogen_in is the mass flow rate of hydrogen in kg/hr + - fuel_cell_efficiency is the efficiency of the fuel cell (0 <= efficiency <= 1) + - HHV_hydrogen is the higher heating value of hydrogen (approximately 142 MJ/kg) + """ + + _time_step_bounds = ( + 3600, + 3600, + ) # (min, max) time step lengths (in seconds) compatible with this model + _control_classifier = "dispatchable" + + def initialize(self): + super().initialize() + self.commodity = "electricity" + self.commodity_rate_units = "kW" + self.commodity_amount_units = "kW*h" + + def setup(self): + super().setup() + + self.config = SO_NG_FuelCellPerformanceConfig.from_dict( + merge_shared_inputs(self.options["tech_config"]["model_inputs"], "performance"), + additional_cls_name=self.__class__.__name__, + ) + + # Add natural gas input, default to 0 --> set using feedstock component + # or upstream hydrogen converter component + self.add_input( + "natural_gas_in", + val=0.0, + shape=self.n_timesteps, + units="kg/h", + ) + + self.add_input( + "oxygen_in", + val=0.0, + shape=self.n_timesteps, + units="kg/h", + ) + + self.add_input( + "stack_temperature", + val=self.config.stack_temperature_K, + units="K", + desc="Operating temperature of the stack", + ) + + # self.add_input( + # "fuel_cell_efficiency", + # val=self.config.fuel_cell_efficiency_hhv, + # units=None, + # desc="HHV efficiency of the fuel cell (0 <= efficiency <= 1)", + # ) + + # Add rated capacity as an input with config value as default + self.add_input( + "system_capacity", + val=self.config.system_capacity_kw, + units="kW", + desc="Capacity of the h2 fuel cell system", + ) + + self.add_output( + "natural_gas_consumed", + val=0.0, + shape=self.n_timesteps, + units="kg/h", + desc="Mass flow rate of natural gas consumed by the fuel cell", + ) + + self.add_output( + "oxygen_consumed", + val=0.0, + shape=self.n_timesteps, + units="kg/h", + desc="Mass flow rate of oxygen consumed by the fuel cell", + ) + + self.add_output( + "water_out", + val=0.0, + shape=self.n_timesteps, + units="kg/h", + desc="Mass flow rate of water produced by the fuel cell", + ) + + self.add_output( + "carbon_dioxide_out", + val=0.0, + shape=self.n_timesteps, + units="kg/h", + desc="Mass flow rate of carbon dioxide produced by the fuel cell", + ) + + # Default the electricity command value input as the rated capacity + self.add_input( + f"{self.commodity}_command_value", + val=self.config.system_capacity_kw, + shape=self.n_timesteps, + units=self.commodity_rate_units, + desc="Electricity command value for SOFC plant", + ) + + def compute(self, inputs, outputs): + """ + Compute electricity output from the fuel cell based on hydrogen input + and fuel cell HHV efficiency. + + Args: + inputs: OpenMDAO inputs object containing natural_gas_in, fuel cell + HHV efficiency, electricity_command_value, and system_capacity. + outputs: OpenMDAO outputs object for electricity_out, + natural_gas_consumed. + """ + + # calculate max input and output + inputs["system_capacity"] # plant capacity in kW + inputs["natural_gas_in"] # kg/h + inputs["oxygen_in"] # kg/h + inputs["stack_temperature"] + # fuel_cell_efficiency = inputs["fuel_cell_efficiency"] + + # Set calculation constants: + self.f_c = 96485.33 # Faraday's constant in A/mol + self.M_H2 = 0.002016 # Molar mass of H2 in kg/mol + self.M_O2 = 0.032 # Molar mass of O2 in kg/mol + self.M_H2O = 0.018 # Molar mass of H2O in kg/mol + self.M_CO2 = 0.044 # Molar mass of CO2 in kg/mol + self.M_CH4 = 0.01604 # Molar mass of CH4 (methane) in kg/mol + self.Tref = 298.15 # Standard room temperature in K [25 deg Celsius] + self.cp_H2 = 14300 # Specific heat of H2 in J/(kg*K) + self.cp_air = 1005 # Specific heat of air in J/(kg*K) + self.cp_H2O = 4184 # Specific heat of water in J/(kg*K) + self.cp_N2 = 1040 # Specific heat of nitrogen in J/(kg*K) + self.cp_O2 = 918 # Specific heat of oxygen in J/(kg*K) + self.hhv_h2 = 141.8 * 1e6 # Higher heating value of hydrogen in J/kg + self.hhv_air = 0 # No higher heating value of air + self.hhv_H2O = 2260 # Higher heating value of water in J/kg + + # Sizing the cells + self.max_cell_power_density = 0.000334 + # is n_cells = N_series? + self.N_series = 1 + self.stack_size = self.config.system_capacity_kw / self.config.n_stacks + self.cell_active_area = 400 # [cm^2] from Battelle (https://www.energy.gov/sites/prod/files/2018/02/f49/fcto_battelle_mfg_cost_analysis_1%20_to_25kw_pp_chp_fc_systems_jan2017_0.pdf) + self.n_cells = round( + self.stack_size / (self.cell_active_area * self.max_cell_power_density) + ) + + # PSUEDO CODE: + """ + 1. Receive power setpoint into fuel cell + 2. Find current with I-V curve - NEED THIS + 3. Calculate power out with current - NEED TO FIND CELL ELECTRICITY CALC + 4. Calculate H2 consumed, O2 consumed, water out + 5. See if H2 in and O2 in can provide this + 6. Repeat step 4 if H2 or O2 limit power out + + """ + + ng_consumed = np.zeros(self.n_timesteps) + o2_consumed = np.zeros(self.n_timesteps) + co2_generated = np.zeros(self.n_timesteps) + h2o_generated = np.zeros(self.n_timesteps) + commodity_out = np.zeros(self.n_timesteps) + + for i in range(self.n_timesteps): + power_reference = inputs[f"{self.commodity}_command_value"][i] + inputs["natural_gas_in"][i] + inputs["oxygen_in"][i] + + # Find current and voltage from IV curve with power setpoint + I_cell, V_cell = calc_current( + power_reference, self.cell_active_area, self.n_cells, self.config.n_stacks + ) + + # Calculate natural gas and oxygen consumed + NG_consumed_rate = ((I_cell * self.N_series * self.M_CH4) / (8.0 * self.f_c)) * ( + self.dt * self.config.n_stacks * self.n_cells + ) # kg/time step + O2_consumed_rate = ((I_cell * self.N_series * self.M_O2) / (4.0 * self.f_c)) * ( + self.dt * self.config.n_stacks * self.n_cells + ) # kg/time step + + # print("NG and O2 consumed per hour", NG_consumed_rate, O2_consumed_rate) + # print(self.stack_size, self.n_cells) + + # TODO: if NG_consumed_rate > NGin or O2_consumed_rate > O2in: + # print("Not enough NG or O2 for this power point") + # implement an adjustment based on H2 & O2 available + + # Compute electricity from the system + electricity_produced = ( + V_cell * I_cell * self.n_cells * self.config.n_stacks / 1e3 + ) # Calculated in watts, convert to kW + + # Compute H2O out + H2O_generated = ( + (I_cell * self.N_series / (2 * self.f_c)) + * self.M_H2O + * (self.dt * self.config.n_stacks * self.n_cells) + ) # in kg/time step + # TODO: check electron numbers in these calculations + # Compute CO2 out + CO2_generated = ( + (I_cell * self.N_series / (8 * self.f_c)) + * self.M_CO2 + * (self.dt * self.config.n_stacks * self.n_cells) + ) # in kg/time step + + ng_consumed[i] = NG_consumed_rate + o2_consumed[i] = O2_consumed_rate + h2o_generated[i] = H2O_generated + co2_generated[i] = CO2_generated + commodity_out[i] = electricity_produced + + # Set Outputs + # clip the electricity output to the system capacity + outputs["electricity_out"] = np.minimum(commodity_out, self.config.system_capacity_kw) + outputs["total_electricity_produced"] = np.sum(outputs["electricity_out"]) * ( + self.dt / 3600 + ) + outputs["rated_electricity_production"] = self.config.system_capacity_kw + outputs["annual_electricity_produced"] = outputs["total_electricity_produced"] * ( + 1 / self.fraction_of_year_simulated + ) + outputs["capacity_factor"] = outputs["total_electricity_produced"] / ( + self.config.system_capacity_kw * self.n_timesteps * (self.dt / 3600) + ) + outputs["natural_gas_consumed"] = ng_consumed + outputs["oxygen_consumed"] = o2_consumed + outputs["water_out"] = h2o_generated + outputs["carbon_dioxide_out"] = co2_generated + + # # conversion factor: kW electricity to kg/h hydrogen, units: (kg/h)/kW + # kw_to_kgh_h2 = (3600.0 * 0.001) / (fuel_cell_efficiency * HHV_H2_MJ_PER_KG) + + # # available feedstock, saturated at maximum system feedstock consumption + # h2_available = np.where( + # inputs["hydrogen_in"] > max_h2_consumption, + # max_h2_consumption, + # inputs["hydrogen_in"], + # ) + + # # h2 consumed is minimum between available feedstock and output demand + # hydrogen_in = np.minimum(h2_available, h2_demand) + + # # make any negative hydrogen input zero + # hydrogen_in = np.maximum(hydrogen_in, 0.0) + + # # calculate electricity output in kW + # electricity_out_kw = hydrogen_in / kw_to_kgh_h2 + + # # clip the electricity output to the system capacity + # outputs["electricity_out"] = np.minimum(electricity_out_kw, system_capacity) + # outputs["total_electricity_produced"] = np.sum(outputs["electricity_out"]) * ( + # self.dt / 3600 + # ) + # outputs["rated_electricity_production"] = system_capacity + # outputs["annual_electricity_produced"] = outputs["total_electricity_produced"] * ( + # 1 / self.fraction_of_year_simulated + # ) + # outputs["capacity_factor"] = outputs["total_electricity_produced"] / ( + # system_capacity * self.n_timesteps * (self.dt / 3600) + # ) + # outputs["hydrogen_consumed"] = outputs["electricity_out"] * kw_to_kgh_h2 + + +@define(kw_only=True) +class SO_NG_FuelCellCostConfig(CostModelBaseConfig): + """Configuration class for the hydrogen fuel cell cost model. + + Fields include `system_capacity_kw`, `capex_stack_per_kw`, `capex_hydrogen_supply_per_kw`, + `capex_air_supply_per_kw`, `capex_cooling_per_kw`, `capex_controls_instrumentation_per_kw`, + `capex_electrical_per_kw`, `capex_assembly_per_kw`, `capex_additional_labor_per_kw`, + and `fixed_opex_per_kw_per_year`. The `cost_year` field is inherited from `CostModelBaseConfig`. + """ + + system_capacity_kw: float = field(validator=gte_zero) + capex_stack_per_kw: float = field(validator=gte_zero) + capex_hydrogen_supply_per_kw: float = field(validator=gte_zero) + capex_air_supply_per_kw: float = field(validator=gte_zero) + capex_cooling_per_kw: float = field(validator=gte_zero) + capex_controls_instrumentation_per_kw: float = field(validator=gte_zero) + capex_electrical_per_kw: float = field(validator=gte_zero) + capex_assembly_per_kw: float = field(validator=gte_zero) + capex_additional_labor_per_kw: float = field(validator=gte_zero) + fixed_opex_per_kw_per_year: float = field(validator=gte_zero) + + +class SO_NG_FuelCellCostModel(CostModelBaseClass): + """ + Cost model for a hydrogen fuel cell system. + + The model calculates capital and fixed operating costs based on system capacity and + specified cost parameters. + """ + + _time_step_bounds = ( + 3600, + 3600, + ) # (min, max) time step lengths (in seconds) compatible with this model + + def setup(self): + self.config = SO_NG_FuelCellCostConfig.from_dict( + merge_shared_inputs(self.options["tech_config"]["model_inputs"], "cost"), + additional_cls_name=self.__class__.__name__, + ) + + super().setup() + + self.add_input( + "system_capacity", + val=self.config.system_capacity_kw, + units="kW", + desc="Capacity of the h2 fuel cell system", + ) + + self.add_input( + "unit_capex", + val=self.config.capex_stack_per_kw + + self.config.capex_hydrogen_supply_per_kw + + self.config.capex_air_supply_per_kw + + self.config.capex_cooling_per_kw + + self.config.capex_controls_instrumentation_per_kw + + self.config.capex_electrical_per_kw + + self.config.capex_assembly_per_kw + + self.config.capex_additional_labor_per_kw, + units="USD/kW", + desc="Capital cost per unit capacity", + ) + + self.add_input( + "fixed_opex_per_year", + val=self.config.fixed_opex_per_kw_per_year, + units="(USD/kW)/year", + desc="Fixed operating expenses per unit capacity per year", + ) + + def compute(self, inputs, outputs, discrete_inputs, discrete_outputs): + """ + Compute capital and fixed operating costs for the fuel cell system. + + Args: + inputs: OpenMDAO inputs object containing system_capacity. + outputs: OpenMDAO outputs object for capital_cost and fixed_operating_cost_per_year. + """ + + system_capacity_kw = inputs["system_capacity"] + + # Calculate capital cost + outputs["CapEx"] = system_capacity_kw * inputs["unit_capex"] + + # Calculate fixed operating cost per year + outputs["OpEx"] = system_capacity_kw * inputs["fixed_opex_per_year"] diff --git a/h2integrate/converters/natural_gas/__init__.py b/h2integrate/converters/natural_gas/__init__.py index 50d2251a8..2050acff8 100644 --- a/h2integrate/converters/natural_gas/__init__.py +++ b/h2integrate/converters/natural_gas/__init__.py @@ -8,3 +8,7 @@ SimpleGasConsumerPerformance, SimpleGasConsumerCost, ) +from h2integrate.converters.natural_gas.SO_NG_fuel_cell import ( + SO_NG_FuelCellPerformanceModel, + SO_NG_FuelCellCostModel, +) diff --git a/h2integrate/core/supported_models.py b/h2integrate/core/supported_models.py index 8f9a78360..610efc5ca 100644 --- a/h2integrate/core/supported_models.py +++ b/h2integrate/core/supported_models.py @@ -85,6 +85,8 @@ def copy(self): "H2FuelCellCostModel": "converters.hydrogen:H2FuelCellCostModel", "SteamMethaneReformerPerformanceModel": "converters.hydrogen:SteamMethaneReformerPerformanceModel", "SteamMethaneReformerCostModel": "converters.hydrogen:SteamMethaneReformerCostModel", + "SO_NG_FuelCellPerformanceModel": "converters.natural_gas:SO_NG_FuelCellPerformanceModel", + "SO_NG_FuelCellCostModel": "converters.natural_gas:SO_NG_FuelCellCostModel", "SimpleASUCostModel": "converters.nitrogen:SimpleASUCostModel", "SimpleASUPerformanceModel": "converters.nitrogen:SimpleASUPerformanceModel", "HOPPComponent": "converters.hopp:HOPPComponent", From e10050b238626d7dcb83319c702b1e22b4bfe26f Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Wed, 24 Jun 2026 16:24:33 -0400 Subject: [PATCH 11/21] Update doc strings --- .../converters/hydrogen/PEM_h2_fuel_cell.py | 123 +++++++++--------- .../converters/natural_gas/SO_NG_fuel_cell.py | 92 +++++-------- 2 files changed, 92 insertions(+), 123 deletions(-) diff --git a/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py b/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py index e5fc269c2..7f37495fe 100644 --- a/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py +++ b/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py @@ -16,8 +16,8 @@ class PEMH2FuelCellPerformanceConfig(BaseConfig): Attributes: system_capacity_kw (float): The capacity of the fuel cell system in kilowatts (kW). - fuel_cell_efficiency_hhv (float): The higher heating value efficiency of the - fuel cell (0 <= efficiency <= 1). + n_stacks (int): The number of stacks in the fuel cell system. + stack_temperature_K (float): The operating temperature of the fuel cell stack in Kelvin (K). """ # TODO: how to size the fuel cell? N_cells + N_stacks? @@ -108,15 +108,25 @@ def calc_current(power_ref, cell_area, n_cells, stack_number): class PEMH2FuelCellPerformanceModel(PerformanceModelBaseClass): """ - Performance model for a hydrogen fuel cell. - - The model implements the relationship: - electricity_out = hydrogen_in * fuel_cell_efficiency_hhv * HHV_hydrogen - - where: - - hydrogen_in is the mass flow rate of hydrogen in kg/hr - - fuel_cell_efficiency is the efficiency of the fuel cell (0 <= efficiency <= 1) - - HHV_hydrogen is the higher heating value of hydrogen (approximately 142 MJ/kg) + Performance model for a PEM hydrogen fuel cell. + + The model simulates electrochemical conversion of hydrogen and oxygen into electricity + and water. It calculates: + - hydrogen and oxygen consumption based on electrochemical reactions + - water production as a byproduct + - electricity output based on system capacity and operational conditions + + Inputs: + - hydrogen_in: mass flow rate of hydrogen (kg/h) + - oxygen_in: mass flow rate of oxygen (kg/h) + - stack_temperature: operating temperature of the fuel cell stack (K) + - system_capacity: rated capacity of the fuel cell system (kW) + + Outputs: + - hydrogen_consumed: hydrogen consumption rate (kg/h) + - oxygen_consumed: oxygen consumption rate (kg/h) + - water_out: water production rate (kg/h) + - electricity_out: electricity output (kW) """ _time_step_bounds = ( @@ -139,8 +149,6 @@ def setup(self): additional_cls_name=self.__class__.__name__, ) - # Add natural gas input, default to 0 --> set using feedstock component - # or upstream hydrogen converter component self.add_input( "hydrogen_in", val=0.0, @@ -198,7 +206,7 @@ def setup(self): val=0.0, shape=self.n_timesteps, units="kg/h", - desc="Mass flow rate of water consumed by the fuel cell", + desc="Mass flow rate of water produced by the fuel cell", ) # Default the electricity command value input as the rated capacity @@ -207,19 +215,22 @@ def setup(self): val=self.config.system_capacity_kw, shape=self.n_timesteps, units=self.commodity_rate_units, - desc="Electricity command value for natural gas plant", + desc="Electricity command value for PEM fuel cell", ) def compute(self, inputs, outputs): """ - Compute electricity output from the fuel cell based on hydrogen input - and fuel cell HHV efficiency. + Compute electricity output from the fuel cell based on hydrogen and oxygen input. + + Uses I-V curve characteristics to calculate fuel cell current and voltage, + then computes hydrogen consumed, oxygen consumed, and water generated for + each timestep based on electrochemical reactions. Args: - inputs: OpenMDAO inputs object containing hydrogen_in, fuel cell - HHV efficiency, electricity_command_value, and system_capacity. - outputs: OpenMDAO outputs object for electricity_out, - hydrogen_consumed. + inputs: OpenMDAO inputs object containing hydrogen_in, oxygen_in, + stack_temperature, electricity_command_value, and system_capacity. + outputs: OpenMDAO outputs object for electricity_out, hydrogen_consumed, + oxygen_consumed, water_out, and various electricity production quantities. """ # calculate max input and output @@ -257,11 +268,12 @@ def compute(self, inputs, outputs): # PSUEDO CODE: """ 1. Receive power setpoint into fuel cell - 2. Find current with I-V curve - NEED THIS - 3. Calculate power out with current - NEED TO FIND CELL ELECTRICITY CALC - 4. Calculate H2 consumed, O2 consumed, water out - 5. See if H2 in and O2 in can provide this - 6. Repeat step 4 if H2 or O2 limit power out + 2. Find current with I-V curve + 3. Calculate H2 consumed and O2 consumed + 4. Check if provided H2 and O2 can meet the demand + 5. If not, adjust current + 6. Calculate power out with current and voltage + 7. Calculate the water produced from the reaction """ @@ -272,8 +284,8 @@ def compute(self, inputs, outputs): for i in range(self.n_timesteps): power_reference = inputs[f"{self.commodity}_command_value"][i] - inputs["hydrogen_in"][i] - inputs["oxygen_in"][i] + H2in = inputs["hydrogen_in"][i] + O2in = inputs["oxygen_in"][i] # Find current and voltage from IV curve with power setpoint I_cell, V_cell = calc_current( @@ -291,9 +303,24 @@ def compute(self, inputs, outputs): # print("H2 and O2 consumed per hour", H2_consumed_rate, O2_consumed_rate) # print(self.stack_size, self.n_cells) - # TODO: if H2_consumed_rate > H2in or O2_consumed_rate > O2in: - # print("Not enough H2 or O2 for this power point") - # implement an adjustment based on H2 & O2 available + # TODO: + if H2_consumed_rate > H2in or O2_consumed_rate > O2in: + # implement an adjustment based on H2 & O2 available + new_i_h2 = ( + H2in + / (self.dt * self.config.n_stacks * self.n_cells) + * (2.0 * self.f_c) + / (self.N_series * self.M_H2) + ) + new_i_o2 = ( + O2in + / (self.dt * self.config.n_stacks * self.n_cells) + * (4.0 * self.f_c) + / (self.N_series * self.M_O2) + ) + I_cell = min(new_i_h2, new_i_o2) + # TODO: recalc voltage based on new current + print("Not enough H2 or O2 for this power point") # Compute electricity from the system electricity_produced = ( @@ -329,38 +356,8 @@ def compute(self, inputs, outputs): outputs["oxygen_consumed"] = o2_consumed outputs["water_out"] = h2o_generated - # # conversion factor: kW electricity to kg/h hydrogen, units: (kg/h)/kW - # kw_to_kgh_h2 = (3600.0 * 0.001) / (fuel_cell_efficiency * HHV_H2_MJ_PER_KG) - - # # available feedstock, saturated at maximum system feedstock consumption - # h2_available = np.where( - # inputs["hydrogen_in"] > max_h2_consumption, - # max_h2_consumption, - # inputs["hydrogen_in"], - # ) - - # # h2 consumed is minimum between available feedstock and output demand - # hydrogen_in = np.minimum(h2_available, h2_demand) - - # # make any negative hydrogen input zero - # hydrogen_in = np.maximum(hydrogen_in, 0.0) - - # # calculate electricity output in kW - # electricity_out_kw = hydrogen_in / kw_to_kgh_h2 - - # # clip the electricity output to the system capacity - # outputs["electricity_out"] = np.minimum(electricity_out_kw, system_capacity) - # outputs["total_electricity_produced"] = np.sum(outputs["electricity_out"]) * ( - # self.dt / 3600 - # ) - # outputs["rated_electricity_production"] = system_capacity - # outputs["annual_electricity_produced"] = outputs["total_electricity_produced"] * ( - # 1 / self.fraction_of_year_simulated - # ) - # outputs["capacity_factor"] = outputs["total_electricity_produced"] / ( - # system_capacity * self.n_timesteps * (self.dt / 3600) - # ) - # outputs["hydrogen_consumed"] = outputs["electricity_out"] * kw_to_kgh_h2 + # TODO: implement a hydrogen and oxygen conversion efficiency based on stack + # temperature and other factors @define(kw_only=True) diff --git a/h2integrate/converters/natural_gas/SO_NG_fuel_cell.py b/h2integrate/converters/natural_gas/SO_NG_fuel_cell.py index 9bef1819f..f9145c405 100644 --- a/h2integrate/converters/natural_gas/SO_NG_fuel_cell.py +++ b/h2integrate/converters/natural_gas/SO_NG_fuel_cell.py @@ -12,12 +12,12 @@ @define(kw_only=True) class SO_NG_FuelCellPerformanceConfig(BaseConfig): - """Configuration class for the hydrogen fuel cell performance model. + """Configuration class for the solid oxide natural gas fuel cell performance model. Attributes: system_capacity_kw (float): The capacity of the fuel cell system in kilowatts (kW). - fuel_cell_efficiency_hhv (float): The higher heating value efficiency of the - fuel cell (0 <= efficiency <= 1). + n_stacks (int): The number of stacks in the fuel cell system. + stack_temperature_K (float): The operating temperature of the fuel cell stack in Kelvin (K). """ # TODO: how to size the fuel cell? N_cells + N_stacks? @@ -108,15 +108,17 @@ def calc_current(power_ref, cell_area, n_cells, stack_number): class SO_NG_FuelCellPerformanceModel(PerformanceModelBaseClass): """ - Performance model for a hydrogen fuel cell. + Performance model for a solid oxide natural gas fuel cell. - The model implements the relationship: - electricity_out = hydrogen_in * fuel_cell_efficiency_hhv * HHV_hydrogen + The model calculates electricity output based on natural gas and oxygen inputs, + with current and voltage determined from power density using IV curves. + Produces water and carbon dioxide as byproducts. where: - - hydrogen_in is the mass flow rate of hydrogen in kg/hr - - fuel_cell_efficiency is the efficiency of the fuel cell (0 <= efficiency <= 1) - - HHV_hydrogen is the higher heating value of hydrogen (approximately 142 MJ/kg) + - natural_gas_in is the mass flow rate of natural gas in kg/hr + - oxygen_in is the mass flow rate of oxygen in kg/hr + - water_out is the mass flow rate of water produced in kg/hr + - carbon_dioxide_out is the mass flow rate of carbon dioxide produced in kg/hr """ _time_step_bounds = ( @@ -220,14 +222,14 @@ def setup(self): def compute(self, inputs, outputs): """ - Compute electricity output from the fuel cell based on hydrogen input - and fuel cell HHV efficiency. + Compute electricity output from the SOFC based on natural gas input, + oxygen availability, and fuel cell electrochemical reactions. Args: - inputs: OpenMDAO inputs object containing natural_gas_in, fuel cell - HHV efficiency, electricity_command_value, and system_capacity. - outputs: OpenMDAO outputs object for electricity_out, - natural_gas_consumed. + inputs: OpenMDAO inputs object containing natural_gas_in, oxygen_in, + stack_temperature, electricity_command_value, and system_capacity. + outputs: OpenMDAO outputs object for electricity_out, natural_gas_consumed, + oxygen_consumed, water_out, and carbon_dioxide_out. """ # calculate max input and output @@ -253,7 +255,7 @@ def compute(self, inputs, outputs): self.hhv_h2 = 141.8 * 1e6 # Higher heating value of hydrogen in J/kg self.hhv_air = 0 # No higher heating value of air self.hhv_H2O = 2260 # Higher heating value of water in J/kg - + self.hhv_CO2 = 0 # No higher heating value of CO2 # Sizing the cells self.max_cell_power_density = 0.000334 # is n_cells = N_series? @@ -267,11 +269,11 @@ def compute(self, inputs, outputs): # PSUEDO CODE: """ 1. Receive power setpoint into fuel cell - 2. Find current with I-V curve - NEED THIS - 3. Calculate power out with current - NEED TO FIND CELL ELECTRICITY CALC - 4. Calculate H2 consumed, O2 consumed, water out - 5. See if H2 in and O2 in can provide this - 6. Repeat step 4 if H2 or O2 limit power out + 2. Find current and voltage from I-V curve based on power setpoint + 3. Calculate natural gas (CH4) and oxygen consumed based on Faraday's law + 4. Calculate electricity produced from cell voltage and current + 5. Calculate H2O and CO2 generated from the reaction + 6. Output electricity clipped to system capacity and other metrics """ @@ -349,45 +351,15 @@ def compute(self, inputs, outputs): outputs["water_out"] = h2o_generated outputs["carbon_dioxide_out"] = co2_generated - # # conversion factor: kW electricity to kg/h hydrogen, units: (kg/h)/kW - # kw_to_kgh_h2 = (3600.0 * 0.001) / (fuel_cell_efficiency * HHV_H2_MJ_PER_KG) - - # # available feedstock, saturated at maximum system feedstock consumption - # h2_available = np.where( - # inputs["hydrogen_in"] > max_h2_consumption, - # max_h2_consumption, - # inputs["hydrogen_in"], - # ) - - # # h2 consumed is minimum between available feedstock and output demand - # hydrogen_in = np.minimum(h2_available, h2_demand) - - # # make any negative hydrogen input zero - # hydrogen_in = np.maximum(hydrogen_in, 0.0) - - # # calculate electricity output in kW - # electricity_out_kw = hydrogen_in / kw_to_kgh_h2 - - # # clip the electricity output to the system capacity - # outputs["electricity_out"] = np.minimum(electricity_out_kw, system_capacity) - # outputs["total_electricity_produced"] = np.sum(outputs["electricity_out"]) * ( - # self.dt / 3600 - # ) - # outputs["rated_electricity_production"] = system_capacity - # outputs["annual_electricity_produced"] = outputs["total_electricity_produced"] * ( - # 1 / self.fraction_of_year_simulated - # ) - # outputs["capacity_factor"] = outputs["total_electricity_produced"] / ( - # system_capacity * self.n_timesteps * (self.dt / 3600) - # ) - # outputs["hydrogen_consumed"] = outputs["electricity_out"] * kw_to_kgh_h2 + # TODO: implement a natural gas and oxygen conversion efficiency based on stack + # temperature and other factors @define(kw_only=True) class SO_NG_FuelCellCostConfig(CostModelBaseConfig): - """Configuration class for the hydrogen fuel cell cost model. + """Configuration class for the solid oxide natural gas fuel cell cost model. - Fields include `system_capacity_kw`, `capex_stack_per_kw`, `capex_hydrogen_supply_per_kw`, + Fields include `system_capacity_kw`, `capex_stack_per_kw`, `capex_fuel_supply_per_kw`, `capex_air_supply_per_kw`, `capex_cooling_per_kw`, `capex_controls_instrumentation_per_kw`, `capex_electrical_per_kw`, `capex_assembly_per_kw`, `capex_additional_labor_per_kw`, and `fixed_opex_per_kw_per_year`. The `cost_year` field is inherited from `CostModelBaseConfig`. @@ -395,7 +367,7 @@ class SO_NG_FuelCellCostConfig(CostModelBaseConfig): system_capacity_kw: float = field(validator=gte_zero) capex_stack_per_kw: float = field(validator=gte_zero) - capex_hydrogen_supply_per_kw: float = field(validator=gte_zero) + capex_fuel_supply_per_kw: float = field(validator=gte_zero) capex_air_supply_per_kw: float = field(validator=gte_zero) capex_cooling_per_kw: float = field(validator=gte_zero) capex_controls_instrumentation_per_kw: float = field(validator=gte_zero) @@ -407,7 +379,7 @@ class SO_NG_FuelCellCostConfig(CostModelBaseConfig): class SO_NG_FuelCellCostModel(CostModelBaseClass): """ - Cost model for a hydrogen fuel cell system. + Cost model for a solid oxide natural gas fuel cell system. The model calculates capital and fixed operating costs based on system capacity and specified cost parameters. @@ -430,13 +402,13 @@ def setup(self): "system_capacity", val=self.config.system_capacity_kw, units="kW", - desc="Capacity of the h2 fuel cell system", + desc="Capacity of the solid oxide natural gas fuel cell system", ) self.add_input( "unit_capex", val=self.config.capex_stack_per_kw - + self.config.capex_hydrogen_supply_per_kw + + self.config.capex_fuel_supply_per_kw + self.config.capex_air_supply_per_kw + self.config.capex_cooling_per_kw + self.config.capex_controls_instrumentation_per_kw @@ -456,7 +428,7 @@ def setup(self): def compute(self, inputs, outputs, discrete_inputs, discrete_outputs): """ - Compute capital and fixed operating costs for the fuel cell system. + Compute capital and fixed operating costs for the solid oxide natural gas fuel cell system. Args: inputs: OpenMDAO inputs object containing system_capacity. From fcb5c96fa050bf12607fe11d55ab5404a16f5c85 Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Wed, 24 Jun 2026 16:32:21 -0400 Subject: [PATCH 12/21] Remove print statements from slc --- .../system_level/cost_minimization_control.py | 2 -- .../system_level/system_level_control_base.py | 3 --- 2 files changed, 5 deletions(-) diff --git a/h2integrate/control/control_strategies/system_level/cost_minimization_control.py b/h2integrate/control/control_strategies/system_level/cost_minimization_control.py index 126e8bd4c..869976d50 100644 --- a/h2integrate/control/control_strategies/system_level/cost_minimization_control.py +++ b/h2integrate/control/control_strategies/system_level/cost_minimization_control.py @@ -74,13 +74,11 @@ def compute(self, inputs, outputs): remaining = np.maximum(demand, 0.0) marginal_costs = self._compute_marginal_costs(inputs) - # print(f"Marginal costs: {marginal_costs}") # TODO: remove debug print # Merit order: sort by mean marginal cost (cheapest first) mean_costs = np.array([mc.mean() for mc in marginal_costs]) dispatch_order = np.argsort(mean_costs) - # jkjk # Initialize all dispatchable set-point outputs to zero for set_point_name in self.dispatchable_set_point_names: outputs[set_point_name] = np.zeros(self.n_timesteps) diff --git a/h2integrate/control/control_strategies/system_level/system_level_control_base.py b/h2integrate/control/control_strategies/system_level/system_level_control_base.py index 493e50ce1..d1b8c0348 100644 --- a/h2integrate/control/control_strategies/system_level/system_level_control_base.py +++ b/h2integrate/control/control_strategies/system_level/system_level_control_base.py @@ -637,9 +637,6 @@ def _compute_marginal_costs(self, inputs): marginal_cost = self._varopex_marginal_cost(inputs, marginal_cost_data) elif marginal_cost_type == "feedstock": marginal_cost = self._feedstock_marginal_cost(inputs, marginal_cost_data) - # print( - # f"Marginal cost for {marginal_cost_type}: {marginal_cost}" - # ) # TODO: remove debug print else: marginal_cost = np.zeros(self.n_timesteps) From 1aa8d934ab124b21d0e234c02b9d4528efdc4448 Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Wed, 24 Jun 2026 16:35:29 -0400 Subject: [PATCH 13/21] Update naming to align with naming conventions --- .../converters/natural_gas/SO_NG_fuel_cell.py | 12 ++++++------ h2integrate/converters/natural_gas/__init__.py | 4 ++-- h2integrate/core/supported_models.py | 4 ++-- 3 files changed, 10 insertions(+), 10 deletions(-) diff --git a/h2integrate/converters/natural_gas/SO_NG_fuel_cell.py b/h2integrate/converters/natural_gas/SO_NG_fuel_cell.py index f9145c405..689c955be 100644 --- a/h2integrate/converters/natural_gas/SO_NG_fuel_cell.py +++ b/h2integrate/converters/natural_gas/SO_NG_fuel_cell.py @@ -11,7 +11,7 @@ @define(kw_only=True) -class SO_NG_FuelCellPerformanceConfig(BaseConfig): +class SONGFuelCellPerformanceConfig(BaseConfig): """Configuration class for the solid oxide natural gas fuel cell performance model. Attributes: @@ -106,7 +106,7 @@ def calc_current(power_ref, cell_area, n_cells, stack_number): return I_cell, V_cell -class SO_NG_FuelCellPerformanceModel(PerformanceModelBaseClass): +class SONGFuelCellPerformanceModel(PerformanceModelBaseClass): """ Performance model for a solid oxide natural gas fuel cell. @@ -136,7 +136,7 @@ def initialize(self): def setup(self): super().setup() - self.config = SO_NG_FuelCellPerformanceConfig.from_dict( + self.config = SONGFuelCellPerformanceConfig.from_dict( merge_shared_inputs(self.options["tech_config"]["model_inputs"], "performance"), additional_cls_name=self.__class__.__name__, ) @@ -356,7 +356,7 @@ def compute(self, inputs, outputs): @define(kw_only=True) -class SO_NG_FuelCellCostConfig(CostModelBaseConfig): +class SONGFuelCellCostConfig(CostModelBaseConfig): """Configuration class for the solid oxide natural gas fuel cell cost model. Fields include `system_capacity_kw`, `capex_stack_per_kw`, `capex_fuel_supply_per_kw`, @@ -377,7 +377,7 @@ class SO_NG_FuelCellCostConfig(CostModelBaseConfig): fixed_opex_per_kw_per_year: float = field(validator=gte_zero) -class SO_NG_FuelCellCostModel(CostModelBaseClass): +class SONGFuelCellCostModel(CostModelBaseClass): """ Cost model for a solid oxide natural gas fuel cell system. @@ -391,7 +391,7 @@ class SO_NG_FuelCellCostModel(CostModelBaseClass): ) # (min, max) time step lengths (in seconds) compatible with this model def setup(self): - self.config = SO_NG_FuelCellCostConfig.from_dict( + self.config = SONGFuelCellCostConfig.from_dict( merge_shared_inputs(self.options["tech_config"]["model_inputs"], "cost"), additional_cls_name=self.__class__.__name__, ) diff --git a/h2integrate/converters/natural_gas/__init__.py b/h2integrate/converters/natural_gas/__init__.py index 2050acff8..d349717f8 100644 --- a/h2integrate/converters/natural_gas/__init__.py +++ b/h2integrate/converters/natural_gas/__init__.py @@ -9,6 +9,6 @@ SimpleGasConsumerCost, ) from h2integrate.converters.natural_gas.SO_NG_fuel_cell import ( - SO_NG_FuelCellPerformanceModel, - SO_NG_FuelCellCostModel, + SONGFuelCellPerformanceModel, + SONGFuelCellCostModel, ) diff --git a/h2integrate/core/supported_models.py b/h2integrate/core/supported_models.py index 610efc5ca..3398e5303 100644 --- a/h2integrate/core/supported_models.py +++ b/h2integrate/core/supported_models.py @@ -85,8 +85,8 @@ def copy(self): "H2FuelCellCostModel": "converters.hydrogen:H2FuelCellCostModel", "SteamMethaneReformerPerformanceModel": "converters.hydrogen:SteamMethaneReformerPerformanceModel", "SteamMethaneReformerCostModel": "converters.hydrogen:SteamMethaneReformerCostModel", - "SO_NG_FuelCellPerformanceModel": "converters.natural_gas:SO_NG_FuelCellPerformanceModel", - "SO_NG_FuelCellCostModel": "converters.natural_gas:SO_NG_FuelCellCostModel", + "SONGFuelCellPerformanceModel": "converters.natural_gas:SONGFuelCellPerformanceModel", + "SONGFuelCellCostModel": "converters.natural_gas:SONGFuelCellCostModel", "SimpleASUCostModel": "converters.nitrogen:SimpleASUCostModel", "SimpleASUPerformanceModel": "converters.nitrogen:SimpleASUPerformanceModel", "HOPPComponent": "converters.hopp:HOPPComponent", From 737fe3514034671c5ee9ea3804a1214fcba3b1d6 Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Wed, 24 Jun 2026 16:37:54 -0400 Subject: [PATCH 14/21] remove leftover example input --- examples/37_pem_fc/hb_inputs.csv | 5 ----- 1 file changed, 5 deletions(-) delete mode 100644 examples/37_pem_fc/hb_inputs.csv diff --git a/examples/37_pem_fc/hb_inputs.csv b/examples/37_pem_fc/hb_inputs.csv deleted file mode 100644 index 64d8ca180..000000000 --- a/examples/37_pem_fc/hb_inputs.csv +++ /dev/null @@ -1,5 +0,0 @@ -Index 0,Index 1,Index 2,Index 3,Index 4,Index 5,Type,Haber Bosch Big,Haber Bosch Small -technologies,ammonia,model_inputs,shared_parameters,production_capacity,,float,100000,10000 -technologies,electrolyzer,model_inputs,performance_parameters,n_clusters,,int,16,2 -technologies,electrolyzer,model_inputs,performance_parameters,include_degradation_penalty,,bool,TRUE,FALSE -technologies,electrolyzer,model_inputs,financial_parameters,capital_items,replacement_cost_percent,float,0.15,0.25 From f1ebbd46c65631d66d3fd0f933b0ae83116a16af Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Mon, 29 Jun 2026 18:09:16 -0400 Subject: [PATCH 15/21] Add documentation and changelog --- CHANGELOG.md | 1 + docs/_toc.yml | 2 + docs/technology_models/PEM_h2_fuel_cell.md | 49 +++++++++++++++++++ docs/technology_models/SO_NG_fuel_cell.md | 49 +++++++++++++++++++ examples/37_pem_fc/tech_config.yaml | 4 +- .../converters/natural_gas/SO_NG_fuel_cell.py | 2 + 6 files changed, 105 insertions(+), 2 deletions(-) create mode 100644 docs/technology_models/PEM_h2_fuel_cell.md create mode 100644 docs/technology_models/SO_NG_fuel_cell.md diff --git a/CHANGELOG.md b/CHANGELOG.md index 01151283a..887c4b0a5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -40,6 +40,7 @@ - Added dynamic operating constraints (turndown, ramping, warm/cold start delays) to `AmmoniaSynLoopPerformanceModel` and split `AmmoniaSynLoopCostModel` into its own module. [PR 770](https://github.com/NatLabRockies/H2Integrate/pull/770) - Speed up the slowest tests in the suite by swapping the Floris wind model for `PYSAMWindPlantPerformanceModel` in examples 01 (`01_onshore_steel_mn`) and 02 (`02_texas_ammonia`), updating the affected `test_steel_example`/`test_simple_ammonia_example` expected values, fixing a pre-existing `cases.sql` cache-path bug and module-scoping the fixtures in `h2integrate/postprocess/test/test_sql_timeseries_to_csv.py` so the example only runs once for all four tests. [PR 782](https://github.com/NatLabRockies/H2Integrate/pull/782) - Exposed `n_timesteps`, `dt`, `plant_life`, and `fraction_of_year_simulated` as attributes on `CostModelBaseClass` (matching `PerformanceModelBaseClass`) and updated all cost and performance model subclasses across `h2integrate/` to use these attributes instead of reading them from `plant_config`, removing redundant boilerplate from individual components. [PR 783](https://github.com/NatLabRockies/H2Integrate/pull/783) +- Added two new fuel cell models: `PEMH2FuelCellPerformanceModel` and `PEMH2FuelCellCostModel` to model a PEM hydrogen fuel cell and `SONGFuelCellPerformanceModel` and `SONGFuelCellCostModel` to model a natural gas solid oxide fuel cell [PR 794](https://github.com/NatLabRockies/H2Integrate/pull/794) ## 0.8 [April 15, 2026] diff --git a/docs/_toc.yml b/docs/_toc.yml index 5a47b5529..2e2d92cbc 100644 --- a/docs/_toc.yml +++ b/docs/_toc.yml @@ -49,6 +49,8 @@ parts: - file: technology_models/iron_dri.md - file: technology_models/steel_eaf.md - file: technology_models/steel_eaf_cmu.md + - file: technology_models/PEM_fuel_cell.md + - file: technology_models/SO_NG_fuel_cell.md - caption: Storage Models chapters: - file: storage/storage_models_index.md diff --git a/docs/technology_models/PEM_h2_fuel_cell.md b/docs/technology_models/PEM_h2_fuel_cell.md new file mode 100644 index 000000000..8f182e2c4 --- /dev/null +++ b/docs/technology_models/PEM_h2_fuel_cell.md @@ -0,0 +1,49 @@ +# PEM Hydrogen Fuel Cell Model + +The PEM hydrogen fuel cell performance model implemented in H2Integrate is an electrochemical model that simulates the conversion of hydrogen and oxygen into electricity and water. The model uses a polynomial fit of an I-V (current-voltage) curve to determine the operating point of each cell based on the input `electricity_set_point`, then computes `hydrogen_consumed`, `oxygen_consumed`, `water_out`, and `electricity_out` as the outputs of the system for each timestep. + +The model is sized by `system_capacity_kw` and `n_stacks`. The number of cells per stack is calculated from the stack size, a fixed cell active area (400 cm²)([1](https://www.energy.gov/sites/prod/files/2018/02/f49/fcto_battelle_mfg_cost_analysis_1%20_to_25kw_pp_chp_fc_systems_jan2017_0.pdf)), and a maximum cell power density. Hydrogen and oxygen availability is checked against the demand at each timestep. If the available oxygen or hydrogen is insuffucient for the requested power, the operating current is reduced to match the lowest available feedstock. Electricity output is clipped to the system capacity. + +There are no non-linear operational considerations in this model such as warm-up delays, degraded performance over operational life, voltage recalculation after current adjustment for limited feedstock supply, or thermal dynamics beyond a constant stack temperature input. + +The PEM hydrogen fuel cell cost model is implemented to use cost values in dollars per kilowatt (or per kilowatt per year), decomposed into capital cost components for the stack, hydrogen supply, air supply, cooling, controls and instrumentation, electrical, assembly, and additional labor. + +```{note} +The I-V curve is currently internally defined in the model and not adjustable. +``` + +## Performance Model + +```{eval-rst} +.. autoclass:: h2integrate.converters.hydrogen.PEM_h2_fuel_cell.PEMH2FuelCellPerformanceConfig + :members: + :undoc-members: + :show-inheritance: + :no-index: +``` + +```{eval-rst} +.. autoclass:: h2integrate.converters.hydrogen.PEM_h2_fuel_cell.PEMH2FuelCellPerformanceModel + :members: + :undoc-members: + :show-inheritance: + :no-index: +``` + +## Cost Model + +```{eval-rst} +.. autoclass:: h2integrate.converters.hydrogen.PEM_h2_fuel_cell.PEMH2FuelCellCostConfig + :members: + :undoc-members: + :show-inheritance: + :no-index: +``` + +```{eval-rst} +.. autoclass:: h2integrate.converters.hydrogen.PEM_h2_fuel_cell.PEMH2FuelCellCostModel + :members: + :undoc-members: + :show-inheritance: + :no-index: +``` diff --git a/docs/technology_models/SO_NG_fuel_cell.md b/docs/technology_models/SO_NG_fuel_cell.md new file mode 100644 index 000000000..c2709d4bc --- /dev/null +++ b/docs/technology_models/SO_NG_fuel_cell.md @@ -0,0 +1,49 @@ +# Solid Oxide Natural Gas Fuel Cell Model + +The solid oxide natural gas (SO NG) fuel cell performance model implemented in H2Integrate is an electrochemical model that simulates the conversion of natural gas (assumed to be methane in the chemical reaction equations), steam, and oxygen into electricity, water, and carbon dioxide. The model uses a polynomial fit of an I-V (current-voltage) curve to determine the operating point of each cell based on the input `electricity_set_point`, then computes `natural_gas_consumed`, `oxygen_consumed`, `water_out`, `carbon_dioxide_out`, and `electricity_out` as the outputs of the system for each timestep. + +The model is sized by `system_capacity_kw` and `n_stacks`. The number of cells per stack is calculated from the stack size, a fixed cell active area (400 cm²), and a maximum cell power density. Natural gas and oxygen consumption are computed from the cell current via Faraday's law assuming complete electrochemical oxidation of methane. The model uses the methane-reforming equation rather than a direct electrochemical methan-combusting reaction. This assumes that the methane is broken down in a steam-reforming process before the feedstock enteres the fuel cell ([1](https://www.sciencedirect.com/science/article/pii/S0360319906005726)). Electricity output is clipped to the system capacity. + +There are no non-linear operational considerations in this model such as warm-up delays, degraded performance over operational life, or thermal dynamics beyond a constant stack temperature input. + +The SO NG fuel cell cost model is implemented to use cost values in dollars per kilowatt (or per kilowatt per year), decomposed into capital cost components for the stack, fuel supply, air supply, cooling, controls and instrumentation, electrical, assembly, and additional labor. + +```{note} +The I-V curve is currently internally defined in the model and not adjustable. +``` + +## Performance Model + +```{eval-rst} +.. autoclass:: h2integrate.converters.natural_gas.SO_NG_fuel_cell.SONGFuelCellPerformanceConfig + :members: + :undoc-members: + :show-inheritance: + :no-index: +``` + +```{eval-rst} +.. autoclass:: h2integrate.converters.natural_gas.SO_NG_fuel_cell.SONGFuelCellPerformanceModel + :members: + :undoc-members: + :show-inheritance: + :no-index: +``` + +## Cost Model + +```{eval-rst} +.. autoclass:: h2integrate.converters.natural_gas.SO_NG_fuel_cell.SONGFuelCellCostConfig + :members: + :undoc-members: + :show-inheritance: + :no-index: +``` + +```{eval-rst} +.. autoclass:: h2integrate.converters.natural_gas.SO_NG_fuel_cell.SONGFuelCellCostModel + :members: + :undoc-members: + :show-inheritance: + :no-index: +``` diff --git a/examples/37_pem_fc/tech_config.yaml b/examples/37_pem_fc/tech_config.yaml index 0edd34f1e..93403c7de 100644 --- a/examples/37_pem_fc/tech_config.yaml +++ b/examples/37_pem_fc/tech_config.yaml @@ -41,7 +41,7 @@ technologies: commodity: hydrogen commodity_rate_units: kg/h performance_parameters: - rated_capacity: 10.0 # kg/h of hydrogen + rated_capacity: 50.0 # kg/h of hydrogen cost_parameters: cost_year: 2022 price: 5.0 @@ -57,7 +57,7 @@ technologies: commodity: oxygen commodity_rate_units: kg/h performance_parameters: - rated_capacity: 100.0 # kg/h of oxygen + rated_capacity: 400.0 # kg/h of oxygen cost_parameters: cost_year: 2022 price: 0.5 # $5/kg of oxygen. price is extremely variable diff --git a/h2integrate/converters/natural_gas/SO_NG_fuel_cell.py b/h2integrate/converters/natural_gas/SO_NG_fuel_cell.py index 689c955be..3b13176f4 100644 --- a/h2integrate/converters/natural_gas/SO_NG_fuel_cell.py +++ b/h2integrate/converters/natural_gas/SO_NG_fuel_cell.py @@ -239,6 +239,8 @@ def compute(self, inputs, outputs): inputs["stack_temperature"] # fuel_cell_efficiency = inputs["fuel_cell_efficiency"] + # Add consumption of water for steam reforming of natural gas to hydrogen + # Set calculation constants: self.f_c = 96485.33 # Faraday's constant in A/mol self.M_H2 = 0.002016 # Molar mass of H2 in kg/mol From 0442a8601701f474a393e10ea032785cfaf3cd68 Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Mon, 29 Jun 2026 18:27:11 -0400 Subject: [PATCH 16/21] Add tests for PEM fuel cell --- .../converters/hydrogen/PEM_h2_fuel_cell.py | 7 + .../hydrogen/test/test_PEM_fuel_cell.py | 270 ++++++++++++++++++ 2 files changed, 277 insertions(+) create mode 100644 h2integrate/converters/hydrogen/test/test_PEM_fuel_cell.py diff --git a/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py b/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py index 7f37495fe..25a01eb8a 100644 --- a/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py +++ b/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py @@ -321,6 +321,13 @@ def compute(self, inputs, outputs): I_cell = min(new_i_h2, new_i_o2) # TODO: recalc voltage based on new current print("Not enough H2 or O2 for this power point") + # Calculate hydrogen and oxygen consumed + H2_consumed_rate = ((I_cell * self.N_series * self.M_H2) / (2.0 * self.f_c)) * ( + self.dt * self.config.n_stacks * self.n_cells + ) # kg/time step + O2_consumed_rate = ((I_cell * self.N_series * self.M_O2) / (4.0 * self.f_c)) * ( + self.dt * self.config.n_stacks * self.n_cells + ) # kg/time step # Compute electricity from the system electricity_produced = ( diff --git a/h2integrate/converters/hydrogen/test/test_PEM_fuel_cell.py b/h2integrate/converters/hydrogen/test/test_PEM_fuel_cell.py new file mode 100644 index 000000000..4f2ffc68f --- /dev/null +++ b/h2integrate/converters/hydrogen/test/test_PEM_fuel_cell.py @@ -0,0 +1,270 @@ +import numpy as np +import pytest +import openmdao.api as om +from pytest import fixture + +from h2integrate.converters.hydrogen.PEM_h2_fuel_cell import ( + PEMH2FuelCellCostModel, + PEMH2FuelCellPerformanceModel, +) + + +@fixture +def plant_config(): + plant_config = { + "plant": { + "plant_life": 30, + "simulation": { + "n_timesteps": 8760, + "dt": 3600, + }, + }, + } + return plant_config + + +@fixture +def tech_config(): + config = { + "model_inputs": { + "performance_parameters": { + "system_capacity_kw": 1500.0, + "n_stacks": 60, + "stack_temperature_K": 278.0, + } + } + } + return config + + +@fixture +def cost_config(): + config = { + "model_inputs": { + "cost_parameters": { + "system_capacity_kw": 1500.0, + "capex_stack_per_kw": 800.0, + "capex_hydrogen_supply_per_kw": 140.67, + "capex_air_supply_per_kw": 138.79, + "capex_cooling_per_kw": 64.14, + "capex_controls_instrumentation_per_kw": 106.52, + "capex_electrical_per_kw": 447.83, + "capex_assembly_per_kw": 63.51, + "capex_additional_labor_per_kw": 137.96, + "fixed_opex_per_kw_per_year": 31.0, + "cost_year": 2026, + } + } + } + return config + + +@pytest.mark.regression +def test_fuel_cell_performance(tech_config, plant_config, subtests): + n_timesteps = int(plant_config["plant"]["simulation"]["n_timesteps"]) + + prob = om.Problem() + + fuel_cell = PEMH2FuelCellPerformanceModel( + plant_config=plant_config, tech_config=tech_config, driver_config={} + ) + + prob.model.add_subsystem("fuel_cell", fuel_cell, promotes=["*"]) + + prob.setup() + + # Provide ample hydrogen and oxygen to run at the default command (rated capacity) + hydrogen_input = np.ones(n_timesteps) * 200.0 # kg/h + oxygen_input = np.ones(n_timesteps) * 2000.0 # kg/h + + prob.set_val("fuel_cell.hydrogen_in", hydrogen_input, units="kg/h") + prob.set_val("fuel_cell.oxygen_in", oxygen_input, units="kg/h") + + prob.run_model() + + electricity_output = prob.get_val("fuel_cell.electricity_out", units="kW") + hydrogen_consumed = prob.get_val("fuel_cell.hydrogen_consumed", units="kg/h") + oxygen_consumed = prob.get_val("fuel_cell.oxygen_consumed", units="kg/h") + water_out = prob.get_val("fuel_cell.water_out", units="kg/h") + + with subtests.test("max electricity output bounded by system capacity"): + assert np.max(electricity_output) <= 1500.0 + 1e-6 + + with subtests.test("electricity output is non-negative"): + assert np.min(electricity_output) >= 0.0 + + with subtests.test("rated_electricity_production"): + assert ( + pytest.approx( + prob.get_val("fuel_cell.rated_electricity_production", units="kW"), rel=1e-6 + ) + == 1500.0 + ) + + with subtests.test("total_electricity_produced matches sum of output"): + assert pytest.approx( + prob.get_val("fuel_cell.total_electricity_produced", units="kW*h"), rel=1e-6 + ) == np.sum(electricity_output) + + with subtests.test("annual_electricity_produced equals total for full-year sim"): + assert pytest.approx( + prob.get_val("fuel_cell.annual_electricity_produced", units="kW*h/year"), rel=1e-6 + ) == np.sum(electricity_output) + + with subtests.test("capacity_factor matches definition"): + assert pytest.approx( + prob.get_val("fuel_cell.capacity_factor", units="unitless"), rel=1e-6 + ) == np.sum(electricity_output) / (1500.0 * n_timesteps) + + with subtests.test("hydrogen consumed is non-negative and bounded by supply"): + assert np.min(hydrogen_consumed) >= 0.0 + assert np.max(hydrogen_consumed) <= 200.0 + 1e-6 + + with subtests.test("oxygen consumed is non-negative and bounded by supply"): + assert np.min(oxygen_consumed) >= 0.0 + assert np.max(oxygen_consumed) <= 2000.0 + 1e-6 + + with subtests.test("water produced is non-negative"): + assert np.min(water_out) >= 0.0 + + with subtests.test("stoichiometric mass balance H2 + O2 -> H2O"): + # Combined reactant mass should equal water produced (mass conservation) + np.testing.assert_allclose( + hydrogen_consumed + oxygen_consumed, + water_out, + rtol=1e-3, + ) + + +@pytest.mark.unit +def test_fuel_cell_demand(tech_config, plant_config, subtests): + n_timesteps = int(plant_config["plant"]["simulation"]["n_timesteps"]) + + prob = om.Problem() + + fuel_cell = PEMH2FuelCellPerformanceModel( + plant_config=plant_config, tech_config=tech_config, driver_config={} + ) + + prob.model.add_subsystem("fuel_cell", fuel_cell, promotes=["*"]) + + prob.setup() + + # Provide ample feedstock for most timesteps; constrain a few to test edge cases + hydrogen_input = np.ones(n_timesteps) * 200.0 # kg/h + oxygen_input = np.ones(n_timesteps) * 2000.0 # kg/h + + # Edge cases for feedstock supply at the first 6 timesteps + hydrogen_input[:6] = ( + 500000000.0, # very high H2 supply + 500000000.0, # very high H2 supply with low set point + 200.0, # ample + 0.0, # zero hydrogen supply + 1.0, # very limited hydrogen supply + 200.0, # ample + ) + oxygen_input[:6] = ( + 2000.0, + 2000.0, + 2000.0, + 2000.0, + 2000.0, + 0.0, # zero oxygen supply + ) + + prob.set_val("fuel_cell.hydrogen_in", hydrogen_input, units="kg/h") + prob.set_val("fuel_cell.oxygen_in", oxygen_input, units="kg/h") + + elec_set_point = np.ones(n_timesteps) * 1500.0 # kW + elec_set_point[:6] = ( + 1500.0, # set point equal to system capacity + 500.0, # set point below system capacity + 1500.0, # set point equal to system capacity + 1500.0, # set point equal to system capacity, no H2 supply + 1500.0, # set point equal to system capacity, limited H2 supply + 0.0, # zero set point + ) + + prob.set_val("fuel_cell.electricity_command_value", elec_set_point, units="kW") + + prob.run_model() + + electricity_output = prob.get_val("fuel_cell.electricity_out", units="kW") + hydrogen_consumed = prob.get_val("fuel_cell.hydrogen_consumed", units="kg/h") + oxygen_consumed = prob.get_val("fuel_cell.oxygen_consumed", units="kg/h") + water_out = prob.get_val("fuel_cell.water_out", units="kg/h") + + with subtests.test("output clipped to system capacity"): + assert electricity_output[0] == pytest.approx(1500.0, rel=1e-4) + + with subtests.test("output follows reduced set point"): + # When set point is below capacity and feedstock is ample, output tracks set point + assert electricity_output[1] == pytest.approx(500.0, rel=1e-2) + + with subtests.test("output non-negative when ample supply"): + assert electricity_output[2] >= 0.0 + assert electricity_output[2] <= 1500.0 + 1e-6 + + with subtests.test("zero hydrogen feedstock supply yields zero output"): + assert electricity_output[3] == pytest.approx(0.0, abs=1e-6) + + with subtests.test("limited hydrogen feedstock supply yields reduced output"): + assert electricity_output[4] < 20.0 + assert electricity_output[4] > 0.0 + + with subtests.test("zero set point yields zero output"): + assert electricity_output[5] == pytest.approx(0.0, abs=1e-6) + + # Test hydrogen_consumed, oxygen_consumed, and water_out for the first 6 timesteps + with subtests.test("hydrogen consumed"): + expected_h2_consumed = [76.589328, 22.920501, 76.589328, 0, 1, 0.0] + np.testing.assert_allclose(hydrogen_consumed[:6], expected_h2_consumed, rtol=1e-4) + + with subtests.test("oxygen consumed"): + expected_o2_consumed = [607.851812, 181.908739, 607.851812, 0, 7.936508, 0.0] + np.testing.assert_allclose(oxygen_consumed[:6], expected_o2_consumed, rtol=1e-4) + + with subtests.test("water out"): + expected_water_out = [683.833289, 204.647332, 683.833289, 0.0, 8.928571, 0.0] + np.testing.assert_allclose(water_out[:6], expected_water_out, rtol=1e-4) + + +@pytest.mark.regression +def test_fuel_cell_cost(cost_config, plant_config, subtests): + prob = om.Problem() + + fuel_cell_cost = PEMH2FuelCellCostModel( + plant_config=plant_config, tech_config=cost_config, driver_config={} + ) + + prob.model.add_subsystem("fuel_cell_cost", fuel_cell_cost, promotes=["*"]) + + prob.setup() + + prob.run_model() + + cp = cost_config["model_inputs"]["cost_parameters"] + expected_unit_capex = ( + cp["capex_stack_per_kw"] + + cp["capex_hydrogen_supply_per_kw"] + + cp["capex_air_supply_per_kw"] + + cp["capex_cooling_per_kw"] + + cp["capex_controls_instrumentation_per_kw"] + + cp["capex_electrical_per_kw"] + + cp["capex_assembly_per_kw"] + + cp["capex_additional_labor_per_kw"] + ) + expected_capex = cp["system_capacity_kw"] * expected_unit_capex + expected_opex = cp["system_capacity_kw"] * cp["fixed_opex_per_kw_per_year"] + + with subtests.test("capex value"): + assert ( + pytest.approx(prob.get_val("fuel_cell_cost.CapEx", units="USD"), rel=1e-6) + == expected_capex + ) + + with subtests.test("opex value"): + assert ( + pytest.approx(prob.get_val("fuel_cell_cost.OpEx", units="USD/year"), rel=1e-6) + == expected_opex + ) From 63ac54284046db4dbefd163b0e5730e6929aec5f Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Tue, 30 Jun 2026 18:06:14 -0400 Subject: [PATCH 17/21] Add SONG tests and model updates --- examples/37_pem_fc/plant_config.yaml | 11 +- examples/37_pem_fc/run_pem_fuel_cell.py | 2 +- examples/37_pem_fc/tech_config.yaml | 2 +- examples/38_song_fc/38_song_fc.yaml | 5 + examples/38_song_fc/driver_config.yaml | 4 + examples/38_song_fc/plant_config.yaml | 60 ++++ examples/38_song_fc/run_song_fuel_cell.py | 18 + examples/38_song_fc/tech_config.yaml | 65 ++++ .../converters/natural_gas/SO_NG_fuel_cell.py | 64 +++- .../natural_gas/test/test_SO_NG_fuel_cell.py | 314 ++++++++++++++++++ 10 files changed, 523 insertions(+), 22 deletions(-) create mode 100644 examples/38_song_fc/38_song_fc.yaml create mode 100644 examples/38_song_fc/driver_config.yaml create mode 100644 examples/38_song_fc/plant_config.yaml create mode 100644 examples/38_song_fc/run_song_fuel_cell.py create mode 100644 examples/38_song_fc/tech_config.yaml create mode 100644 h2integrate/converters/natural_gas/test/test_SO_NG_fuel_cell.py diff --git a/examples/37_pem_fc/plant_config.yaml b/examples/37_pem_fc/plant_config.yaml index 7b0bf2776..1237e0198 100644 --- a/examples/37_pem_fc/plant_config.yaml +++ b/examples/37_pem_fc/plant_config.yaml @@ -10,10 +10,10 @@ sites: # this will naturally grow as we mature the interconnected tech technology_interconnections: # connect feedstocks to fuel cell - - [h2_feedstock, h2_fuel_cell, hydrogen, pipe] - - [o2_feedstock, h2_fuel_cell, oxygen, pipe] + - [h2_feedstock, PEM_fuel_cell, hydrogen, pipe] + - [o2_feedstock, PEM_fuel_cell, oxygen, pipe] # connect fuel cell to demand component - - [h2_fuel_cell, electricity_load_demand, electricity, cable] + - [PEM_fuel_cell, electricity_load_demand, electricity, cable] plant: plant_life: 2 finance_parameters: @@ -53,5 +53,8 @@ finance_parameters: technologies: [o2_feedstock] electricity: commodity: electricity + commodity_stream: PEM_fuel_cell technologies: - - h2_fuel_cell + - PEM_fuel_cell + - h2_feedstock + - o2_feedstock diff --git a/examples/37_pem_fc/run_pem_fuel_cell.py b/examples/37_pem_fc/run_pem_fuel_cell.py index 6cbc09505..daa9ce6f3 100644 --- a/examples/37_pem_fc/run_pem_fuel_cell.py +++ b/examples/37_pem_fc/run_pem_fuel_cell.py @@ -11,7 +11,7 @@ # Set fuel cell demand profile demand_profile = np.ones(8760) * 1000 -model.prob.set_val("h2_fuel_cell.electricity_set_point", demand_profile, units="kW") +model.prob.set_val("PEM_fuel_cell.electricity_set_point", demand_profile, units="kW") # Run model model.run() diff --git a/examples/37_pem_fc/tech_config.yaml b/examples/37_pem_fc/tech_config.yaml index 93403c7de..955059111 100644 --- a/examples/37_pem_fc/tech_config.yaml +++ b/examples/37_pem_fc/tech_config.yaml @@ -1,7 +1,7 @@ name: technology_config description: This hybrid plant produces ammonia technologies: - h2_fuel_cell: + PEM_fuel_cell: performance_model: model: PEMH2FuelCellPerformanceModel cost_model: diff --git a/examples/38_song_fc/38_song_fc.yaml b/examples/38_song_fc/38_song_fc.yaml new file mode 100644 index 000000000..a2e06f587 --- /dev/null +++ b/examples/38_song_fc/38_song_fc.yaml @@ -0,0 +1,5 @@ +name: H2Integrate_config +system_summary: This reference contains a hydrogen fuel cell meeting an electricity demand +driver_config: driver_config.yaml +technology_config: tech_config.yaml +plant_config: plant_config.yaml diff --git a/examples/38_song_fc/driver_config.yaml b/examples/38_song_fc/driver_config.yaml new file mode 100644 index 000000000..e6b823fec --- /dev/null +++ b/examples/38_song_fc/driver_config.yaml @@ -0,0 +1,4 @@ +name: driver_config +description: This analysis runs a hybrid plant to match the first example in H2Integrate +general: + folder_output: outputs diff --git a/examples/38_song_fc/plant_config.yaml b/examples/38_song_fc/plant_config.yaml new file mode 100644 index 000000000..45bf0db2a --- /dev/null +++ b/examples/38_song_fc/plant_config.yaml @@ -0,0 +1,60 @@ +name: plant_config +description: This plant is located in MN, USA... +sites: + site: + latitude: 32.31714 + longitude: -100.18 +# array of arrays containing left-to-right technology +# interconnections; can support bidirectional connections +# with the reverse definition. +# this will naturally grow as we mature the interconnected tech +technology_interconnections: + # connect feedstocks to fuel cell + - [ng_feedstock, SONG_fuel_cell, natural_gas, pipe] + - [o2_feedstock, SONG_fuel_cell, oxygen, pipe] + # connect fuel cell to demand component + - [SONG_fuel_cell, electricity_load_demand, electricity, cable] +plant: + plant_life: 2 +finance_parameters: + finance_groups: + finance_model: ProFastLCO + model_inputs: + params: + analysis_start_year: 2032 + installation_time: 36 # months + inflation_rate: 0.0 # 0 for real analysis + discount_rate: 0.06 # nominal return based on 2024 ATB baseline workbook for land-based wind + debt_equity_ratio: 0.724 # 2024 ATB uses 72.4% debt for land-based wind + property_tax_and_insurance: 0.025 # percent of CAPEX estimated based on https://www.nlr.gov/docs/fy25osti/91775.pdf https://www.house.mn.gov/hrd/issinfo/clsrates.aspx + total_income_tax_rate: 0.2574 # 0.257 tax rate in 2024 atb baseline workbook, value here is based on federal (21%) and state in MN (9.8) + capital_gains_tax_rate: 0.15 # H2FAST default + sales_tax_rate: 0.0 # average combined state and local sales tax https://taxfoundation.org/location/texas/ + debt_interest_rate: 0.07 # based on 2024 ATB nominal interest rate for land-based wind + debt_type: Revolving debt # can be "Revolving debt" or "One time loan". Revolving debt is H2FAST default and leads to much lower LCOH + loan_period_if_used: 0 # H2FAST default, not used for revolving debt + cash_onhand_months: 1 # H2FAST default + admin_expense: 0.00 # percent of sales H2FAST default + capital_items: + depr_type: MACRS # can be "MACRS" or "Straight line" + depr_period: 7 # 5 years - for clean energy facilities as specified by the IRS MACRS schedule https://www.irs.gov/publications/p946#en_US_2020_publink1000107507 + refurb: [0.] + cost_adjustment_parameters: + cost_year_adjustment_inflation: 0.025 # used to adjust modeled costs to target_dollar_year + target_dollar_year: 2022 + finance_subgroups: + # h2: + # commodity: hydrogen + # commodity_stream: h2_feedstock + # technologies: [h2_feedstock] + # o2: + # commodity: oxygen + # commodity_stream: o2_feedstock + # technologies: [o2_feedstock] + electricity: + commodity: electricity + commodity_stream: SONG_fuel_cell + technologies: + - SONG_fuel_cell + - ng_feedstock + - o2_feedstock diff --git a/examples/38_song_fc/run_song_fuel_cell.py b/examples/38_song_fc/run_song_fuel_cell.py new file mode 100644 index 000000000..38ee3bb26 --- /dev/null +++ b/examples/38_song_fc/run_song_fuel_cell.py @@ -0,0 +1,18 @@ +import numpy as np + +from h2integrate import H2IntegrateModel + + +# Create a H2Integrate model +model = H2IntegrateModel("38_song_fc.yaml") + +# Setup the model +model.setup() + +# Set fuel cell demand profile +demand_profile = np.ones(8760) * 1000 +model.prob.set_val("SONG_fuel_cell.electricity_set_point", demand_profile, units="kW") + +# Run model +model.run() +model.post_process() diff --git a/examples/38_song_fc/tech_config.yaml b/examples/38_song_fc/tech_config.yaml new file mode 100644 index 000000000..fa4a67e0a --- /dev/null +++ b/examples/38_song_fc/tech_config.yaml @@ -0,0 +1,65 @@ +name: technology_config +description: This hybrid plant produces ammonia +technologies: + SONG_fuel_cell: + performance_model: + model: SONGFuelCellPerformanceModel + cost_model: + model: SONGFuelCellCostModel + model_inputs: + shared_parameters: + system_capacity_kw: 1500 + performance_parameters: + n_stacks: 60 + stack_temperature_K: 1073 # Not used yet + cost_parameters: + capex_stack_per_kw: 800 + capex_fuel_supply_per_kw: 140.67 + capex_air_supply_per_kw: 138.79 + capex_cooling_per_kw: 64.14 + capex_controls_instrumentation_per_kw: 106.52 + capex_electrical_per_kw: 447.83 + capex_assembly_per_kw: 63.51 + capex_additional_labor_per_kw: 137.96 + fixed_opex_per_kw_per_year: 31 + cost_year: 2026 + electricity_load_demand: + performance_model: + model: GenericDemandComponent + model_inputs: + performance_parameters: + commodity: electricity + commodity_rate_units: kW + demand_profile: 1000 + ng_feedstock: + performance_model: + model: FeedstockPerformanceModel + cost_model: + model: FeedstockCostModel + model_inputs: + shared_parameters: + commodity: natural_gas + commodity_rate_units: MMBtu/h + performance_parameters: + rated_capacity: 20.0 # MMBtu/h of natural gas + cost_parameters: + cost_year: 2022 + price: 0.1 + annual_cost: 0. + start_up_cost: 0.0 + o2_feedstock: + performance_model: + model: FeedstockPerformanceModel + cost_model: + model: FeedstockCostModel + model_inputs: + shared_parameters: + commodity: oxygen + commodity_rate_units: kg/h + performance_parameters: + rated_capacity: 400.0 # kg/h of oxygen + cost_parameters: + cost_year: 2022 + price: 0.5 # $0.5/kg of oxygen. price is extremely variable + annual_cost: 0. + start_up_cost: 0.0 diff --git a/h2integrate/converters/natural_gas/SO_NG_fuel_cell.py b/h2integrate/converters/natural_gas/SO_NG_fuel_cell.py index 3b13176f4..57b0c6379 100644 --- a/h2integrate/converters/natural_gas/SO_NG_fuel_cell.py +++ b/h2integrate/converters/natural_gas/SO_NG_fuel_cell.py @@ -97,10 +97,11 @@ def calc_current(power_ref, cell_area, n_cells, stack_number): # convert power_ref to Watts power_ref = power_ref * 1e3 + # clip very high power_ref to max power in curve power_density = power_ref / cell_area / stack_number / n_cells - # print("Power density", power_density) + power_density_clip = min(power_density, max(power_curve)) - I_cell = max(power_I_curve(power_density), 0) + I_cell = max(power_I_curve(power_density_clip), 0) V_cell = V_I_curve(I_cell) I_cell = I_cell * cell_area return I_cell, V_cell @@ -113,6 +114,7 @@ class SONGFuelCellPerformanceModel(PerformanceModelBaseClass): The model calculates electricity output based on natural gas and oxygen inputs, with current and voltage determined from power density using IV curves. Produces water and carbon dioxide as byproducts. + Possible source: https://www.pnnl.gov/main/publications/external/technical_reports/PNNL-18338.pdf where: - natural_gas_in is the mass flow rate of natural gas in kg/hr @@ -147,7 +149,7 @@ def setup(self): "natural_gas_in", val=0.0, shape=self.n_timesteps, - units="kg/h", + units="MMBtu/h", ) self.add_input( @@ -183,7 +185,7 @@ def setup(self): "natural_gas_consumed", val=0.0, shape=self.n_timesteps, - units="kg/h", + units="MMBtu/h", desc="Mass flow rate of natural gas consumed by the fuel cell", ) @@ -234,11 +236,23 @@ def compute(self, inputs, outputs): # calculate max input and output inputs["system_capacity"] # plant capacity in kW - inputs["natural_gas_in"] # kg/h - inputs["oxygen_in"] # kg/h - inputs["stack_temperature"] + natural_gas_in = inputs["natural_gas_in"] # MMBtu/h + oxygen_in = inputs["oxygen_in"] # kg/h + stack_temperature = inputs["stack_temperature"] # fuel_cell_efficiency = inputs["fuel_cell_efficiency"] + # Convert natural gas input from MMBtu/h to kg/h for CH4 + # Assume 1 MMBtu ~ 1 MCF + # Use ideal gas law to account for temperature difference + stack_temperature = inputs["stack_temperature"] + ng_density_cold = 0.717 # kg/m^3 at 25 deg C and 1 atm + ng_density_hot = ng_density_cold * (298.15 / stack_temperature) # kg/m^3 at stack temp + + # Use Mass = volume * density, and 1 MCF = 1000 ft^3 = 28.3168 m^3 + mcf_to_kg = 28.3168 * ng_density_hot # kg + print(natural_gas_in, mcf_to_kg) + natural_gas_kg_hr = natural_gas_in * mcf_to_kg # Convert MMBtu/h to kg/h for CH4 + # Add consumption of water for steam reforming of natural gas to hydrogen # Set calculation constants: @@ -287,8 +301,8 @@ def compute(self, inputs, outputs): for i in range(self.n_timesteps): power_reference = inputs[f"{self.commodity}_command_value"][i] - inputs["natural_gas_in"][i] - inputs["oxygen_in"][i] + natural_gas_in_loop = natural_gas_kg_hr[i] + oxygen_in_loop = oxygen_in[i] # Find current and voltage from IV curve with power setpoint I_cell, V_cell = calc_current( @@ -303,12 +317,30 @@ def compute(self, inputs, outputs): self.dt * self.config.n_stacks * self.n_cells ) # kg/time step - # print("NG and O2 consumed per hour", NG_consumed_rate, O2_consumed_rate) - # print(self.stack_size, self.n_cells) - - # TODO: if NG_consumed_rate > NGin or O2_consumed_rate > O2in: - # print("Not enough NG or O2 for this power point") - # implement an adjustment based on H2 & O2 available + if NG_consumed_rate > natural_gas_in_loop or O2_consumed_rate > oxygen_in_loop: + # implement an adjustment based on H2 & O2 available + new_i_ng = ( + natural_gas_in_loop + / (self.dt * self.config.n_stacks * self.n_cells) + * (8.0 * self.f_c) + / (self.N_series * self.M_CH4) + ) + new_i_o2 = ( + oxygen_in_loop + / (self.dt * self.config.n_stacks * self.n_cells) + * (4.0 * self.f_c) + / (self.N_series * self.M_O2) + ) + I_cell = min(new_i_ng, new_i_o2) + # TODO: recalc voltage based on new current + print("Not enough NG or O2 for this power point") + # Calculate natural gas and oxygen consumed + NG_consumed_rate = ((I_cell * self.N_series * self.M_CH4) / (8.0 * self.f_c)) * ( + self.dt * self.config.n_stacks * self.n_cells + ) # kg/time step + O2_consumed_rate = ((I_cell * self.N_series * self.M_O2) / (4.0 * self.f_c)) * ( + self.dt * self.config.n_stacks * self.n_cells + ) # kg/time step # Compute electricity from the system electricity_produced = ( @@ -348,7 +380,7 @@ def compute(self, inputs, outputs): outputs["capacity_factor"] = outputs["total_electricity_produced"] / ( self.config.system_capacity_kw * self.n_timesteps * (self.dt / 3600) ) - outputs["natural_gas_consumed"] = ng_consumed + outputs["natural_gas_consumed"] = ng_consumed / mcf_to_kg # Convert back to MMBtu/h outputs["oxygen_consumed"] = o2_consumed outputs["water_out"] = h2o_generated outputs["carbon_dioxide_out"] = co2_generated diff --git a/h2integrate/converters/natural_gas/test/test_SO_NG_fuel_cell.py b/h2integrate/converters/natural_gas/test/test_SO_NG_fuel_cell.py new file mode 100644 index 000000000..b713cc323 --- /dev/null +++ b/h2integrate/converters/natural_gas/test/test_SO_NG_fuel_cell.py @@ -0,0 +1,314 @@ +import numpy as np +import pytest +import openmdao.api as om +from pytest import fixture + +from h2integrate.converters.natural_gas.SO_NG_fuel_cell import ( + SONGFuelCellCostModel, + SONGFuelCellPerformanceModel, +) + + +@fixture +def plant_config(): + plant_config = { + "plant": { + "plant_life": 1, + "simulation": { + "n_timesteps": 48, + "dt": 3600, + }, + }, + } + return plant_config + + +@fixture +def tech_config(): + config = { + "model_inputs": { + "performance_parameters": { + "system_capacity_kw": 1500.0, + "n_stacks": 60, + "stack_temperature_K": 1073.0, + } + } + } + return config + + +@fixture +def cost_config(): + config = { + "model_inputs": { + "cost_parameters": { + "system_capacity_kw": 1500.0, + "capex_stack_per_kw": 800.0, + "capex_fuel_supply_per_kw": 140.67, + "capex_air_supply_per_kw": 138.79, + "capex_cooling_per_kw": 64.14, + "capex_controls_instrumentation_per_kw": 106.52, + "capex_electrical_per_kw": 447.83, + "capex_assembly_per_kw": 63.51, + "capex_additional_labor_per_kw": 137.96, + "fixed_opex_per_kw_per_year": 31.0, + "cost_year": 2026, + } + } + } + return config + + +@pytest.mark.regression +def test_fuel_cell_performance(tech_config, plant_config, subtests): + n_timesteps = int(plant_config["plant"]["simulation"]["n_timesteps"]) + + prob = om.Problem() + + fuel_cell = SONGFuelCellPerformanceModel( + plant_config=plant_config, tech_config=tech_config, driver_config={} + ) + + prob.model.add_subsystem("fuel_cell", fuel_cell, promotes=["*"]) + + prob.setup() + + # Provide ample natural gas and oxygen to run at the default command (rated capacity) + natural_gas_input = np.ones(n_timesteps) * 20.0 # MMBtu/h + oxygen_input = np.ones(n_timesteps) * 2000.0 # kg/h + + prob.set_val("fuel_cell.natural_gas_in", natural_gas_input, units="MMBtu/h") + prob.set_val("fuel_cell.oxygen_in", oxygen_input, units="kg/h") + prob.set_val("fuel_cell.electricity_command_value", np.ones(n_timesteps) * 1000.0, units="kW") + + prob.run_model() + + with subtests.test("max electricity output"): + assert ( + pytest.approx(np.max(prob.get_val("fuel_cell.electricity_out", units="kW")), rel=1e-2) + == 1000.0 + ) + + with subtests.test("electricity out"): + assert ( + pytest.approx(np.sum(prob.get_val("fuel_cell.electricity_out", units="kW")), rel=1e-6) + == 48209.2 + ) + + with subtests.test("capacity_factor"): + assert ( + pytest.approx(prob.get_val("fuel_cell.capacity_factor", units="unitless"), rel=1e-2) + == 0.669 + ) + + with subtests.test("annual_electricity_production"): + assert ( + pytest.approx( + prob.get_val("fuel_cell.annual_electricity_produced", units="kW*h/year"), rel=1e-2 + ) + == 8760000.86 + ) + + with subtests.test("rated_electricity_production"): + assert ( + pytest.approx( + prob.get_val("fuel_cell.rated_electricity_production", units="kW"), rel=1e-6 + ) + == 1500.0 + ) + + with subtests.test("total_electricity_produced"): + assert ( + pytest.approx( + prob.get_val("fuel_cell.total_electricity_produced", units="kW*h"), rel=1e-6 + ) + == 48209.2 + ) + + with subtests.test("natural gas consumed"): + assert ( + pytest.approx( + np.sum(prob.get_val("fuel_cell.natural_gas_consumed", units="MMBtu/h")), rel=1e-6 + ) + == 807.926382702 + ) + + with subtests.test("oxygen consumed"): + assert ( + pytest.approx(np.sum(prob.get_val("fuel_cell.oxygen_consumed", units="kg/h")), rel=1e-6) + == 18186.36156 + ) + + with subtests.test("water out"): + assert ( + pytest.approx(np.sum(prob.get_val("fuel_cell.water_out", units="kg/h")), rel=1e-6) + == 20459.65675 + ) + + with subtests.test("carbon dioxide out"): + assert ( + pytest.approx( + np.sum(prob.get_val("fuel_cell.carbon_dioxide_out", units="kg/h")), rel=1e-6 + ) + == 12503.12357 + ) + + +@pytest.mark.unit +def test_fuel_cell_demand(tech_config, plant_config, subtests): + n_timesteps = int(plant_config["plant"]["simulation"]["n_timesteps"]) + + prob = om.Problem() + + fuel_cell = SONGFuelCellPerformanceModel( + plant_config=plant_config, tech_config=tech_config, driver_config={} + ) + + prob.model.add_subsystem("fuel_cell", fuel_cell, promotes=["*"]) + + prob.setup() + + # Provide ample feedstock supply for most timesteps; constrain a few to test + # feedstock-limited dynamics. + natural_gas_input = np.ones(n_timesteps) * 30.0 # MMBtu/h + oxygen_input = np.ones(n_timesteps) * 2000.0 # kg/h + + # Edge cases for feedstock supply at timesteps 5-7 + natural_gas_input[5:8] = ( + 0.0, # zero NG supply -> output should collapse to zero + 1, # severely limited NG supply -> output reduced + 20.0, # ample NG supply, but O2 will be limited below + ) + oxygen_input[5:8] = ( + 2000.0, + 2000.0, + 0.0, # zero O2 supply -> output should collapse to zero + ) + + prob.set_val("fuel_cell.natural_gas_in", natural_gas_input, units="MMBtu/h") + prob.set_val("fuel_cell.oxygen_in", oxygen_input, units="kg/h") + + elec_set_point = np.ones(n_timesteps) * 1500.0 # kW + + # First 5 timesteps test set-point edge cases (with ample feedstock). + # Timesteps 5-7 test feedstock-limited dynamics at rated set point. + elec_set_point[:5] = ( + 1500.0, # set point equal to system capacity + 500.0, # set point below system capacity + 2500.0, # very high set point (should be clipped to system capacity) + 0.0, # zero set point + 750.0, # set point at half of system capacity + ) + + prob.set_val("fuel_cell.electricity_command_value", elec_set_point, units="kW") + + prob.run_model() + + electricity_output = prob.get_val("fuel_cell.electricity_out", units="kW") + ng_consumed = prob.get_val("fuel_cell.natural_gas_consumed", units="MMBtu/h") + o2_consumed = prob.get_val("fuel_cell.oxygen_consumed", units="kg/h") + water_out = prob.get_val("fuel_cell.water_out", units="kg/h") + co2_out = prob.get_val("fuel_cell.carbon_dioxide_out", units="kg/h") + + with subtests.test("output bounded by system capacity"): + assert np.max(electricity_output) <= 1500.0 + 1e-6 + + with subtests.test("output non-negative"): + assert np.min(electricity_output) >= 0.0 + + with subtests.test("output clipped to system capacity at rated set point"): + assert electricity_output[0] == pytest.approx(1500.0, rel=1e-2) + + with subtests.test("output follows reduced set point"): + assert electricity_output[1] == pytest.approx(500.0, rel=1e-2) + + with subtests.test("very high set point clipped to system capacity"): + assert electricity_output[2] == pytest.approx(1500.0, rel=1e-2) + + with subtests.test("zero set point yields zero output"): + assert electricity_output[3] == pytest.approx(0.0, abs=1e-6) + + with subtests.test("half-rated set point yields ~half output"): + assert electricity_output[4] == pytest.approx(750.0, rel=5e-2) + + with subtests.test("natural gas consumed non-negative"): + assert np.min(ng_consumed) >= 0.0 + + with subtests.test("oxygen consumed non-negative"): + assert np.min(o2_consumed) >= 0.0 + + with subtests.test("water produced non-negative"): + assert np.min(water_out) >= 0.0 + + with subtests.test("CO2 produced non-negative"): + assert np.min(co2_out) >= 0.0 + + with subtests.test("zero set point yields zero consumption and byproducts"): + assert ng_consumed[3] == pytest.approx(0.0, abs=1e-6) + assert o2_consumed[3] == pytest.approx(0.0, abs=1e-6) + assert water_out[3] == pytest.approx(0.0, abs=1e-6) + assert co2_out[3] == pytest.approx(0.0, abs=1e-6) + + # Feedstock-limited dynamics: + + with subtests.test("zero NG supply collapses output to zero"): + assert electricity_output[5] == pytest.approx(0.0, abs=1e-6) + assert ng_consumed[5] == pytest.approx(0.0, abs=1e-6) + assert o2_consumed[5] == pytest.approx(0.0, abs=1e-6) + assert water_out[5] == pytest.approx(0.0, abs=1e-6) + assert co2_out[5] == pytest.approx(0.0, abs=1e-6) + + with subtests.test("limited NG supply reduces electricity output and byproducts"): + assert electricity_output[6] == pytest.approx(55.547789, rel=1e-2) + assert ng_consumed[6] == pytest.approx(1.0, rel=1e-2) + assert o2_consumed[6] == pytest.approx(22.509924, rel=1e-2) + assert water_out[6] == pytest.approx(25.323664, rel=1e-2) + assert co2_out[6] == pytest.approx(15.475572, rel=1e-2) + + with subtests.test("zero O2 supply collapses output to zero"): + assert electricity_output[7] == pytest.approx(0.0, abs=1e-6) + assert ng_consumed[7] == pytest.approx(0.0, abs=1e-6) + assert o2_consumed[7] == pytest.approx(0.0, abs=1e-6) + assert water_out[7] == pytest.approx(0.0, abs=1e-6) + assert co2_out[7] == pytest.approx(0.0, abs=1e-6) + + +@pytest.mark.regression +def test_fuel_cell_cost(cost_config, plant_config, subtests): + prob = om.Problem() + + fuel_cell_cost = SONGFuelCellCostModel( + plant_config=plant_config, tech_config=cost_config, driver_config={} + ) + + prob.model.add_subsystem("fuel_cell_cost", fuel_cell_cost, promotes=["*"]) + + prob.setup() + + prob.run_model() + + cp = cost_config["model_inputs"]["cost_parameters"] + expected_unit_capex = ( + cp["capex_stack_per_kw"] + + cp["capex_fuel_supply_per_kw"] + + cp["capex_air_supply_per_kw"] + + cp["capex_cooling_per_kw"] + + cp["capex_controls_instrumentation_per_kw"] + + cp["capex_electrical_per_kw"] + + cp["capex_assembly_per_kw"] + + cp["capex_additional_labor_per_kw"] + ) + expected_capex = cp["system_capacity_kw"] * expected_unit_capex + expected_opex = cp["system_capacity_kw"] * cp["fixed_opex_per_kw_per_year"] + + with subtests.test("capex value"): + assert ( + pytest.approx(prob.get_val("fuel_cell_cost.CapEx", units="USD"), rel=1e-6) + == expected_capex + ) + + with subtests.test("opex value"): + assert ( + pytest.approx(prob.get_val("fuel_cell_cost.OpEx", units="USD/year"), rel=1e-6) + == expected_opex + ) From 7e6d3625cd761b2d166b9e655647ea86ac8a33c1 Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Mon, 6 Jul 2026 16:50:37 -0400 Subject: [PATCH 18/21] Update PEM regression tests --- .../hydrogen/test/test_PEM_fuel_cell.py | 88 ++++++++++++------- .../converters/natural_gas/SO_NG_fuel_cell.py | 31 +++---- 2 files changed, 70 insertions(+), 49 deletions(-) diff --git a/h2integrate/converters/hydrogen/test/test_PEM_fuel_cell.py b/h2integrate/converters/hydrogen/test/test_PEM_fuel_cell.py index 4f2ffc68f..5c4735682 100644 --- a/h2integrate/converters/hydrogen/test/test_PEM_fuel_cell.py +++ b/h2integrate/converters/hydrogen/test/test_PEM_fuel_cell.py @@ -15,7 +15,7 @@ def plant_config(): "plant": { "plant_life": 30, "simulation": { - "n_timesteps": 8760, + "n_timesteps": 48, "dt": 3600, }, }, @@ -79,13 +79,14 @@ def test_fuel_cell_performance(tech_config, plant_config, subtests): prob.set_val("fuel_cell.hydrogen_in", hydrogen_input, units="kg/h") prob.set_val("fuel_cell.oxygen_in", oxygen_input, units="kg/h") + prob.set_val("fuel_cell.electricity_command_value", np.ones(n_timesteps) * 1000.0, units="kW") prob.run_model() electricity_output = prob.get_val("fuel_cell.electricity_out", units="kW") - hydrogen_consumed = prob.get_val("fuel_cell.hydrogen_consumed", units="kg/h") - oxygen_consumed = prob.get_val("fuel_cell.oxygen_consumed", units="kg/h") - water_out = prob.get_val("fuel_cell.water_out", units="kg/h") + prob.get_val("fuel_cell.hydrogen_consumed", units="kg/h") + prob.get_val("fuel_cell.oxygen_consumed", units="kg/h") + prob.get_val("fuel_cell.water_out", units="kg/h") with subtests.test("max electricity output bounded by system capacity"): assert np.max(electricity_output) <= 1500.0 + 1e-6 @@ -106,33 +107,60 @@ def test_fuel_cell_performance(tech_config, plant_config, subtests): prob.get_val("fuel_cell.total_electricity_produced", units="kW*h"), rel=1e-6 ) == np.sum(electricity_output) - with subtests.test("annual_electricity_produced equals total for full-year sim"): - assert pytest.approx( - prob.get_val("fuel_cell.annual_electricity_produced", units="kW*h/year"), rel=1e-6 - ) == np.sum(electricity_output) + with subtests.test("electricity out"): + assert ( + pytest.approx(np.sum(prob.get_val("fuel_cell.electricity_out", units="kW")), rel=1e-6) + == 48209.20157 + ) - with subtests.test("capacity_factor matches definition"): - assert pytest.approx( - prob.get_val("fuel_cell.capacity_factor", units="unitless"), rel=1e-6 - ) == np.sum(electricity_output) / (1500.0 * n_timesteps) - - with subtests.test("hydrogen consumed is non-negative and bounded by supply"): - assert np.min(hydrogen_consumed) >= 0.0 - assert np.max(hydrogen_consumed) <= 200.0 + 1e-6 - - with subtests.test("oxygen consumed is non-negative and bounded by supply"): - assert np.min(oxygen_consumed) >= 0.0 - assert np.max(oxygen_consumed) <= 2000.0 + 1e-6 - - with subtests.test("water produced is non-negative"): - assert np.min(water_out) >= 0.0 - - with subtests.test("stoichiometric mass balance H2 + O2 -> H2O"): - # Combined reactant mass should equal water produced (mass conservation) - np.testing.assert_allclose( - hydrogen_consumed + oxygen_consumed, - water_out, - rtol=1e-3, + with subtests.test("capacity_factor"): + assert ( + pytest.approx(prob.get_val("fuel_cell.capacity_factor", units="unitless"), rel=1e-2) + == 0.669 + ) + + with subtests.test("annual_electricity_production"): + assert ( + pytest.approx( + prob.get_val("fuel_cell.annual_electricity_produced", units="kW*h/year"), rel=1e-2 + ) + == 8760000.86 + ) + + with subtests.test("rated_electricity_production"): + assert ( + pytest.approx( + prob.get_val("fuel_cell.rated_electricity_production", units="kW"), rel=1e-6 + ) + == 1500.0 + ) + + with subtests.test("total_electricity_produced"): + assert ( + pytest.approx( + prob.get_val("fuel_cell.total_electricity_produced", units="kW*h"), rel=1e-6 + ) + == 48209.20157 + ) + + with subtests.test("hydrogen consumed"): + assert ( + pytest.approx( + np.sum(prob.get_val("fuel_cell.hydrogen_consumed", units="kg/h")), rel=1e-6 + ) + == 2291.481556 + ) + + with subtests.test("oxygen consumed"): + assert ( + pytest.approx(np.sum(prob.get_val("fuel_cell.oxygen_consumed", units="kg/h")), rel=1e-6) + == 18186.361561 + ) + + with subtests.test("water out"): + assert ( + pytest.approx(np.sum(prob.get_val("fuel_cell.water_out", units="kg/h")), rel=1e-6) + == 20459.65675 ) diff --git a/h2integrate/converters/natural_gas/SO_NG_fuel_cell.py b/h2integrate/converters/natural_gas/SO_NG_fuel_cell.py index 57b0c6379..b7b347d70 100644 --- a/h2integrate/converters/natural_gas/SO_NG_fuel_cell.py +++ b/h2integrate/converters/natural_gas/SO_NG_fuel_cell.py @@ -241,6 +241,10 @@ def compute(self, inputs, outputs): stack_temperature = inputs["stack_temperature"] # fuel_cell_efficiency = inputs["fuel_cell_efficiency"] + ############################################################################ + # Can convert from MMBtu to Joules in openmdao + # Then can convert from Joules to kg using some heating value or enthalpy + # Convert natural gas input from MMBtu/h to kg/h for CH4 # Assume 1 MMBtu ~ 1 MCF # Use ideal gas law to account for temperature difference @@ -248,7 +252,7 @@ def compute(self, inputs, outputs): ng_density_cold = 0.717 # kg/m^3 at 25 deg C and 1 atm ng_density_hot = ng_density_cold * (298.15 / stack_temperature) # kg/m^3 at stack temp - # Use Mass = volume * density, and 1 MCF = 1000 ft^3 = 28.3168 m^3 + # Use Mass = volume * density, and 1 MCF = 1e3 ft^3 = 28.3168 m^3 mcf_to_kg = 28.3168 * ng_density_hot # kg print(natural_gas_in, mcf_to_kg) natural_gas_kg_hr = natural_gas_in * mcf_to_kg # Convert MMBtu/h to kg/h for CH4 @@ -393,21 +397,15 @@ def compute(self, inputs, outputs): class SONGFuelCellCostConfig(CostModelBaseConfig): """Configuration class for the solid oxide natural gas fuel cell cost model. - Fields include `system_capacity_kw`, `capex_stack_per_kw`, `capex_fuel_supply_per_kw`, - `capex_air_supply_per_kw`, `capex_cooling_per_kw`, `capex_controls_instrumentation_per_kw`, - `capex_electrical_per_kw`, `capex_assembly_per_kw`, `capex_additional_labor_per_kw`, - and `fixed_opex_per_kw_per_year`. The `cost_year` field is inherited from `CostModelBaseConfig`. + Fields include `system_capacity_kw`, `capex_stack_per_kw`, + `capex_bop`, `capex_indirect_costs_per_kw`, and `fixed_opex_per_kw_per_year`. + The `cost_year` field is inherited from `CostModelBaseConfig`. """ system_capacity_kw: float = field(validator=gte_zero) capex_stack_per_kw: float = field(validator=gte_zero) - capex_fuel_supply_per_kw: float = field(validator=gte_zero) - capex_air_supply_per_kw: float = field(validator=gte_zero) - capex_cooling_per_kw: float = field(validator=gte_zero) - capex_controls_instrumentation_per_kw: float = field(validator=gte_zero) - capex_electrical_per_kw: float = field(validator=gte_zero) - capex_assembly_per_kw: float = field(validator=gte_zero) - capex_additional_labor_per_kw: float = field(validator=gte_zero) + capex_bop: float = field(validator=gte_zero) + capex_indirect_costs_per_kw: float = field(validator=gte_zero) fixed_opex_per_kw_per_year: float = field(validator=gte_zero) @@ -442,13 +440,8 @@ def setup(self): self.add_input( "unit_capex", val=self.config.capex_stack_per_kw - + self.config.capex_fuel_supply_per_kw - + self.config.capex_air_supply_per_kw - + self.config.capex_cooling_per_kw - + self.config.capex_controls_instrumentation_per_kw - + self.config.capex_electrical_per_kw - + self.config.capex_assembly_per_kw - + self.config.capex_additional_labor_per_kw, + + self.config.capex_bop + + self.config.capex_indirect_costs_per_kw, units="USD/kW", desc="Capital cost per unit capacity", ) From a1f97ea0d6efd3ca1aea27de9ef89791a8fcd458 Mon Sep 17 00:00:00 2001 From: genevievestarke <103534902+genevievestarke@users.noreply.github.com> Date: Tue, 7 Jul 2026 13:19:25 -0400 Subject: [PATCH 19/21] Update h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py Co-authored-by: elenya-grant <116225007+elenya-grant@users.noreply.github.com> --- h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py b/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py index 25a01eb8a..2ea47017d 100644 --- a/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py +++ b/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py @@ -259,7 +259,7 @@ def compute(self, inputs, outputs): self.max_cell_power_density = 0.000334 # is n_cells = N_series? self.N_series = 1 - self.stack_size = self.config.system_capacity_kw / self.config.n_stacks + self.stack_size = inputs["system_capacity"] / self.config.n_stacks self.cell_active_area = 400 # [cm^2] from Battelle (https://www.energy.gov/sites/prod/files/2018/02/f49/fcto_battelle_mfg_cost_analysis_1%20_to_25kw_pp_chp_fc_systems_jan2017_0.pdf) self.n_cells = round( self.stack_size / (self.cell_active_area * self.max_cell_power_density) From af7944aa1a8c621264bb90a10b9a672a2cc42c19 Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Fri, 10 Jul 2026 19:08:13 -0400 Subject: [PATCH 20/21] Add PEM FC updates from other branch --- .../converters/hydrogen/PEM_h2_fuel_cell.py | 274 ++++++------------ .../hydrogen/test/test_PEM_fuel_cell.py | 26 +- 2 files changed, 102 insertions(+), 198 deletions(-) diff --git a/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py b/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py index 2ea47017d..f31887497 100644 --- a/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py +++ b/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py @@ -3,6 +3,7 @@ from h2integrate.core.utilities import BaseConfig, merge_shared_inputs from h2integrate.core.validators import gte_zero +from h2integrate.tools.constants import H_MW, O2_MW, faraday from h2integrate.core.model_baseclasses import ( CostModelBaseClass, CostModelBaseConfig, @@ -30,80 +31,46 @@ class PEMH2FuelCellPerformanceConfig(BaseConfig): # fuel_cell_efficiency_hhv: float = field(validator=range_val(0, 1)) -def calc_current(power_ref, cell_area, n_cells, stack_number): +def calc_current(system_power_reference, cell_area, n_cells, n_stacks): + """_summary_ + + Args: + system_power_reference (np.ndarray): power demanded of the entire system in W + cell_area (float): cell active area in cm^2 + n_cells (int): number of cells per stack + n_stacks (int): number of stacks in the system + + Returns: + tuple(np.ndarray, np.poly1d): stack current and + function to convert from current density to voltage + """ # Calculates the current and voltage from IV curve based on power reference - current_curve = [ - 0.0356, - 0.05413333, - 0.0796, - 0.11366667, - 0.244, - 0.454, - # 0.70366667, - # 0.96933333, - # 1.24, - # 1.52666667, - # 1.80333333, - # 2.07, - # 2.32, - # 2.54333333, - # 2.73666667, - # 2.9, - ] # in A - voltage_curve = [ - 0.987, - 0.936, - 0.884, - 0.838, - 0.786, - 0.736, - # 0.686, - # 0.636, - # 0.586, - # 0.53566667, - # 0.486, - # 0.436, - # 0.386, - # 0.33533333, - # 0.286, - # 0.236, - ] # in V - power_curve = [ - 35.16666667, - 50.53333333, - 70.33333333, - 95.46666667, - 191.66666667, - 334.33333333, - # 482.66666667, - # 616.66666667, - # 729.0, - # 817.0, - # 875.33333333, - # 902.0, - # 895.0, - # 854.33333333, - # 782.66666667, - # 684.33333333, - ] - - # Change power from mW to W - power_curve = [x / 1e3 for x in power_curve] - - power_coefs = np.polyfit(power_curve, current_curve, 5) - power_I_curve = np.poly1d(power_coefs) - V_coefs = np.polyfit(current_curve, voltage_curve, 5) - V_I_curve = np.poly1d(V_coefs) - - # convert power_ref to Watts - power_ref = power_ref * 1e3 - power_density = power_ref / cell_area / stack_number / n_cells - # print("Power density", power_density) - - I_cell = max(power_I_curve(power_density), 0) - V_cell = V_I_curve(I_cell) - I_cell = I_cell * cell_area - return I_cell, V_cell + J_curve = np.array([0.0356, 0.05413333, 0.0796, 0.11366667, 0.244, 0.454]) # in A/cm^2 + voltage_curve = np.array([0.987, 0.936, 0.884, 0.838, 0.786, 0.736]) # in V + power_curve = ( + np.array([35.16666667, 50.53333333, 70.33333333, 95.46666667, 191.66666667, 334.33333333]) + / 1e3 + ) # in W/cm^2 + + # function to calculate voltage from current density + V_coefs = np.polyfit(J_curve, voltage_curve, 5) + V_J_curve = np.poly1d(V_coefs) + + # function to calculate current density from power + # I_curve = J_curve * cell_area + # P_curve = I_curve * (n_cells * voltage_curve) + stack_P_curve = power_curve * cell_area * n_cells + J_coefs = np.polyfit(stack_P_curve, J_curve, 5) + J_P_curve = np.poly1d(J_coefs) + + # Calculate power per stack and power density + power_per_stack = system_power_reference / n_stacks + + stack_current_density = J_P_curve(power_per_stack) # convert to kW for the curve + stack_current = stack_current_density * cell_area * n_cells # in A + stack_current = np.clip(stack_current, a_min=0.0, a_max=None) # clip negative values + + return stack_current, V_J_curve class PEMH2FuelCellPerformanceModel(PerformanceModelBaseClass): @@ -233,131 +200,76 @@ def compute(self, inputs, outputs): oxygen_consumed, water_out, and various electricity production quantities. """ - # calculate max input and output - inputs["system_capacity"] # plant capacity in kW - inputs["hydrogen_in"] # kg/h - inputs["oxygen_in"] # kg/h - inputs["stack_temperature"] - # fuel_cell_efficiency = inputs["fuel_cell_efficiency"] - # Set calculation constants: - self.f_c = 96485.33 # Faraday's constant in A/mol - self.M_H2 = 0.002016 # Molar mass of H2 in kg/mol - self.M_O2 = 0.032 # Molar mass of O2 in kg/mol - self.M_H2O = 0.018 # Molar mass of H2O in kg/mol - self.Tref = 298.15 # Standard room temperature in K [25 deg Celsius] - self.cp_H2 = 14300 # Specific heat of H2 in J/(kg*K) - self.cp_air = 1005 # Specific heat of air in J/(kg*K) - self.cp_H2O = 4184 # Specific heat of water in J/(kg*K) - self.cp_N2 = 1040 # Specific heat of nitrogen in J/(kg*K) - self.cp_O2 = 918 # Specific heat of oxygen in J/(kg*K) - self.hhv_h2 = 141.8 * 1e6 # Higher heating value of hydrogen in J/kg - self.hhv_air = 0 # No higher heating value of air - self.hhv_H2O = 2260 # Higher heating value of water in J/kg + M_H2 = H_MW * 2 / 1000 # Molar mass of H2 in kg/mol + M_O2 = O2_MW / 1000 # Molar mass of O2 in kg/mol + M_H2O = M_H2 + M_O2 / 2 # Molar mass of H2O in kg/mol # Sizing the cells - self.max_cell_power_density = 0.000334 - # is n_cells = N_series? - self.N_series = 1 - self.stack_size = inputs["system_capacity"] / self.config.n_stacks - self.cell_active_area = 400 # [cm^2] from Battelle (https://www.energy.gov/sites/prod/files/2018/02/f49/fcto_battelle_mfg_cost_analysis_1%20_to_25kw_pp_chp_fc_systems_jan2017_0.pdf) - self.n_cells = round( - self.stack_size / (self.cell_active_area * self.max_cell_power_density) + max_cell_power_density = 0.000334 # in W/cm^2 + stack_size = inputs["system_capacity"][0] / self.config.n_stacks + cell_active_area = 400 # [cm^2] from Battelle (https://www.energy.gov/sites/prod/files/2018/02/f49/fcto_battelle_mfg_cost_analysis_1%20_to_25kw_pp_chp_fc_systems_jan2017_0.pdf) + n_cells = round(stack_size / (cell_active_area * max_cell_power_density)) + # Recalculate the rated power production based on final fuel cell sizing + rated_power_production = ( + max_cell_power_density * n_cells * cell_active_area * self.config.n_stacks ) - # PSUEDO CODE: - """ - 1. Receive power setpoint into fuel cell - 2. Find current with I-V curve - 3. Calculate H2 consumed and O2 consumed - 4. Check if provided H2 and O2 can meet the demand - 5. If not, adjust current - 6. Calculate power out with current and voltage - 7. Calculate the water produced from the reaction + ################## Model Calculations ################## + # 1. Receive power setpoint into fuel cell + power_reference = np.clip( + inputs[f"{self.commodity}_command_value"], a_min=0.0, a_max=rated_power_production + ) - """ + # 2. Find stack current from power reference + commanded_I_stack, V_J_curve = calc_current( + power_reference * 1e3, cell_active_area, n_cells, self.config.n_stacks + ) - h2_consumed = np.zeros(self.n_timesteps) - o2_consumed = np.zeros(self.n_timesteps) - h2o_generated = np.zeros(self.n_timesteps) - commodity_out = np.zeros(self.n_timesteps) - - for i in range(self.n_timesteps): - power_reference = inputs[f"{self.commodity}_command_value"][i] - H2in = inputs["hydrogen_in"][i] - O2in = inputs["oxygen_in"][i] - - # Find current and voltage from IV curve with power setpoint - I_cell, V_cell = calc_current( - power_reference, self.cell_active_area, self.n_cells, self.config.n_stacks - ) - - # Calculate hydrogen and oxygen consumed - H2_consumed_rate = ((I_cell * self.N_series * self.M_H2) / (2.0 * self.f_c)) * ( - self.dt * self.config.n_stacks * self.n_cells - ) # kg/time step - O2_consumed_rate = ((I_cell * self.N_series * self.M_O2) / (4.0 * self.f_c)) * ( - self.dt * self.config.n_stacks * self.n_cells - ) # kg/time step - - # print("H2 and O2 consumed per hour", H2_consumed_rate, O2_consumed_rate) - # print(self.stack_size, self.n_cells) - - # TODO: - if H2_consumed_rate > H2in or O2_consumed_rate > O2in: - # implement an adjustment based on H2 & O2 available - new_i_h2 = ( - H2in - / (self.dt * self.config.n_stacks * self.n_cells) - * (2.0 * self.f_c) - / (self.N_series * self.M_H2) - ) - new_i_o2 = ( - O2in - / (self.dt * self.config.n_stacks * self.n_cells) - * (4.0 * self.f_c) - / (self.N_series * self.M_O2) - ) - I_cell = min(new_i_h2, new_i_o2) - # TODO: recalc voltage based on new current - print("Not enough H2 or O2 for this power point") - # Calculate hydrogen and oxygen consumed - H2_consumed_rate = ((I_cell * self.N_series * self.M_H2) / (2.0 * self.f_c)) * ( - self.dt * self.config.n_stacks * self.n_cells - ) # kg/time step - O2_consumed_rate = ((I_cell * self.N_series * self.M_O2) / (4.0 * self.f_c)) * ( - self.dt * self.config.n_stacks * self.n_cells - ) # kg/time step - - # Compute electricity from the system - electricity_produced = ( - V_cell * I_cell * self.n_cells * self.config.n_stacks / 1e3 - ) # Calculated in watts, convert to kW - - # Compute H2O out - H2O_generated = ( - (I_cell * self.N_series / (2 * self.f_c)) - * self.M_H2O - * (self.dt * self.config.n_stacks * self.n_cells) - ) # in kg/time step - - h2_consumed[i] = H2_consumed_rate - o2_consumed[i] = O2_consumed_rate - h2o_generated[i] = H2O_generated - commodity_out[i] = electricity_produced + # 3. Find available hydrogen and oxygen for each timestep + h2_in_kg_per_s = inputs["hydrogen_in"] / 3600 # convert from kg/h to kg/s + o2_in_kg_per_s = inputs["oxygen_in"] / 3600 # convert from kg/h to kg/s + + # convert from kg/s to A per stack - current that feedstocks in can support + I_stack_from_h2 = (h2_in_kg_per_s * 2 * faraday) / (M_H2 * self.config.n_stacks) + I_stack_from_o2 = (o2_in_kg_per_s * 4 * faraday) / (M_O2 * self.config.n_stacks) + + # 4. Take minimum current from power reference, hydrogen available, and oxygen available + I_stack = np.minimum(commanded_I_stack, np.minimum(I_stack_from_h2, I_stack_from_o2)) + + # 5. Calculate current density and voltage from I-V curve + J_cell = I_stack / (cell_active_area * n_cells) # in A/cm^2 + I_cell = J_cell * cell_active_area # in A + V_cell = V_J_curve(J_cell) # in V + + # 6. Calculate power output from current and voltage + power_out = V_cell * I_cell * n_cells * self.config.n_stacks / 1e3 # in kW + + # 7. Calculate hydrogen and oxygen consumed and water produced + # based on electrochemical reactions + h2_consumed = ((I_cell * M_H2) / (2.0 * faraday)) * ( + self.dt * self.config.n_stacks * n_cells + ) # kg/time step + o2_consumed = ((I_cell * M_O2) / (4.0 * faraday)) * ( + self.dt * self.config.n_stacks * n_cells + ) # kg/time step + h2o_generated = ((I_cell * M_H2O) / (2 * faraday)) * ( + self.dt * self.config.n_stacks * n_cells + ) # kg/time step # Set Outputs # clip the electricity output to the system capacity - outputs["electricity_out"] = np.minimum(commodity_out, self.config.system_capacity_kw) + outputs["rated_electricity_production"] = rated_power_production + + outputs["electricity_out"] = np.minimum(power_out, rated_power_production) outputs["total_electricity_produced"] = np.sum(outputs["electricity_out"]) * ( self.dt / 3600 ) - outputs["rated_electricity_production"] = self.config.system_capacity_kw outputs["annual_electricity_produced"] = outputs["total_electricity_produced"] * ( 1 / self.fraction_of_year_simulated ) outputs["capacity_factor"] = outputs["total_electricity_produced"] / ( - self.config.system_capacity_kw * self.n_timesteps * (self.dt / 3600) + rated_power_production * self.n_timesteps * (self.dt / 3600) ) outputs["hydrogen_consumed"] = h2_consumed outputs["oxygen_consumed"] = o2_consumed diff --git a/h2integrate/converters/hydrogen/test/test_PEM_fuel_cell.py b/h2integrate/converters/hydrogen/test/test_PEM_fuel_cell.py index 5c4735682..e67288920 100644 --- a/h2integrate/converters/hydrogen/test/test_PEM_fuel_cell.py +++ b/h2integrate/converters/hydrogen/test/test_PEM_fuel_cell.py @@ -94,14 +94,6 @@ def test_fuel_cell_performance(tech_config, plant_config, subtests): with subtests.test("electricity output is non-negative"): assert np.min(electricity_output) >= 0.0 - with subtests.test("rated_electricity_production"): - assert ( - pytest.approx( - prob.get_val("fuel_cell.rated_electricity_production", units="kW"), rel=1e-6 - ) - == 1500.0 - ) - with subtests.test("total_electricity_produced matches sum of output"): assert pytest.approx( prob.get_val("fuel_cell.total_electricity_produced", units="kW*h"), rel=1e-6 @@ -130,9 +122,9 @@ def test_fuel_cell_performance(tech_config, plant_config, subtests): with subtests.test("rated_electricity_production"): assert ( pytest.approx( - prob.get_val("fuel_cell.rated_electricity_production", units="kW"), rel=1e-6 + prob.get_val("fuel_cell.rated_electricity_production", units="kW"), rel=1e-4 ) - == 1500.0 + == 1498.9 ) with subtests.test("total_electricity_produced"): @@ -154,13 +146,13 @@ def test_fuel_cell_performance(tech_config, plant_config, subtests): with subtests.test("oxygen consumed"): assert ( pytest.approx(np.sum(prob.get_val("fuel_cell.oxygen_consumed", units="kg/h")), rel=1e-6) - == 18186.361561 + == 18185.22451853 ) with subtests.test("water out"): assert ( pytest.approx(np.sum(prob.get_val("fuel_cell.water_out", units="kg/h")), rel=1e-6) - == 20459.65675 + == 20476.70602546 ) @@ -223,7 +215,7 @@ def test_fuel_cell_demand(tech_config, plant_config, subtests): water_out = prob.get_val("fuel_cell.water_out", units="kg/h") with subtests.test("output clipped to system capacity"): - assert electricity_output[0] == pytest.approx(1500.0, rel=1e-4) + assert electricity_output[0] == pytest.approx(1500.0, rel=1e-3) with subtests.test("output follows reduced set point"): # When set point is below capacity and feedstock is ample, output tracks set point @@ -237,7 +229,7 @@ def test_fuel_cell_demand(tech_config, plant_config, subtests): assert electricity_output[3] == pytest.approx(0.0, abs=1e-6) with subtests.test("limited hydrogen feedstock supply yields reduced output"): - assert electricity_output[4] < 20.0 + assert electricity_output[4] < 30.0 assert electricity_output[4] > 0.0 with subtests.test("zero set point yields zero output"): @@ -245,15 +237,15 @@ def test_fuel_cell_demand(tech_config, plant_config, subtests): # Test hydrogen_consumed, oxygen_consumed, and water_out for the first 6 timesteps with subtests.test("hydrogen consumed"): - expected_h2_consumed = [76.589328, 22.920501, 76.589328, 0, 1, 0.0] + expected_h2_consumed = [76.501247, 22.920501, 76.501247, 0, 1, 0.0] np.testing.assert_allclose(hydrogen_consumed[:6], expected_h2_consumed, rtol=1e-4) with subtests.test("oxygen consumed"): - expected_o2_consumed = [607.851812, 181.908739, 607.851812, 0, 7.936508, 0.0] + expected_o2_consumed = [607.114803, 181.897366, 607.114803, 0, 7.936012, 0.0] np.testing.assert_allclose(oxygen_consumed[:6], expected_o2_consumed, rtol=1e-4) with subtests.test("water out"): - expected_water_out = [683.833289, 204.647332, 683.833289, 0.0, 8.928571, 0.0] + expected_water_out = [683.61605, 204.817867, 683.61605, 0.0, 8.936012, 0.0] np.testing.assert_allclose(water_out[:6], expected_water_out, rtol=1e-4) From 59001b20377ada59ce97ec9db2b567a7dc3ba3e1 Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Fri, 10 Jul 2026 19:27:20 -0400 Subject: [PATCH 21/21] Update SO NG FC --- .../converters/natural_gas/SO_NG_fuel_cell.py | 301 +++++++----------- .../natural_gas/test/test_SO_NG_fuel_cell.py | 32 +- 2 files changed, 120 insertions(+), 213 deletions(-) diff --git a/h2integrate/converters/natural_gas/SO_NG_fuel_cell.py b/h2integrate/converters/natural_gas/SO_NG_fuel_cell.py index b7b347d70..1bdfbf9e6 100644 --- a/h2integrate/converters/natural_gas/SO_NG_fuel_cell.py +++ b/h2integrate/converters/natural_gas/SO_NG_fuel_cell.py @@ -3,6 +3,7 @@ from h2integrate.core.utilities import BaseConfig, merge_shared_inputs from h2integrate.core.validators import gte_zero +from h2integrate.tools.constants import H_MW, O2_MW, CH4_MW, CO2_MW, faraday from h2integrate.core.model_baseclasses import ( CostModelBaseClass, CostModelBaseConfig, @@ -30,81 +31,46 @@ class SONGFuelCellPerformanceConfig(BaseConfig): # fuel_cell_efficiency_hhv: float = field(validator=range_val(0, 1)) -def calc_current(power_ref, cell_area, n_cells, stack_number): +def calc_current(system_power_reference, cell_area, n_cells, n_stacks): + """_summary_ + + Args: + system_power_reference (np.ndarray): power demanded of the entire system in W + cell_area (float): cell active area in cm^2 + n_cells (int): number of cells per stack + n_stacks (int): number of stacks in the system + + Returns: + tuple(np.ndarray, np.poly1d): stack current and + function to convert from current density to voltage + """ # Calculates the current and voltage from IV curve based on power reference - current_curve = [ - 0.0356, - 0.05413333, - 0.0796, - 0.11366667, - 0.244, - 0.454, - # 0.70366667, - # 0.96933333, - # 1.24, - # 1.52666667, - # 1.80333333, - # 2.07, - # 2.32, - # 2.54333333, - # 2.73666667, - # 2.9, - ] # in A - voltage_curve = [ - 0.987, - 0.936, - 0.884, - 0.838, - 0.786, - 0.736, - # 0.686, - # 0.636, - # 0.586, - # 0.53566667, - # 0.486, - # 0.436, - # 0.386, - # 0.33533333, - # 0.286, - # 0.236, - ] # in V - power_curve = [ - 35.16666667, - 50.53333333, - 70.33333333, - 95.46666667, - 191.66666667, - 334.33333333, - # 482.66666667, - # 616.66666667, - # 729.0, - # 817.0, - # 875.33333333, - # 902.0, - # 895.0, - # 854.33333333, - # 782.66666667, - # 684.33333333, - ] - - # Change power from mW to W - power_curve = [x / 1e3 for x in power_curve] - - power_coefs = np.polyfit(power_curve, current_curve, 5) - power_I_curve = np.poly1d(power_coefs) - V_coefs = np.polyfit(current_curve, voltage_curve, 5) - V_I_curve = np.poly1d(V_coefs) - - # convert power_ref to Watts - power_ref = power_ref * 1e3 - # clip very high power_ref to max power in curve - power_density = power_ref / cell_area / stack_number / n_cells - power_density_clip = min(power_density, max(power_curve)) - - I_cell = max(power_I_curve(power_density_clip), 0) - V_cell = V_I_curve(I_cell) - I_cell = I_cell * cell_area - return I_cell, V_cell + J_curve = np.array([0.0356, 0.05413333, 0.0796, 0.11366667, 0.244, 0.454]) # in A/cm^2 + voltage_curve = np.array([0.987, 0.936, 0.884, 0.838, 0.786, 0.736]) # in V + power_curve = ( + np.array([35.16666667, 50.53333333, 70.33333333, 95.46666667, 191.66666667, 334.33333333]) + / 1e3 + ) # in W/cm^2 + + # function to calculate voltage from current density + V_coefs = np.polyfit(J_curve, voltage_curve, 5) + V_J_curve = np.poly1d(V_coefs) + + # function to calculate current density from power + # I_curve = J_curve * cell_area + # P_curve = I_curve * (n_cells * voltage_curve) + stack_P_curve = power_curve * cell_area * n_cells + J_coefs = np.polyfit(stack_P_curve, J_curve, 5) + J_P_curve = np.poly1d(J_coefs) + + # Calculate power per stack and power density + power_per_stack = system_power_reference / n_stacks + + stack_current_density = J_P_curve(power_per_stack) # convert to kW for the curve + stack_current = stack_current_density * cell_area * n_cells # in A + stack_current = np.clip(stack_current, a_min=0.0, a_max=None) # clip negative values + + return stack_current, V_J_curve class SONGFuelCellPerformanceModel(PerformanceModelBaseClass): @@ -233,12 +199,17 @@ def compute(self, inputs, outputs): outputs: OpenMDAO outputs object for electricity_out, natural_gas_consumed, oxygen_consumed, water_out, and carbon_dioxide_out. """ + # Set calculation constants: + M_H2 = H_MW * 2 / 1000 # Molar mass of H2 in kg/mol + M_O2 = O2_MW / 1000 # Molar mass of O2 in kg/mol + M_H2O = M_H2 + M_O2 / 2 # Molar mass of H2O in kg/mol + M_CH4 = CH4_MW / 1000 # Molar mass of CH4 in kg/mol + M_CO2 = CO2_MW / 1000 # Molar mass of CO2 in kg/mol # calculate max input and output - inputs["system_capacity"] # plant capacity in kW - natural_gas_in = inputs["natural_gas_in"] # MMBtu/h - oxygen_in = inputs["oxygen_in"] # kg/h + inputs["oxygen_in"] # kg/h stack_temperature = inputs["stack_temperature"] + natural_gas_in = inputs["natural_gas_in"] # MMBtu/h # fuel_cell_efficiency = inputs["fuel_cell_efficiency"] ############################################################################ @@ -248,136 +219,78 @@ def compute(self, inputs, outputs): # Convert natural gas input from MMBtu/h to kg/h for CH4 # Assume 1 MMBtu ~ 1 MCF # Use ideal gas law to account for temperature difference - stack_temperature = inputs["stack_temperature"] ng_density_cold = 0.717 # kg/m^3 at 25 deg C and 1 atm ng_density_hot = ng_density_cold * (298.15 / stack_temperature) # kg/m^3 at stack temp # Use Mass = volume * density, and 1 MCF = 1e3 ft^3 = 28.3168 m^3 mcf_to_kg = 28.3168 * ng_density_hot # kg - print(natural_gas_in, mcf_to_kg) natural_gas_kg_hr = natural_gas_in * mcf_to_kg # Convert MMBtu/h to kg/h for CH4 # Add consumption of water for steam reforming of natural gas to hydrogen - # Set calculation constants: - self.f_c = 96485.33 # Faraday's constant in A/mol - self.M_H2 = 0.002016 # Molar mass of H2 in kg/mol - self.M_O2 = 0.032 # Molar mass of O2 in kg/mol - self.M_H2O = 0.018 # Molar mass of H2O in kg/mol - self.M_CO2 = 0.044 # Molar mass of CO2 in kg/mol - self.M_CH4 = 0.01604 # Molar mass of CH4 (methane) in kg/mol - self.Tref = 298.15 # Standard room temperature in K [25 deg Celsius] - self.cp_H2 = 14300 # Specific heat of H2 in J/(kg*K) - self.cp_air = 1005 # Specific heat of air in J/(kg*K) - self.cp_H2O = 4184 # Specific heat of water in J/(kg*K) - self.cp_N2 = 1040 # Specific heat of nitrogen in J/(kg*K) - self.cp_O2 = 918 # Specific heat of oxygen in J/(kg*K) - self.hhv_h2 = 141.8 * 1e6 # Higher heating value of hydrogen in J/kg - self.hhv_air = 0 # No higher heating value of air - self.hhv_H2O = 2260 # Higher heating value of water in J/kg - self.hhv_CO2 = 0 # No higher heating value of CO2 # Sizing the cells - self.max_cell_power_density = 0.000334 - # is n_cells = N_series? - self.N_series = 1 - self.stack_size = self.config.system_capacity_kw / self.config.n_stacks - self.cell_active_area = 400 # [cm^2] from Battelle (https://www.energy.gov/sites/prod/files/2018/02/f49/fcto_battelle_mfg_cost_analysis_1%20_to_25kw_pp_chp_fc_systems_jan2017_0.pdf) - self.n_cells = round( - self.stack_size / (self.cell_active_area * self.max_cell_power_density) + max_cell_power_density = 0.000334 # in W/cm^2 + stack_size = inputs["system_capacity"][0] / self.config.n_stacks + cell_active_area = 400 # [cm^2] from Battelle (https://www.energy.gov/sites/prod/files/2018/02/f49/fcto_battelle_mfg_cost_analysis_1%20_to_25kw_pp_chp_fc_systems_jan2017_0.pdf) + n_cells = round(stack_size / (cell_active_area * max_cell_power_density)) + # Recalculate the rated power production based on final fuel cell sizing + rated_power_production = ( + max_cell_power_density * n_cells * cell_active_area * self.config.n_stacks ) - # PSUEDO CODE: - """ - 1. Receive power setpoint into fuel cell - 2. Find current and voltage from I-V curve based on power setpoint - 3. Calculate natural gas (CH4) and oxygen consumed based on Faraday's law - 4. Calculate electricity produced from cell voltage and current - 5. Calculate H2O and CO2 generated from the reaction - 6. Output electricity clipped to system capacity and other metrics + ################## Model Calculations ################## + # 1. Receive power setpoint into fuel cell + power_reference = np.clip( + inputs[f"{self.commodity}_command_value"], a_min=0.0, a_max=rated_power_production + ) - """ + # 2. Find stack current from power reference + commanded_I_stack, V_J_curve = calc_current( + power_reference * 1e3, cell_active_area, n_cells, self.config.n_stacks + ) - ng_consumed = np.zeros(self.n_timesteps) - o2_consumed = np.zeros(self.n_timesteps) - co2_generated = np.zeros(self.n_timesteps) - h2o_generated = np.zeros(self.n_timesteps) - commodity_out = np.zeros(self.n_timesteps) - - for i in range(self.n_timesteps): - power_reference = inputs[f"{self.commodity}_command_value"][i] - natural_gas_in_loop = natural_gas_kg_hr[i] - oxygen_in_loop = oxygen_in[i] - - # Find current and voltage from IV curve with power setpoint - I_cell, V_cell = calc_current( - power_reference, self.cell_active_area, self.n_cells, self.config.n_stacks - ) - - # Calculate natural gas and oxygen consumed - NG_consumed_rate = ((I_cell * self.N_series * self.M_CH4) / (8.0 * self.f_c)) * ( - self.dt * self.config.n_stacks * self.n_cells - ) # kg/time step - O2_consumed_rate = ((I_cell * self.N_series * self.M_O2) / (4.0 * self.f_c)) * ( - self.dt * self.config.n_stacks * self.n_cells - ) # kg/time step - - if NG_consumed_rate > natural_gas_in_loop or O2_consumed_rate > oxygen_in_loop: - # implement an adjustment based on H2 & O2 available - new_i_ng = ( - natural_gas_in_loop - / (self.dt * self.config.n_stacks * self.n_cells) - * (8.0 * self.f_c) - / (self.N_series * self.M_CH4) - ) - new_i_o2 = ( - oxygen_in_loop - / (self.dt * self.config.n_stacks * self.n_cells) - * (4.0 * self.f_c) - / (self.N_series * self.M_O2) - ) - I_cell = min(new_i_ng, new_i_o2) - # TODO: recalc voltage based on new current - print("Not enough NG or O2 for this power point") - # Calculate natural gas and oxygen consumed - NG_consumed_rate = ((I_cell * self.N_series * self.M_CH4) / (8.0 * self.f_c)) * ( - self.dt * self.config.n_stacks * self.n_cells - ) # kg/time step - O2_consumed_rate = ((I_cell * self.N_series * self.M_O2) / (4.0 * self.f_c)) * ( - self.dt * self.config.n_stacks * self.n_cells - ) # kg/time step - - # Compute electricity from the system - electricity_produced = ( - V_cell * I_cell * self.n_cells * self.config.n_stacks / 1e3 - ) # Calculated in watts, convert to kW - - # Compute H2O out - H2O_generated = ( - (I_cell * self.N_series / (2 * self.f_c)) - * self.M_H2O - * (self.dt * self.config.n_stacks * self.n_cells) - ) # in kg/time step - # TODO: check electron numbers in these calculations - # Compute CO2 out - CO2_generated = ( - (I_cell * self.N_series / (8 * self.f_c)) - * self.M_CO2 - * (self.dt * self.config.n_stacks * self.n_cells) - ) # in kg/time step - - ng_consumed[i] = NG_consumed_rate - o2_consumed[i] = O2_consumed_rate - h2o_generated[i] = H2O_generated - co2_generated[i] = CO2_generated - commodity_out[i] = electricity_produced + # 3. Find available hydrogen and oxygen for each timestep + ng_in_kg_per_s = natural_gas_kg_hr / 3600 # convert from kg/h to kg/s + o2_in_kg_per_s = inputs["oxygen_in"] / 3600 # convert from kg/h to kg/s + + # convert from kg/s to A per stack - current that feedstocks in can support + I_stack_from_ng = (ng_in_kg_per_s * 8 * faraday) / (M_CH4 * self.config.n_stacks) + I_stack_from_o2 = (o2_in_kg_per_s * 4 * faraday) / (M_O2 * self.config.n_stacks) + + # 4. Take minimum current from power reference, hydrogen available, and oxygen available + I_stack = np.minimum(commanded_I_stack, np.minimum(I_stack_from_ng, I_stack_from_o2)) + + # 5. Calculate current density and voltage from I-V curve + J_cell = I_stack / (cell_active_area * n_cells) # in A/cm^2 + I_cell = J_cell * cell_active_area # in A + V_cell = V_J_curve(J_cell) # in V + + # 6. Calculate power output from current and voltage + power_out = V_cell * I_cell * n_cells * self.config.n_stacks / 1e3 # in kW + + # 7. Calculate hydrogen and oxygen consumed and water produced + # based on electrochemical reactions + ng_consumed = ((I_cell * M_CH4) / (8.0 * faraday)) * ( + self.dt * self.config.n_stacks * n_cells + ) # kg/time step + o2_consumed = ((I_cell * M_O2) / (4.0 * faraday)) * ( + self.dt * self.config.n_stacks * n_cells + ) # kg/time step + h2o_generated = ((I_cell * M_H2O) / (2 * faraday)) * ( + self.dt * self.config.n_stacks * n_cells + ) # kg/time step + + co2_generated = ((I_cell * M_CO2) / (8 * faraday)) * ( + self.dt * self.config.n_stacks * n_cells + ) # kg/time step # Set Outputs # clip the electricity output to the system capacity - outputs["electricity_out"] = np.minimum(commodity_out, self.config.system_capacity_kw) + outputs["rated_electricity_production"] = rated_power_production + outputs["electricity_out"] = np.minimum(power_out, rated_power_production) outputs["total_electricity_produced"] = np.sum(outputs["electricity_out"]) * ( self.dt / 3600 ) - outputs["rated_electricity_production"] = self.config.system_capacity_kw outputs["annual_electricity_produced"] = outputs["total_electricity_produced"] * ( 1 / self.fraction_of_year_simulated ) @@ -397,14 +310,20 @@ def compute(self, inputs, outputs): class SONGFuelCellCostConfig(CostModelBaseConfig): """Configuration class for the solid oxide natural gas fuel cell cost model. - Fields include `system_capacity_kw`, `capex_stack_per_kw`, - `capex_bop`, `capex_indirect_costs_per_kw`, and `fixed_opex_per_kw_per_year`. + Attributes: + system_capacity_kw (float): The capacity of the fuel cell system in kilowatts (kW). + capex_stack_per_kw (float): Capital expenditure for the stack per kilowatt ($/kW). + capex_bop (float): Capital expenditure for balance of plant ($/kW). + capex_indirect_costs_per_kw (float): Indirect capital costs per kilowatt ($/kW). + fixed_opex_per_kw_per_year (float): Fixed operating expenditure per + kilowatt per year ($/kW/year). + The `cost_year` field is inherited from `CostModelBaseConfig`. """ system_capacity_kw: float = field(validator=gte_zero) capex_stack_per_kw: float = field(validator=gte_zero) - capex_bop: float = field(validator=gte_zero) + capex_bop_per_kw: float = field(validator=gte_zero) capex_indirect_costs_per_kw: float = field(validator=gte_zero) fixed_opex_per_kw_per_year: float = field(validator=gte_zero) @@ -440,7 +359,7 @@ def setup(self): self.add_input( "unit_capex", val=self.config.capex_stack_per_kw - + self.config.capex_bop + + self.config.capex_bop_per_kw + self.config.capex_indirect_costs_per_kw, units="USD/kW", desc="Capital cost per unit capacity", diff --git a/h2integrate/converters/natural_gas/test/test_SO_NG_fuel_cell.py b/h2integrate/converters/natural_gas/test/test_SO_NG_fuel_cell.py index b713cc323..e919e7bff 100644 --- a/h2integrate/converters/natural_gas/test/test_SO_NG_fuel_cell.py +++ b/h2integrate/converters/natural_gas/test/test_SO_NG_fuel_cell.py @@ -43,14 +43,9 @@ def cost_config(): "model_inputs": { "cost_parameters": { "system_capacity_kw": 1500.0, - "capex_stack_per_kw": 800.0, - "capex_fuel_supply_per_kw": 140.67, - "capex_air_supply_per_kw": 138.79, - "capex_cooling_per_kw": 64.14, - "capex_controls_instrumentation_per_kw": 106.52, - "capex_electrical_per_kw": 447.83, - "capex_assembly_per_kw": 63.51, - "capex_additional_labor_per_kw": 137.96, + "capex_stack_per_kw": 500.0, + "capex_bop_per_kw": 1000, + "capex_indirect_costs_per_kw": 1500, "fixed_opex_per_kw_per_year": 31.0, "cost_year": 2026, } @@ -112,9 +107,9 @@ def test_fuel_cell_performance(tech_config, plant_config, subtests): with subtests.test("rated_electricity_production"): assert ( pytest.approx( - prob.get_val("fuel_cell.rated_electricity_production", units="kW"), rel=1e-6 + prob.get_val("fuel_cell.rated_electricity_production", units="kW"), rel=1e-4 ) - == 1500.0 + == 1498.9 ) with subtests.test("total_electricity_produced"): @@ -136,13 +131,13 @@ def test_fuel_cell_performance(tech_config, plant_config, subtests): with subtests.test("oxygen consumed"): assert ( pytest.approx(np.sum(prob.get_val("fuel_cell.oxygen_consumed", units="kg/h")), rel=1e-6) - == 18186.36156 + == 18185.2245185 ) with subtests.test("water out"): assert ( pytest.approx(np.sum(prob.get_val("fuel_cell.water_out", units="kg/h")), rel=1e-6) - == 20459.65675 + == 20476.7060254 ) with subtests.test("carbon dioxide out"): @@ -150,7 +145,7 @@ def test_fuel_cell_performance(tech_config, plant_config, subtests): pytest.approx( np.sum(prob.get_val("fuel_cell.carbon_dioxide_out", units="kg/h")), rel=1e-6 ) - == 12503.12357 + == 12505.9649206 ) @@ -259,7 +254,7 @@ def test_fuel_cell_demand(tech_config, plant_config, subtests): assert co2_out[5] == pytest.approx(0.0, abs=1e-6) with subtests.test("limited NG supply reduces electricity output and byproducts"): - assert electricity_output[6] == pytest.approx(55.547789, rel=1e-2) + assert electricity_output[6] == pytest.approx(79.43635323, rel=1e-2) assert ng_consumed[6] == pytest.approx(1.0, rel=1e-2) assert o2_consumed[6] == pytest.approx(22.509924, rel=1e-2) assert water_out[6] == pytest.approx(25.323664, rel=1e-2) @@ -289,14 +284,7 @@ def test_fuel_cell_cost(cost_config, plant_config, subtests): cp = cost_config["model_inputs"]["cost_parameters"] expected_unit_capex = ( - cp["capex_stack_per_kw"] - + cp["capex_fuel_supply_per_kw"] - + cp["capex_air_supply_per_kw"] - + cp["capex_cooling_per_kw"] - + cp["capex_controls_instrumentation_per_kw"] - + cp["capex_electrical_per_kw"] - + cp["capex_assembly_per_kw"] - + cp["capex_additional_labor_per_kw"] + cp["capex_stack_per_kw"] + cp["capex_bop_per_kw"] + cp["capex_indirect_costs_per_kw"] ) expected_capex = cp["system_capacity_kw"] * expected_unit_capex expected_opex = cp["system_capacity_kw"] * cp["fixed_opex_per_kw_per_year"]