diff --git a/CHANGELOG.md b/CHANGELOG.md index 79e906e30..da312bc6c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -86,6 +86,7 @@ - Adds `feedstock_dir` to the EIA natural gas retrieval to align the downloading or loading of the feedstock data with the resource data methodology [PR 801](https://github.com/NatLabRockies/H2Integrate/pull/801). - Add support for slice notation in technology connections to allow users to connect between variables of different shapes. [PR 774](https://github.com/NatLabRockies/H2Integrate/pull/774) - Updated edge attribute `commodity` of in `H2Integrate.create_technology_graph` to use lists instead of strings to account for systems with multiple commodities connected between two technologies [PR 823](https://github.com/NatLabRockies/H2Integrate/pull/823) +- Added two new fuel cell models: `PEMH2FuelCellPerformanceModel` to model a PEM hydrogen fuel cell and `SONGFuelCellPerformanceModel` to model a natural gas solid oxide fuel cell [PR 794](https://github.com/NatLabRockies/H2Integrate/pull/794) ### Fixes - Bug fix so multi-level output path won't throw an error; updated test for EIA API handling. [PR 820](https://github.com/NatLabRockies/H2Integrate/pull/820) diff --git a/docs/_toc.yml b/docs/_toc.yml index e9b00f934..71c0548fc 100644 --- a/docs/_toc.yml +++ b/docs/_toc.yml @@ -51,6 +51,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_h2_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..0272c9646 --- /dev/null +++ b/docs/technology_models/PEM_h2_fuel_cell.md @@ -0,0 +1,32 @@ +# 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 IV curve was extracted from the data in this GitHub repository [1](https://github.com/ECSIM/pem-dataset1). The model uses chemical reaction equations to calculate the consumed and produced commodities from the fuel cell reaction ([3](https://www.sciencedirect.com/science/article/pii/S0360319906005726), [4](https://www.sciencedirect.com/science/article/pii/S0360319924051577)). + +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²)([2](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. + + +```{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: +``` 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..1d384b058 --- /dev/null +++ b/docs/technology_models/SO_NG_fuel_cell.md @@ -0,0 +1,31 @@ +# 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 IV curve was extracted from the data in this GitHub repository [1](https://github.com/ECSIM/pem-dataset1). The model uses chemical reaction equations to calculate the consumed and produced commodities from the fuel cell reaction ([3](https://www.sciencedirect.com/science/article/pii/S0360319906005726), [4](https://www.sciencedirect.com/science/article/pii/S0360319924051577)). + +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 ([2](https://www.sciencedirect.com/science/article/pii/S0360319906005726)). 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. 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. + +```{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: +``` diff --git a/docs/user_guide/model_overview.md b/docs/user_guide/model_overview.md index 1d54273b1..a24b031a8 100644 --- a/docs/user_guide/model_overview.md +++ b/docs/user_guide/model_overview.md @@ -155,6 +155,7 @@ auto-generated API page. + {py:class}`~h2integrate.converters.hydrogen.pem_electrolyzer.ECOElectrolyzerPerformanceModel` - An OpenMDAO component that wraps the PEM electrolyzer model. + {py:class}`~h2integrate.converters.hydrogen.htse_electrolyzer.HTSEPerformanceModel` - A simplified high-temperature steam electrolysis (HTSE) model. + {py:class}`~h2integrate.converters.hydrogen.h2_fuel_cell.LinearH2FuelCellPerformanceModel` - Performance model for a hydrogen fuel cell. + + {py:class}`~h2integrate.converters.hydrogen.PEM_h2_fuel_cell.PEMH2FuelCellPerformanceModel` - Performance model for a PEM hydrogen fuel cell. + {py:class}`~h2integrate.converters.hydrogen.steam_methane_reformer.SteamMethaneReformerPerformanceModel` - Performance model for steam methane reforming (SMR) hydrogen production plants. - cost models: + {py:class}`~h2integrate.converters.hydrogen.basic_cost_model.BasicElectrolyzerCostModel` - An OpenMDAO component that computes the cost of a PEM electrolyzer. @@ -197,6 +198,7 @@ auto-generated API page. - `natural_gas`: natural gas combined-cycle and combustion turbines - performance models: + {py:class}`~h2integrate.converters.natural_gas.natural_gas_cc_ct.NaturalGasPerformanceModel` - Performance model for natural gas power plants. + + {py:class}`~h2integrate.converters.natural_gas.SO_NG_fuel_cell.SONGFuelCellPerformanceModel` - Performance model for a solid oxide natural gas fuel cell. - cost models: + {py:class}`~h2integrate.converters.natural_gas.natural_gas_cc_ct.NaturalGasCostModel` - Cost model for natural gas power plants. 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/plant_config.yaml b/examples/37_pem_fc/plant_config.yaml new file mode 100644 index 000000000..1237e0198 --- /dev/null +++ b/examples/37_pem_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 + - [h2_feedstock, PEM_fuel_cell, hydrogen, pipe] + - [o2_feedstock, PEM_fuel_cell, oxygen, pipe] + # connect fuel cell to demand component + - [PEM_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: PEM_fuel_cell + technologies: + - 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 new file mode 100644 index 000000000..daa9ce6f3 --- /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) * 1000 +model.prob.set_val("PEM_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..619a18180 --- /dev/null +++ b/examples/37_pem_fc/tech_config.yaml @@ -0,0 +1,62 @@ +name: technology_config +description: This hybrid plant produces ammonia +technologies: + PEM_fuel_cell: + performance_model: + model: PEMH2FuelCellPerformanceModel + cost_model: + model: GenericConverterCostModel + model_inputs: + performance_parameters: + system_capacity_kw: 1500 + n_stacks: 60 + stack_temperature_K: 278 # Not used yet + hhv: 39.4 # kWh/kg + cost_parameters: + commodity: electricity + commodity_rate_units: kW + commodity_amount_units: kW*h + unit_capex: 1999.6 # USD/kW + unit_opex: 31 # USD/kW/year + unit_varopex: 0.0 # USD/kWh + cost_year: 2026 + electricity_load_demand: + performance_model: + model: GenericDemandComponent + model_inputs: + performance_parameters: + commodity: electricity + commodity_rate_units: kW + demand_profile: 1000 + 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: 400.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/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..549744544 --- /dev/null +++ b/examples/38_song_fc/tech_config.yaml @@ -0,0 +1,62 @@ +name: technology_config +description: This hybrid plant produces ammonia +technologies: + SONG_fuel_cell: + performance_model: + model: SONGFuelCellPerformanceModel + cost_model: + model: GenericConverterCostModel + model_inputs: + performance_parameters: + system_capacity_kw: 1500 + n_stacks: 60 + stack_temperature_K: 1073 + hhv: 15.4 # kWh/kg --> based on methane + cost_parameters: + commodity: electricity + commodity_rate_units: kW + commodity_amount_units: kW*h + unit_capex: 3174 # USD/kW + unit_opex: 35 # USD/kW/year + unit_varopex: 0.0 # USD/kWh + 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: 30.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: 650.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/hydrogen/PEM_h2_fuel_cell.py b/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py new file mode 100644 index 000000000..238c30e01 --- /dev/null +++ b/h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py @@ -0,0 +1,345 @@ +import numpy as np +from attrs import field, define, validators + +from h2integrate.core.utilities import BaseConfig, merge_shared_inputs +from h2integrate.tools.constants import H_MW, O2_MW, faraday +from h2integrate.core.model_baseclasses import 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). + 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). + """ + + system_capacity_kw: float = field(validator=validators.ge(0)) + n_stacks: int = field(validator=validators.ge(0)) + stack_temperature_K: float = field(validator=validators.ge(0)) + hhv: float = field(validator=validators.ge(0)) + # min_system_power_fraction_kw: float + # fuel_cell_efficiency_hhv: float = field(validator=range_val(0, 1)) + + +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 + # These current, voltage and power values are from the fuel cell data collected here: https://github.com/ECSIM/pem-dataset1 + # Using the data from the "Activation Test MEA Standard Protocol (Repeat)" case + 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 + 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 # in Watts + + # Create power/current density relationship curve + stack_current_density = J_P_curve(power_per_stack) + stack_current = stack_current_density * cell_area # 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): + """ + 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 = ( + 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 = PEMH2FuelCellPerformanceConfig.from_dict( + merge_shared_inputs(self.options["tech_config"]["model_inputs"], "performance"), + additional_cls_name=self.__class__.__name__, + ) + + 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", + ) + + # 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="Rated electricity production of the PEM fuel cell system", + ) + + self.add_output( + "hydrogen_consumed", + val=0.0, + shape=self.n_timesteps, + units=f"kg/({self.dt}*s)", + desc="Mass flow rate of hydrogen consumed by the fuel cell", + ) + + self.add_output( + "oxygen_consumed", + val=0.0, + shape=self.n_timesteps, + units=f"kg/({self.dt}*s)", + desc="Mass flow rate of oxygen consumed by the fuel cell", + ) + + self.add_output( + "water_out", + val=0.0, + shape=self.n_timesteps, + units=f"kg/({self.dt}*s)", + desc="Mass flow rate of water produced by the fuel cell", + ) + + self.add_output( + "heat_out", + val=0.0, + shape=self.n_timesteps, + units=self.commodity_rate_units, + desc="Heat generated by the fuel cell", + ) + + self.add_output( + "rated_h2_consumed", + val=0.0, + units=f"kg/({self.dt}*s)", + desc="Rated hydrogen consumed by the fuel cell", + ) + + self.add_output( + "rated_o2_consumed", + val=0.0, + units=f"kg/({self.dt}*s)", + desc="Rated oxygen consumed by the fuel cell", + ) + + self.add_output( + "rated_water_production", + val=0.0, + units=f"kg/({self.dt}*s)", + desc="Rated water produced by the fuel cell", + ) + + self.add_output( + "rated_heat_production", + val=0.0, + units=self.commodity_rate_units, + desc="Rated heat generated by the fuel cell", + ) + + self.add_output( + "rated_stack_efficiency", + val=0.0, + desc="Rated stack efficiency of 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 PEM fuel cell", + ) + + def compute(self, inputs, outputs): + """ + 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, 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. + """ + + # 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 + # Electron transfer constants + n_h2 = 2 # number of electrons transferred per mole of H2 + n_o2 = 4 # number of electrons transferred per mole of O2 + n_h2o = 2 # number of electrons transferred per mole of H2O + + # Sizing the cells + max_cell_power_density = 0.000334 # in kW/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 + ) + + # Calculate the rated outputs of the system + rated_I_stack, _ = calc_current( + rated_power_production * 1e3, cell_active_area, n_cells, self.config.n_stacks + ) # in A per stack + rated_h2_consumed = ((rated_I_stack * M_H2 * n_cells) / (n_h2 * faraday)) * ( + self.dt * self.config.n_stacks + ) # kg/time step + rated_o2_consumed = ((rated_I_stack * M_O2 * n_cells) / (n_o2 * faraday)) * ( + self.dt * self.config.n_stacks + ) # kg/time step + rated_h2o_out = ((rated_I_stack * M_H2O * n_cells) / (n_h2o * faraday)) * ( + self.dt * self.config.n_stacks + ) # kg/time step + rated_heat_out = rated_h2_consumed * self.config.hhv - rated_power_production # in kW + + ################## 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 + ) # in commodity rate units (kW) + + # 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 + ) # in A + + # 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 * n_h2 * faraday) / ( + M_H2 * self.config.n_stacks * n_cells + ) + I_stack_from_o2 = (o2_in_kg_per_s * n_o2 * faraday) / ( + M_O2 * self.config.n_stacks * n_cells + ) + + # 4. Take minimum current from power reference, hydrogen available, and oxygen available + I_stack = np.minimum.reduce([commanded_I_stack, I_stack_from_h2, I_stack_from_o2]) + # in Amps per stack + + # 5. Calculate current density and voltage from I-V curve, all of these are per stack + J_cell = I_stack / (cell_active_area) # in A/cm^2 + V_cell = V_J_curve(J_cell) # in V + + # 6. Calculate power output from current and voltage + power_out = V_cell * I_stack * n_cells * self.config.n_stacks / 1e3 + # Calculated in Watts, then converted to kW + + # 7. Calculate hydrogen and oxygen consumed and water produced + # based on electrochemical reactions + h2_consumed = ((I_stack * M_H2 * n_cells) / (n_h2 * faraday)) * ( + self.dt * self.config.n_stacks + ) # kg/time step + o2_consumed = ((I_stack * M_O2 * n_cells) / (n_o2 * faraday)) * ( + self.dt * self.config.n_stacks + ) # kg/time step + h2o_generated = ((I_stack * M_H2O * n_cells) / (n_h2o * faraday)) * ( + self.dt * self.config.n_stacks + ) # kg/time step + + # 8. Calculate heat generation from fuel cell + heat_generated = h2_consumed * self.config.hhv - power_out # in kW + + # Set Outputs + # clip the electricity output to the system capacity + 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["annual_electricity_produced"] = outputs["total_electricity_produced"] * ( + 1 / self.fraction_of_year_simulated + ) + outputs["capacity_factor"] = outputs["total_electricity_produced"] / ( + rated_power_production * self.n_timesteps * (self.dt / 3600) + ) + outputs["hydrogen_consumed"] = h2_consumed + outputs["oxygen_consumed"] = o2_consumed + outputs["water_out"] = h2o_generated + outputs["heat_out"] = heat_generated + + # Set rated outputs + outputs["rated_h2_consumed"] = rated_h2_consumed + outputs["rated_o2_consumed"] = rated_o2_consumed + outputs["rated_water_production"] = rated_h2o_out + outputs["rated_heat_production"] = rated_heat_out + outputs["rated_stack_efficiency"] = rated_power_production / ( + rated_heat_out + rated_power_production + ) + + # TODO: implement a hydrogen and oxygen conversion efficiency based on stack + # temperature and other factors diff --git a/h2integrate/converters/hydrogen/__init__.py b/h2integrate/converters/hydrogen/__init__.py index 0bf80327d..cd0e54670 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, 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..68d5a4fc3 --- /dev/null +++ b/h2integrate/converters/hydrogen/test/test_PEM_fuel_cell.py @@ -0,0 +1,260 @@ +import numpy as np +import pytest +import openmdao.api as om +from pytest import fixture + +from h2integrate.converters.hydrogen.PEM_h2_fuel_cell import PEMH2FuelCellPerformanceModel + + +@fixture +def plant_config(): + plant_config = { + "plant": { + "plant_life": 30, + "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": 278.0, + "hhv": 39.4, # kWh/kg + } + } + } + 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.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") + 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 + + with subtests.test("electricity output is non-negative"): + assert np.min(electricity_output) >= 0.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("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"): + 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-4 + ) + == 1498.9 + ) + + 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) + == 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) + == 20476.70602546 + ) + + with subtests.test("rated h2 consumed"): + assert ( + pytest.approx(prob.get_val("fuel_cell.rated_h2_consumed", units="kg/h"), rel=1e-6) + == 76.5012465 + ) + + with subtests.test("rated o2 consumed"): + assert ( + pytest.approx(prob.get_val("fuel_cell.rated_o2_consumed", units="kg/h"), rel=1e-6) + == 607.11480 + ) + + with subtests.test("rated water out"): + assert ( + pytest.approx(prob.get_val("fuel_cell.rated_water_production", units="kg/h"), rel=1e-6) + == 683.6160498 + ) + + with subtests.test("rated heat out"): + assert ( + pytest.approx(prob.get_val("fuel_cell.rated_heat_production", units="kW"), rel=1e-6) + == 1515.1571138 + ) + + with subtests.test("rated stack efficiency"): + assert ( + pytest.approx(prob.get_val("fuel_cell.rated_stack_efficiency"), rel=1e-4) + == 0.49732 # Example calculation + ) + + +@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") + heat_out = prob.get_val("fuel_cell.heat_out", units="kW") + + with subtests.test("output clipped to system capacity"): + 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 + 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] < 30.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.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.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.61605, 204.817867, 683.61605, 0.0, 8.936012, 0.0] + np.testing.assert_allclose(water_out[:6], expected_water_out, rtol=1e-4) + + with subtests.test("heat out"): + expected_heat_out = [1513.436595, 404.4013, 1513.436595, 0.0, 10.155, 0.0] + np.testing.assert_allclose(heat_out[:6], expected_heat_out, rtol=1e-4) 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..700f70add --- /dev/null +++ b/h2integrate/converters/natural_gas/SO_NG_fuel_cell.py @@ -0,0 +1,391 @@ +import numpy as np +from attrs import field, define, validators +from openmdao.utils import units + +from h2integrate.core.utilities import BaseConfig, merge_shared_inputs +from h2integrate.tools.constants import H_MW, O2_MW, CH4_MW, CO2_MW, LHV_CH4_MJ_PER_KG, faraday +from h2integrate.core.model_baseclasses import PerformanceModelBaseClass + + +@define(kw_only=True) +class SONGFuelCellPerformanceConfig(BaseConfig): + """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). + 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). + hhv (float): higher heating value of the fuel in kWh/kg + """ + + system_capacity_kw: float = field(validator=validators.ge(0)) + n_stacks: int = field(validator=validators.ge(0)) + stack_temperature_K: float = field(validator=validators.ge(0)) + hhv: float = field(validator=validators.ge(0)) + # min_system_power_fraction_kw: float + # fuel_cell_efficiency_hhv: float = field(validator=range_val(0, 1)) + + +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 + # These current, voltage and power values are from the fuel cell data collected here: https://github.com/ECSIM/pem-dataset1 + # Using the data from the "Activation Test MEA Standard Protocol (Repeat)" case + 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 + 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 # in Watts + + # Create power/current density relationship curve + stack_current_density = J_P_curve(power_per_stack) + stack_current = stack_current_density * cell_area # 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): + """ + Performance model for a solid oxide natural gas fuel cell. + + 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 + - 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 + - co2_out is the mass flow rate of carbon dioxide produced in kg/hr + """ + + _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 = SONGFuelCellPerformanceConfig.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="MMBtu/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", + ) + + # 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="Rated electricity production of the SOFC system", + ) + + self.add_output( + "natural_gas_consumed", + val=0.0, + shape=self.n_timesteps, + units="MMBtu/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=f"kg/({self.dt}*s)", + desc="Mass flow rate of oxygen consumed by the fuel cell", + ) + + self.add_output( + "water_out", + val=0.0, + shape=self.n_timesteps, + units=f"kg/({self.dt}*s)", + desc="Mass flow rate of water produced by the fuel cell", + ) + + self.add_output( + "co2_out", + val=0.0, + shape=self.n_timesteps, + units=f"kg/({self.dt}*s)", + desc="Mass flow rate of carbon dioxide produced by the fuel cell", + ) + + self.add_output( + "heat_out", + val=0.0, + shape=self.n_timesteps, + units=self.commodity_rate_units, + desc="Heat generated by the fuel cell", + ) + + self.add_output( + "rated_natural_gas_consumed", + val=0.0, + units="MMBtu/h", + desc="Rated natural gas consumed by the fuel cell", + ) + + self.add_output( + "rated_oxygen_consumed", + val=0.0, + units=f"kg/({self.dt}*s)", + desc="Rated oxygen consumed by the fuel cell", + ) + + self.add_output( + "rated_water_production", + val=0.0, + units=f"kg/({self.dt}*s)", + desc="Rated water produced by the fuel cell", + ) + + self.add_output( + "rated_co2_production", + val=0.0, + units=f"kg/({self.dt}*s)", + desc="Rated carbon dioxide produced by the fuel cell", + ) + + self.add_output( + "rated_heat_production", + val=0.0, + units=self.commodity_rate_units, + desc="Rated heat generated by the fuel cell", + ) + + self.add_output( + "rated_stack_efficiency", + val=0.0, + desc="Rated stack efficiency of 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 SOFC based on natural gas input, + oxygen availability, and fuel cell electrochemical reactions. + + Args: + 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 co2_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 + # Electron transfer constants + n_ng = 8 # number of electrons transferred per mole of CH4 + n_o2 = 4 # number of electrons transferred per mole of O2 + n_h2o = 2 # number of electrons transferred per mole of H2O + n_co2 = 8 # number of electrons transferred per mole of CO2 + + # calculate max input and output + inputs["stack_temperature"] + natural_gas_in = inputs["natural_gas_in"] # MMBtu/h + + ############################################################################ + # Convert from MMBtu/h to kg/h for CH4 using LHV + natural_gas_in_MJ_per_hr = units.convert_units(natural_gas_in, "MMBtu/h", "MJ/h") + natural_gas_kg_hr = natural_gas_in_MJ_per_hr / LHV_CH4_MJ_PER_KG + # Convert MJ/h to kg/h for CH4 + + # TODO: Add consumption of water for steam reforming of natural gas to hydrogen + + # Sizing the cells + max_cell_power_density = 0.000334 # in kW/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 + ) + + # Calculate the rated outputs of the system + rated_I_stack, _ = calc_current( + rated_power_production * 1e3, cell_active_area, n_cells, self.config.n_stacks + ) # in A per stack + rated_ng_consumed = ((rated_I_stack * M_CH4 * n_cells) / (n_ng * faraday)) * ( + self.dt * self.config.n_stacks + ) # kg/time step + rated_o2_consumed = ((rated_I_stack * M_O2 * n_cells) / (n_o2 * faraday)) * ( + self.dt * self.config.n_stacks + ) # kg/time step + rated_h2o_out = ((rated_I_stack * M_H2O * n_cells) / (n_h2o * faraday)) * ( + self.dt * self.config.n_stacks + ) # kg/time step + rated_co2_out = ((rated_I_stack * M_CO2 * n_cells) / (n_co2 * faraday)) * ( + self.dt * self.config.n_stacks + ) # kg/time step + # Convert rated natural gas from kg/time step to MMBtu/h for output units + rated_ng_consumed_MJ_per_hr = ( + rated_ng_consumed * LHV_CH4_MJ_PER_KG * (3600 / self.dt) + ) # Convert kg/time step to MJ/h + rated_ng_consumed_mmbtu_per_hr = units.convert_units( + rated_ng_consumed_MJ_per_hr, "MJ/h", "MMBtu/h" + ) # Convert MJ/h to MMBtu/h + rated_heat_out = ( + rated_ng_consumed * (3600 / self.dt) * self.config.hhv - rated_power_production + ) + # in kW + + ################## 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 + ) # in commodity rate units (kW) + + # 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 + ) # current in Amps per stack + + # 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 * n_ng * faraday) / ( + M_CH4 * self.config.n_stacks * n_cells + ) + I_stack_from_o2 = (o2_in_kg_per_s * n_o2 * faraday) / ( + M_O2 * self.config.n_stacks * n_cells + ) + + # 4. Take minimum current from power reference, hydrogen available, and oxygen available + I_stack = np.minimum.reduce([commanded_I_stack, I_stack_from_ng, I_stack_from_o2]) + # in Amps per stack + + # 5. Calculate current density and voltage from I-V curve, all of these are per stack + J_cell = I_stack / (cell_active_area) # in A/cm^2 + V_cell = V_J_curve(J_cell) # in V + + # 6. Calculate power output from current and voltage + power_out = V_cell * I_stack * n_cells * self.config.n_stacks / 1e3 + # Calculated in Watts, then converted to kW + + # 7. Calculate hydrogen and oxygen consumed and water produced + # based on electrochemical reactions + ng_consumed = ((I_stack * M_CH4 * n_cells) / (n_ng * faraday)) * ( + self.dt * self.config.n_stacks + ) # kg/time step + o2_consumed = ((I_stack * M_O2 * n_cells) / (n_o2 * faraday)) * ( + self.dt * self.config.n_stacks + ) # kg/time step + h2o_generated = ((I_stack * M_H2O * n_cells) / (n_h2o * faraday)) * ( + self.dt * self.config.n_stacks + ) # kg/time step + co2_generated = ((I_stack * M_CO2 * n_cells) / (n_co2 * faraday)) * ( + self.dt * self.config.n_stacks + ) # kg/time step + + # 8. Calculate heat generated by the fuel cell + heat_generated = ng_consumed * (3600 / self.dt) * self.config.hhv - power_out # in kW + + # Set Outputs + # clip the electricity output to the system capacity + 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["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) + ) + ng_consumed_MJ_per_hr = ( + ng_consumed * LHV_CH4_MJ_PER_KG * (3600 / self.dt) + ) # Convert kg/time step to MJ/h + ng_consumed_mmbtu_per_hr = units.convert_units( + ng_consumed_MJ_per_hr, "MJ/h", "MMBtu/h" + ) # Convert MJ/h to MMBtu/h + outputs["natural_gas_consumed"] = ng_consumed_mmbtu_per_hr # Convert back to MMBtu/h + outputs["oxygen_consumed"] = o2_consumed + outputs["water_out"] = h2o_generated + outputs["co2_out"] = co2_generated + outputs["heat_out"] = heat_generated + + # Set rated outputs + outputs["rated_natural_gas_consumed"] = rated_ng_consumed_mmbtu_per_hr + outputs["rated_oxygen_consumed"] = rated_o2_consumed + outputs["rated_water_production"] = rated_h2o_out + outputs["rated_co2_production"] = rated_co2_out + outputs["rated_heat_production"] = rated_heat_out + outputs["rated_stack_efficiency"] = rated_power_production / ( + rated_heat_out + rated_power_production + ) + + # TODO: implement a natural gas and oxygen conversion efficiency based on stack + # temperature and other factors diff --git a/h2integrate/converters/natural_gas/__init__.py b/h2integrate/converters/natural_gas/__init__.py index 50d2251a8..3dfb0567e 100644 --- a/h2integrate/converters/natural_gas/__init__.py +++ b/h2integrate/converters/natural_gas/__init__.py @@ -8,3 +8,6 @@ SimpleGasConsumerPerformance, SimpleGasConsumerCost, ) +from h2integrate.converters.natural_gas.SO_NG_fuel_cell import ( + SONGFuelCellPerformanceModel, +) 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..bd9c6140f --- /dev/null +++ b/h2integrate/converters/natural_gas/test/test_SO_NG_fuel_cell.py @@ -0,0 +1,285 @@ +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 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, + "hhv": 15.4, # kWh/kg --> based on methane + } + } + } + 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) * 10.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-4 + ) + == 1498.9 + ) + + 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 + ) + == 216.00547288 + ) + + with subtests.test("oxygen consumed"): + assert ( + pytest.approx(np.sum(prob.get_val("fuel_cell.oxygen_consumed", units="kg/h")), rel=1e-6) + == 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) + == 20476.7060254 + ) + + with subtests.test("co2 out"): + assert ( + pytest.approx(np.sum(prob.get_val("fuel_cell.co2_out", units="kg/h")), rel=1e-6) + == 12505.9649206 + ) + + with subtests.test("rated natural gas consumed"): + assert ( + pytest.approx( + prob.get_val("fuel_cell.rated_natural_gas_consumed", units="MMBtu/h"), rel=1e-6 + ) + == 7.2113556 + ) + + with subtests.test("rated oxygen consumed"): + assert ( + pytest.approx(prob.get_val("fuel_cell.rated_oxygen_consumed", units="kg/h"), rel=1e-6) + == 607.11480 + ) + + with subtests.test("rated water out"): + assert ( + pytest.approx(prob.get_val("fuel_cell.rated_water_production", units="kg/h"), rel=1e-6) + == 683.61604984 + ) + + with subtests.test("rated co2 out"): + assert ( + pytest.approx(prob.get_val("fuel_cell.rated_co2_production", units="kg/h"), rel=1e-6) + == 417.512383 + ) + + with subtests.test("rated heat out"): + assert ( + pytest.approx(prob.get_val("fuel_cell.rated_heat_production", units="kW"), rel=1e-6) + == 844.389934 # TODO: update expected value + ) + + with subtests.test("rated stack efficiency"): + assert ( + pytest.approx(prob.get_val("fuel_cell.rated_stack_efficiency"), rel=1e-6) + == 0.63967037 # Example calculation + ) + + +@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.co2_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(258.4037579, rel=1e-2) + assert ng_consumed[6] == pytest.approx(1.0, rel=1e-2) + assert o2_consumed[6] == pytest.approx(84.18872, rel=1e-2) + assert water_out[6] == pytest.approx(94.7971630, rel=1e-2) + assert co2_out[6] == pytest.approx(57.8965188, 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) diff --git a/h2integrate/core/supported_models.py b/h2integrate/core/supported_models.py index a99961239..8d9beef62 100644 --- a/h2integrate/core/supported_models.py +++ b/h2integrate/core/supported_models.py @@ -85,9 +85,11 @@ def copy(self): "CustomElectrolyzerCostModel": "converters.hydrogen:CustomElectrolyzerCostModel", "WOMBATElectrolyzerModel": "converters.hydrogen:WOMBATElectrolyzerModel", "LinearH2FuelCellPerformanceModel": "converters.hydrogen:LinearH2FuelCellPerformanceModel", + "PEMH2FuelCellPerformanceModel": "converters.hydrogen:PEMH2FuelCellPerformanceModel", "H2FuelCellCostModel": "converters.hydrogen:H2FuelCellCostModel", "SteamMethaneReformerPerformanceModel": "converters.hydrogen:SteamMethaneReformerPerformanceModel", "SteamMethaneReformerCostModel": "converters.hydrogen:SteamMethaneReformerCostModel", + "SONGFuelCellPerformanceModel": "converters.natural_gas:SONGFuelCellPerformanceModel", "SimpleASUCostModel": "converters.nitrogen:SimpleASUCostModel", "SimpleASUPerformanceModel": "converters.nitrogen:SimpleASUPerformanceModel", "HOPPComponent": "converters.hopp:HOPPComponent",