Add new PEM and NG solid oxide fuel cell models by genevievestarke · Pull Request #794 · NatLabRockies/H2Integrate

genevievestarke · 2026-06-24T20:29:02Z

Add new PEM and NG fuel cell models

This PR adds two new fuel cell models to H2I:

A PEM fuel cell model
A natural gas solid oxide fuel cell model

Each of these models include an I-V curve and basic reaction modeling to determine the feedstocks consumed.

Section 1: Type of Contribution

Section 2: Draft PR Checklist

Open draft PR
Describe the feature that will be added
Fill out TODO list steps
Describe requested feedback from reviewers on draft PR
Complete Section 7: New Model Checklist (if applicable)

TODO:

Implement reaction efficiencies?

Type of Reviewer Feedback Requested (on Draft PR)

Structural feedback:
There is some potential overlap here with the steam natural gas reforming model that already exists in the hydrogen folder. In the SO NG fuel cell model, the fuel cell itself still uses a reaction of hydrogen and oxygen, but the natural gas is put through a steam reforming process before the deconstructed H_2 and CO are fed into the fuel cell. The addition of other things in the fuel cell input other than pure hydrogen does affect the fuel cell performance, so I think it's a worthwhile addition to the code on its own. In the model added, natural gas goes into the model, and electricity, water, and carbon dioxide come out.
Let me know your thoughts on possible overlap issues!

Implementation feedback:
Are these models where you would expect them to be in the code?
Thoughts on Solid Oxide Natural Gas model (SONG)? Should it be NGSO (natural gas solid oxide) instead?

Other feedback:
Are these being worked on/needed in/for other projects?

Section 3: General PR Checklist

Section 4: Related Issues

Section 5: Impacted Areas of the Software

Section 5.1: New Files

path/to/file.extension
- method1: What and why something was changed in one sentence or less.

Section 5.2: Modified Files

path/to/file.extension
- method1: What and why something was changed in one sentence or less.

Section 6: Additional Supporting Information

Section 7: Test Results, if applicable

Section 8 (Optional): New Model Checklist

elenya-grant

hi Gen! I really only looked at the PEM_h2_fuel_cell.py method but wanted to give some preliminary feedback. My main thoughts are:

I think keeping track of what's calculated at the cell level vs stack level vs the system level is a bit confusing in the compute() method and the calc_current function. I think some extra in-line comments would be helpful for improving clarity.
I think that there are some mis-matches on units (kg/hr vs kg/dt)
I think the looping in the compute() method could be removed.
using I for current and J for current density would help me follow the logic!

Overall - this is a great start! I haven't dug into the doc page or the examples or tests yet but would love to on a re-review! Let me know if you want to chat about any of my comments sometime! Happy to hop on a call!

elenya-grant · 2026-07-01T19:22:45Z

+        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


some of these already exist in h2integrate/tools/constants.py. I think you could remove these as attributes and import the constants instead. Some of the constants that exist in constants.py are:

constants.py var = your variable equivalent

f_c = faraday

H_MW = M_H2*1e3/2

O2_MW = M_O2*1e3

H_MW*2 + (O2_MW/2) = M_H2O/1e3

HHV_H2_MJ_PER_KG = hhv_h2/1e6

I think the other constants you have here could be moved to constants.py also!

elenya-grant · 2026-07-01T19:53:21Z

+    # fuel_cell_efficiency_hhv: float = field(validator=range_val(0, 1))
+
+
+def calc_current(power_ref, cell_area, n_cells, stack_number):


I may not be fully appreciating what this function does, but I don't think the power_curve part is needed, power = current*voltage. I have some ideas on possible revisions but add those in later.

I think what's partially confusing about this method is what scale it should be on. like power_ref is the power to the entire system. Then power is divided among stacks. Every cell in the stack has the same current and voltage. below is what I tried to do instead and it gave similar results as your method:

def calc_current_refactor(power_to_system, cell_area, n_cells, n_stacks): # note: power_to_system should be in W J_curve = np.array([ 0.0356, 0.05413333, 0.0796, 0.11366667, 0.244, 0.454, ]) # A/cm2 voltage_curve = np.array([ 0.987, 0.936, 0.884, 0.838, 0.786, 0.736]) # 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) J_coefs = np.polyfit(P_curve, J_curve, 5) J_P_curve = np.poly1d(J_coefs) power_per_stack = power_to_system/n_stacks stack_current_density = J_P_curve(power_per_stack) stack_current = stack_current_density*cell_area cell_voltage = V_J_curve(stack_current_density) return stack_current, cell_voltage, V_J_curve

elenya-grant · 2026-07-01T19:56:15Z

