# 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:

1. Water loop priming and purge
2. Temperature ramp to allowable operating window
3. Pressure stabilization and differential-pressure check
4. 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.*