System-Level and Stack Modeling

System-Level and Stack Modeling#

Multicell Stack Models and Assumptions#

A stack model aggregates electrochemical, thermal, and hydraulic behavior across \(N_c\) cells while preserving the dominant non-uniformities introduced by manifolds and manufacturing dispersion.

Stack-level electrical balance#

For series-connected cells, stack voltage is the sum of cell voltages:

\[ U_{\text{stack}} = \sum_{i=1}^{N_c} U_i \]

with each cell represented by:

\[ U_i = E_{\text{rev},i} + \eta_{\text{act},i} + \eta_{\Omega,i} + \eta_{\text{mt},i} \]

All cells share the same stack current \(I\) (series connection), while local current density differs if active area, contact resistance, or local conditions vary:

\[ j_i = \frac{I}{A_{\text{act},i}} \]

Species and mass balances#

At stack scale, hydrogen and oxygen production follow Faraday’s law:

\[ \dot n_{\mathrm{H_2}} = \eta_F\,\frac{N_c I}{2F}, \qquad \dot n_{\mathrm{O_2}} = \eta_F\,\frac{N_c I}{4F} \]

For manifolded channels, a common reduced form for each cell control volume is:

\[ \frac{d m_i}{dt} = \dot m_{i,\text{in}} - \dot m_{i,\text{out}} + \dot m_{i,\text{gen}} - \dot m_{i,\text{cross}} \]

where generation terms come from electrochemistry and crossover terms capture membrane permeation and leakage effects.

Stack thermal balance#

A lumped stack thermal model is often written as:

\[ C_{\text{th,stack}}\frac{dT_{\text{stack}}}{dt} = \dot Q_{\text{irr}} - \dot Q_{\text{cool}} - \dot Q_{\text{loss}} \]

with irreversible heat approximated by:

\[ \dot Q_{\text{irr}} \approx N_c I\left(U_{\text{cell,avg}} - E_{\text{th}}\right) \]

In higher-fidelity models, \(T_i\) is tracked per cell or per stack segment to capture thermal gradients along flow direction.

Manifolding assumptions#

To keep computation tractable, stack models typically adopt one of these manifold assumptions:

Uniform distribution assumption: equal inlet flow to each cell; suitable for control-oriented models.
Hydraulic network assumption: manifold and channel pressure drops computed explicitly; better for diagnosing maldistribution.
Quasi-steady hydraulics: pressure-flow algebraic equations coupled to slower electro-thermal dynamics.

Pressure-drop relations often use Darcy-Weisbach-type forms in manifolds and effective porous resistance in PTLs/channels.

Cell-to-cell variability modeling#

Real stacks exhibit dispersion in kinetics, resistance, and transport. A practical parameterization is:

\[ \theta_i = \bar\theta\,(1 + \delta_i), \qquad \delta_i \sim \mathcal{N}(0,\sigma_\theta^2) \]

for parameters such as \(j_{0,i}\), \(R_{\Omega,i}\), and \(j_{\text{lim},i}\). Spatially correlated variability (for example, edge vs center cells) is often more realistic than purely independent noise.

Typical observable consequences are:

Cell-voltage spread under load
Local hot spots and unequal aging rates
Different gas purity and crossover margins per cell

These effects motivate distributed sensing and balancing/control strategies (flow balancing, thermal balancing, current derating, and fault-tolerant operation).

Thermal and Fluidic Modeling#

Thermal and fluidic models describe how coolant, water, and gases are distributed through channels and manifolds, and how this distribution feeds back into electrochemical performance.

Coolant-side modeling#

At stack scale, coolant temperature is often modeled with 1D energy balances along the coolant path. For a control volume \(k\):

\[ \dot m_{\text{cool}} c_p \left(T_{k+1}-T_k\right) = \dot Q_{\text{gen},k} - \dot Q_{\text{loss},k} \]

where \(\dot Q_{\text{gen},k}\) is local heat generated by irreversible electrochemical losses and \(\dot Q_{\text{loss},k}\) is heat rejected to ambient or adjacent components.

In reduced-order dynamic form:

\[ C_{\text{th},k}\frac{dT_k}{dt} = \dot Q_{\text{gen},k} - h_k A_k\left(T_k - T_{\text{cool},k}\right) \]

This captures the key coupling between current density, temperature gradients, and local cell voltage.

Water feed and recirculation management#

Water management is primarily anode-side in PEM electrolysis. Models typically track:

Feed flow rate and stoichiometric excess
Local two-phase distribution (liquid water + oxygen bubbles)
Recirculation and separator return flow

Per-cell anode mass balance can be expressed as:

\[ \frac{dm_{\text{liq},i}}{dt} = \dot m_{\text{in},i} - \dot m_{\text{out},i} - \dot m_{\text{cons},i} + \dot m_{\text{return},i} \]