+    I_cell = max(power_I_curve(power_density), 0)
+    V_cell = V_I_curve(I_cell)
+    I_cell = I_cell * cell_area


I may need to work out the units for myself to check the logic here but I think that if you're doing this:

I_cell = I_cell * cell_area

that makes me think that current_curve is actually current density (A/cm2), even though the comment says A. could you clarify whether the I_cell being returned is current or current density and similar question for current_curve?

You're right, the current curve is for current density! I'll clarify that!

elenya-grant · 2026-07-01T19:58:59Z

+            "oxygen_consumed",
+            val=0.0,
+            shape=self.n_timesteps,
+            units="kg/h",


I think the units are f"kg/({self.dt}*s)" for the hydrogen consumed, oxygen consumed, and the water output.

elenya-grant · 2026-07-01T20:21:26Z

+            )
+
+            # Calculate hydrogen and oxygen consumed
+            H2_consumed_rate = ((I_cell * self.N_series * self.M_H2) / (2.0 * self.f_c)) * (


I think it may be nice to split this line up a bit:

h2_consumption_rate_stack =(n_cells * I_cell) / (2 * faradays) # mol/s h2_consumption_rate_system = n_stacks*h2_consumption_rate_stack*H2_MW/1e3 # kg/s h2_consumed = h2_consumption_rate_system*self.dt

elenya-grant · 2026-07-01T20:27:37Z

+        h2o_generated = np.zeros(self.n_timesteps)
+        commodity_out = np.zeros(self.n_timesteps)
+
+        for i in range(self.n_timesteps):


I think this loop could be removed?

here's a potentially alternative approach to the looping

power_demanded = np.clip(inputs[f"{self.commodity}_command_value"], a_min=0.0, a_max=inputs["system_capacity"]) i_stack, v_cell, V_J_curve = calc_current_refactor(power_demanded*1e3, cell_active_area, n_cells, self.config.n_stacks) h2_in_kg_per_s = inputs["hydrogen_in"]/3600 o2_in_kg_per_s = inputs["oxygen_in"]/3600 # stack current i_from_h2 = (h2_in_kg_per_s*2*f_c)/(M_H2*self.config.n_stacks*n_cells) i_from_o2 = (o2_in_kg_per_s*4*f_c)/(M_O2*self.config.n_stacks*n_cells) i_op = np.minimum.reduce([i_stack, i_from_h2, i_from_o2]) j_op = i_op/cell_active_area v_op = V_J_curve(j_op) power_out = v_op*n_cells*i_op*self.config.n_stacks/1e3 # kW H2_consumed_rate = ((i_op * M_H2) / (2.0 * f_c)) * ( 3600 * self.config.n_stacks * n_cells ) # kg/hr O2_consumed_rate = ((i_op * M_O2) / (4.0 * f_c)) * ( 3600 * self.config.n_stacks * n_cells ) # kg/hr

note that calc_current_refactor is a function that I put in an above comment.

elenya-grant · 2026-07-01T20:30:04Z

+        # 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)


should cell active area be recalculated so that stack_size = n_cells*cell_active_area*power_dens?

elenya-grant · 2026-07-01T20:52:05Z

+            # print(self.stack_size, self.n_cells)
+
+            # TODO:
+            if H2_consumed_rate > H2in or O2_consumed_rate > O2in:


H2in and H2_consumed_rate won't be in the same units for non-hourly timesteps. H2_consumed_rate is in kg/dt but H2in is in kg/h.

Co-authored-by: elenya-grant <116225007+elenya-grant@users.noreply.github.com>

elenya-grant

Howdy Gen! This is looking great! I had a few small comments but nothing too serious. I looked primarily at the SO_NG_fuel_cell - but some of the comments are relevant to both models. I'm curious if theres any way to utilize some base-class here, but it's tough since they're in different folders. Anyway - I think it's looking good and would be happy to chat through any comments if ya want! thanks!

I just want to double-check some units-stuff (it never hurts to have multiple eyes on units conversions) but once I do that - I'll approve this! I imagine that it just needs a few small updates, some commented out code removed, then it's good to go!

elenya-grant · 2026-08-10T17:23:51Z

+        )
+
+        self.add_output(
+            "carbon_dioxide_out",


some of the carbon capture technologies (in converters/co2/marine/) use co2 in the variable names instead of carbon_dioxide. I think that it'd be nice for us to use the same naming convention for the same commodities. Anyone have a preference on carbon_dioxide vs co2? (@kbrunik)

elenya-grant · 2026-08-10T17:24:40Z

+        M_CO2 = CO2_MW / 1000  # Molar mass of CO2 in kg/mol
+
+        # calculate max input and output
+        inputs["oxygen_in"]  # kg/h


reminder to remove this line

elenya-grant · 2026-08-10T17:28:02Z

+        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)


Could you add a comment noting what the 8 and 4 here represent? I'm assuming its number of moles or electrons transferred?

elenya-grant · 2026-08-10T17:29:02Z

+
+        # 7. Calculate hydrogen and oxygen consumed and water produced
+        #       based on electrochemical reactions
+        ng_consumed = ((I_cell * M_CH4) / (8.0 * faraday)) * (


maybe 8 and 4 should be set to variables? Like n_ng = 8 and n_o2 = 4 and n_h2o=2? just so that its easier to see where they're used?

elenya-grant · 2026-08-10T17:31:34Z

+        self.add_input(
+            "system_capacity",
+            val=self.config.system_capacity_kw,
+            units="kW",
+            desc="Capacity of the solid oxide natural gas fuel cell system",
+        )


what if this input was rated_electricity_production? as the output from the performance model? I don't think thats strictly necessary but the rated_electricity_production output follows some standard naming convention whereas the system_capacity input is semi-unique to this performance model (aka - we haven't standardized input design variable naming convention for converters so it could be nice to use the output names)

elenya-grant · 2026-08-10T17:36:39Z

+        outputs["water_out"] = h2o_generated
+        outputs["carbon_dioxide_out"] = co2_generated


the fuel cell primary output commodity is electricity, but also produces water and carbon dioxide. If someone wanted to (for whatever reason) calculate the LCO of the carbon dioxide or the water production from the fuel cell, the fuel cell model would also have to output rated_carbon_dioxide_production and rated_water_production (those outputs may be nice to add just in case)

elenya-grant · 2026-08-10T17:39:05Z

+        self.add_input(
+            "system_capacity",
+            val=self.config.system_capacity_kw,
+            units="kW",
+            desc="Capacity of the h2 fuel cell system",
+        )


similar comment here as the SO NG fuel cell cost model about changing this input to rated_electricity_production instead?

Genevieve Starke added 19 commits May 20, 2026 14:10

Initial fuel cell model with psuedo code

4551f87

Initial fuel cell model

9115bbe

fuel cell example - draft

9529cf6

Update example

040cf3b

updates

7e2c693

Remove print statements

0352910

Merge in develop

f7e8b90

add new model to hydrogen init file

6bb5751

Updates to model

bee8478

Merge branch 'develop' into feature/add_new_fc

4bff858

Add more detailed FC cost model, update for passthrough controller

30205b3

Update interpolation model and add NG FC model

8b47b45

Update doc strings

e10050b

Remove print statements from slc

fcb5c96

Update naming to align with naming conventions

1aa8d93

remove leftover example input

737fe35

Add documentation and changelog

f1ebbd4

Add tests for PEM fuel cell

0442a86

Add SONG tests and model updates

63ac542

elenya-grant reviewed Jul 1, 2026

View reviewed changes

Genevieve Starke and others added 2 commits July 6, 2026 16:50

Update PEM regression tests

7e6d362

Update h2integrate/converters/hydrogen/PEM_h2_fuel_cell.py

a1f97ea

Co-authored-by: elenya-grant <116225007+elenya-grant@users.noreply.github.com>

elenya-grant added the needs modifications This PR has been reviewed, at least partially, and is ready for PR author response label Jul 8, 2026

Genevieve Starke added 2 commits July 10, 2026 19:08

Add PEM FC updates from other branch

af7944a

Update SO NG FC

59001b2

elenya-grant reviewed Aug 10, 2026

View reviewed changes

		# fuel_cell_efficiency_hhv: float = field(validator=range_val(0, 1))


		def calc_current(power_ref, cell_area, n_cells, stack_number):

		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)

		outputs["water_out"] = h2o_generated
		outputs["carbon_dioxide_out"] = co2_generated

Uh oh!

Conversation

genevievestarke commented Jun 24, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Add new PEM and NG fuel cell models

Section 1: Type of Contribution

Section 2: Draft PR Checklist

TODO:

Type of Reviewer Feedback Requested (on Draft PR)

Section 3: General PR Checklist

Section 4: Related Issues

Section 5: Impacted Areas of the Software

Section 5.1: New Files

Section 5.2: Modified Files

Section 6: Additional Supporting Information

Section 7: Test Results, if applicable

Section 8 (Optional): New Model Checklist

Uh oh!

elenya-grant left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

elenya-grant Jul 1, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

elenya-grant Jul 1, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

elenya-grant left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

genevievestarke commented Jun 24, 2026 •

edited

Loading

elenya-grant Jul 1, 2026 •

edited

Loading

elenya-grant Jul 1, 2026 •

edited

Loading