with consumption linked to electrolysis rate. Maintaining sufficient liquid inventory avoids local dry-out and excessive ohmic rise.

Gas management (H2 and O2 paths)#

Gas management models account for production, transport, separation, and pressure control:

O2 generation and disengagement on anode side
H2 production, transport, and drying on cathode side
Gas crossover and purity constraints

For each side, species balances combine generation and convective transport:

\[ \frac{dn_{s,i}}{dt} = \dot n_{s,\text{in},i} - \dot n_{s,\text{out},i} + \dot n_{s,\text{gen},i} - \dot n_{s,\text{cross},i} \]

where \(s\in\{\mathrm{H_2},\mathrm{O_2},\mathrm{H_2O}\}\).

Pressure drop in manifolds and channels#

Pressure non-uniformity drives flow maldistribution and therefore cell-to-cell variability. A common manifold/channel pressure model uses Darcy-Weisbach:

\[ \Delta p = f\,\frac{L}{D_h}\,\frac{\rho u^2}{2} \]

where \(f\) is friction factor, \(D_h\) hydraulic diameter, and \(u\) mean velocity.

For porous transport layers, pressure loss is often represented by Darcy’s law:

\[ \Delta p = \frac{\mu L}{K}\,u \]

with permeability \(K\). In two-phase regions, effective friction/permeability corrections are applied to account for gas-liquid interactions.

Coupled electro-thermal-fluid feedbacks#

The stack exhibits strong bidirectional coupling:

Higher current increases heat generation and gas void fraction.
Temperature rise changes membrane conductivity and kinetics.
Flow maldistribution alters local reactant availability and pressure.
Local pressure/phase changes modify crossover and transport losses.

Therefore, practical stack simulators solve electrical, thermal, and hydraulic balances together, either in fully coupled form or through staggered time-step updates.

Typical modeling assumptions for control-oriented models#

To reduce complexity while keeping predictive value, control-oriented models commonly assume:

Quasi-steady pressure field at each control step
Lumped coolant temperature per stack segment
Effective two-phase correction factors instead of full CFD
Aggregate separator dynamics for recirculation loops

These assumptions usually provide a good trade-off between computational cost and fidelity for control design, digital twins, and real-time estimation.

Dynamic Behavior and Transients#

Dynamic models are required when the stack is operated with variable renewable sources, where current commands, thermal inertia, and fluid delays interact on different time scales.

Load-change dynamics#

When current is ramped from \(I_1\) to \(I_2\), the cell voltage response is typically multi-time-scale:

Fast electrical response (milliseconds to seconds): activation and ohmic terms
Intermediate transport response (seconds): gas fraction, flow redistribution, pressure adaptation
Slow thermal response (tens of seconds to minutes): stack and coolant temperature evolution

A compact transient representation is:

\[ U_{\text{stack}}(t) = U_{\text{qs}}\big(I(t),T(t),p(t)\big) + \Delta U_{\text{dyn}}(t) \]

where \(U_{\text{qs}}\) is quasi-steady voltage and \(\Delta U_{\text{dyn}}\) captures delay effects (double-layer, transport, and manifold filling).

Ramp-rate limits are usually imposed to avoid excessive pressure/temperature gradients and local starvation.

Start-up sequence modeling#

Start-up is constrained by hydration, thermal conditioning, and gas safety. A simplified sequence is:

Water loop priming and purge
Temperature ramp to allowable operating window
Pressure stabilization and differential-pressure check
Current ramp from near-zero to target operating point

A state-machine model is often used in controls:

Idle -> Pre-fill -> Warm-up -> Pressurize -> Power-producing

Each transition is guarded by measured thresholds (flow, conductivity, pressure difference, temperature).

Shutdown and emergency transients#

Shutdown logic must prevent hazardous gas mixing and thermal shock. Typical controlled shutdown includes:

Current ramp-down with pressure coordination
Purge/vent sequence
Pump and valve post-run to remove residual heat

Emergency shutdown (trip) models are more abrupt and should account for residual heat and trapped gases after electrical disconnection.

Interactions with external power systems#

PEM stacks are usually connected through power converters, so dynamic behavior depends on converter control mode:

Current-controlled mode for hydrogen production tracking
Power-controlled mode for grid services or renewable-following

Grid-side disturbances (voltage sag, frequency events, converter limits) propagate to stack current and therefore to gas production and thermal states.

A simple coupling model is:

\[ P_{\text{el}}(t) = U_{\text{stack}}(t)\,I(t), \qquad \tau_I\,\frac{dI}{dt} + I = I_{\text{ref}}(t) \]

where \(\tau_I\) represents effective converter/current-loop dynamics.

Key transient indicators#

For model validation and control tuning, common transient KPIs include:

Voltage overshoot/undershoot during ramps
Time to reach new steady state after a load step
Maximum stack temperature and gradient during transients
Minimum local flow/pressure margins (starvation risk)
Temporary gas purity degradation during start/stop events

These metrics help define safe and efficient operating envelopes for renewable-coupled electrolyzer plants.

BoP models#

Balance-of-plant (BoP) models extend stack equations with auxiliary subsystems that dominate transient constraints and parasitic power. A control-oriented BoP model usually includes pumps, valves, pressure regulation, thermal loop dynamics, and gas-liquid separation.

Water pump model#

Water-feed and recirculation pumps are commonly represented by a first-order actuator plus a pump head-flow relation (Gulich, 2019; Karassik et al., 2008).

Actuator dynamics:

\[ \tau_p\,\frac{d\omega_p}{dt} + \omega_p = \omega_{p,\text{ref}} \]

Flow approximation (near operating point):

\[ \dot m_w \approx k_{p1}\,\omega_p - k_{p2}\,\Delta p \]

This linearized form is typically obtained from manufacturer pump curves around a nominal operating point and is well suited for control-oriented electrolyzer models (Gulich, 2019; ISO 9906:2012).

Hydraulic power and electrical consumption:

\[ P_{\text{pump}} = \frac{\Delta p\,\dot V}{\eta_p} \]

where \(\dot V\) is volumetric flow and \(\eta_p\) is pump efficiency. This term is needed to compute system efficiency, not only stack efficiency. BoP-oriented dynamic and efficiency analyses for water electrolysis systems also emphasize this auxiliary consumption term (Ursua et al., 2012).

Pressure regulation model#

Pressure is controlled using back-pressure valves, regulator valves, and optional compressor coordination. A simple manifold pressure state is:

\[ \frac{dp}{dt} = \frac{RT}{V}\left(\dot n_{\text{in}} - \dot n_{\text{out}} + \dot n_{\text{gen}}\right) \]

Valve flow can be represented with an orifice-type relation:

\[ \dot m_v = C_v\,u_v\,\sqrt{\max(p_{\text{up}}-p_{\text{down}},0)} \]

with valve command \(u_v\in[0,1]\). Differential-pressure constraints are enforced to protect membrane integrity:

\[ \left|p_{\text{anode}}-p_{\text{cathode}}\right| \le \Delta p_{\max} \]

Temperature regulation model#

Temperature regulation combines coolant flow control and heat-exchanger dynamics. A practical two-state form is:

\[ C_{\text{st}}\frac{dT_{\text{st}}}{dt} = \dot Q_{\text{irr}} - UA\left(T_{\text{st}}-T_{\text{cool}}\right) \]

\[ C_{\text{cool}}\frac{dT_{\text{cool}}}{dt} = \dot m_{\text{cool}}c_p\left(T_{\text{in}}-T_{\text{cool}}\right) + UA\left(T_{\text{st}}-T_{\text{cool}}\right) \]

Cooling pump and valve commands regulate \(\dot m_{\text{cool}}\) and effective \(UA\) to keep stack temperature within safe operating range.

Gas-liquid separator and recirculation model#

Separators are often modeled with holdup states and phase split efficiency:

\[ \frac{dm_{\ell}}{dt} = \dot m_{\ell,\text{in}} - \dot m_{\ell,\text{out}} - \dot m_{\ell\to g} \]

\[ \frac{dm_g}{dt} = \dot m_{g,\text{in}} - \dot m_{g,\text{out}} + \dot m_{\ell\to g} \]

Return flow from liquid holdup closes the water loop, while gas outlet quality influences downstream drying and purity metrics.

Compressor and downstream pressure interface (optional)#

When high-pressure hydrogen delivery is required, compressor dynamics can be approximated by:

\[ \tau_c\,\frac{d\dot m_c}{dt} + \dot m_c = \dot m_{c,\text{ref}} \]

with power consumption estimated from pressure ratio and isentropic efficiency. This subsystem can dominate plant-level energy consumption at high delivery pressure.

BoP control loops and coupling with stack model#

Typical coordinated loops are:

Current/power loop (converter)
Water-flow loop (pump speed)
Pressure loop (back-pressure valve)
Temperature loop (coolant flow or heat-exchanger bypass)

Because these loops are coupled through pressure, temperature, and two-phase behavior, model predictive or gain-scheduled control is often preferred over independent PI tuning at wide operating range.

BoP-oriented performance indicators#

To evaluate BoP model quality and control design, useful KPIs include:

Auxiliary power fraction: \(P_{\text{aux}}/P_{\text{el}}\)
Pressure-tracking error and differential-pressure violations
Temperature overshoot and settling time
Water-starvation margin and separator carryover risk
Net system efficiency (including BoP consumption)

Case Study: Stack Simulation Under Load Variation#

Present a representative transient simulation and key performance indicators.