Abstract

In a SOFC stack, manufacturing tolerances – cell warp, interconnect flatness, stack preload – set per-cell voltage and lifetime, but operators see only the stack terminal voltage, behind which the per-cell variability stays hidden. We resolve it with a 1D finite-volume cell solver, a phenomenological contact-resistance closure, and a multi-cell stack solver.

The model resolves per-cell behaviour that the stack terminal voltage hides. From an a-priori QC dispersion it projects a peak-to-min contact-overpotential spread of ca. 26 mV across a 30-cell stack at 0.5 A/cm² (\(10\)–\(52\) mV over the dispersion envelope), and a worst cell of median 38 mV rising to 76 mV at 1.0 A/cm² – a first-cell-to-fail margin invisible to terminal-voltage monitoring and resolved at real-time rates.

The cell electrochemistry is validated against published data – 8 of 12 experimental I–V cases within a 50 mV RMSE band from one parameter set, no curve fitting – while the stack-level spread is a projection from QC-anchored dispersion, with hardware validation pending. The solver runs at ca. 30 times real time per cell and 4.7 times for a coupled 30-cell stack inside a 300 ms hardware-in-the-loop budget. To our knowledge, no prior SOFC stack model couples first-principles cell physics, per-cell stack resolution, real-time HIL feasibility, and a-priori QC dispersion in one end-to-end chain.

Introduction

A digital twin of a solid oxide fuel cell (SOFC) stack must connect two domains that are normally modeled in isolation: the as-built mechanical state fixed at manufacture, and the electrochemical performance observed in operation. The bridging quantity is a per-cell ohmic-loss prediction grounded in the physics that governs both. No published model chain maps per-cell manufacturing dispersion – the spread captured by production quality control (QC) – to real-time operational voltage at hardware-in-the-loop (HIL) rates in a deployable form. The present work addresses this gap.

This paper reports an in-house multi-scale SOFC stack model. The architecture is deliberately hybrid: first-principles physics where it is available (cell electrochemistry, transport, reforming kinetics), and a calibrated phenomenological closure where it is not (the mechanical-to-electrical bridge that maps manufactured geometry to per-cell contact resistance). The two layers are kept separate so the provenance of each downstream prediction stays visible.

The contribution has three parts. First, a 1D finite-volume cell solver built on established physics: a Nernst potential from NASA Glenn 7-coefficient thermodynamics, symmetric Butler-Volmer activation (solved in inverse-sinh form), concentration-overpotential closures (limiting-current and Hafsi binary-diffusion [1]), and Xu-Froment internal reforming [2]. Second, the mechanical-bridge submodel: a geometric thresholding expression for contact area, coupled to a Holm constriction model, calibrated against simulation data, not a first-principles prediction. Third, a multi-cell stack solver that propagates a per-cell area-specific resistance (ASR) distribution – sampled from production-grade QC dispersion – to per-cell voltage at HIL-relevant rates.

Fig. 1 sketches the chain end to end. Validated on a steady-state I–V suite (8 of 12 experimental cases within a 50 mV band) and exercised on dynamic-transient and internal-reforming cases, the chain predicts a plate-dominated cell-to-cell spread of ca. 26 mV across a 30-cell stack – per-cell variation invisible at the terminal.

The same solver runs the validation suite and the assembly-twin scenarios reported here, and a service layer exposes it to an operator dashboard – a direct path from validation to deployment.

1 Scientific Approach

The assembly twin chain is: component QC data (interconnect (IC) plate flatness, cell warp, assembly tolerance) \(\to\) mechanical submodel \(\to\) per-cell contact-resistance map \(\to\) SOFC stack solver \(\to\) per-cell voltage \(V\), power \(P\), and fuel-utilization \(\eta_{\mathrm{FU}}\) spreads. The forward direction – predicting batch behavior from production tolerances – is the subject of this paper; the same mapping can be read backward to set tolerance budgets from end-use performance targets.

We bridge manufacturing and operation through the ohmic-loss chain because it is the chain most sensitive to mechanical variability. Cell warp, IC topography, and assembly-tolerance dispersion show up first as variation in real contact area between the cell electrode and the interconnect rib; contact area then maps linearly to contact resistance via the Holm constriction model [3]. Thermal coupling between cells is a secondary effect: once a cell sees an unusual current density due to its \(\ensuremath{R_{\text{contact}}}\), its operating temperature follows, and the stack-thermal coupling then redistributes heat.

2 Methods

This section documents the three layers of the model: the cell-level solver, the mechanical submodel, and the multi-cell stack solver. The validation suite is described last.

2.1 Cell-level solver

The 1D axial cell solver uses \(N_z\) control volumes spanning fuel and air inlets. For each control volume, seven gas-phase species (H\(_2\), H\(_2\)O, O\(_2\), N\(_2\), CO, CO\(_2\), CH\(_4\)) are tracked per channel; the balance equations are convection-reaction with electrochemical sources from Butler-Volmer kinetics and chemical sources from reforming.

\[\frac{\partial c_i}{\partial t} + \frac{\partial (u\, c_i)}{\partial z} = S_i^{\mathrm{ec}} + S_i^{\mathrm{ref}}, \label{eq:species}\]

with electrochemical source \(S_i^{\mathrm{ec}} = \pm j / (n_e F)\) on the relevant side of the electrolyte and aggregated reforming source \(S_i^{\mathrm{ref}}\). Charge transfer is carried by H\(_2\) only; CO contributes through the water-gas shift rather than direct electro-oxidation.

The solid energy balance is solved on a three-layer positive-electrolyte-negative (PEN) thermal structure (anode, electrolyte, cathode; the anode functional sublayer is resolved electrochemically but lumped thermally), with axial conduction, convective coupling to both channels, and convective/radiative losses at the stack-boundary cells. Its heat sources are ohmic (Joule) dissipation, the activation and concentration irreversibilities, the reforming reaction enthalpy, and the reversible entropic (\(T\Delta S\)) heat of the cell reaction.

The cell voltage is uniform along \(z\) and satisfies the local charge-balance constraint (Fig. 2):

\[\ensuremath{V_{\text{cell}}}= V_{\mathrm{Nernst}}(z) - \eta_{\mathrm{ohm}}(z) - \eta_{\mathrm{act}}(z) - \eta_{\mathrm{conc}}(z) - \ensuremath{R_{\text{contact}}}\cdot j(z), \label{eq:charge}\]

with each term closed as follows:

• \(V_{\mathrm{Nernst}}\): NASA-Glenn 7-coefficient thermochemical polynomials (McBride 1993; Burcat-Ruscic compilation), giving reaction-specific \(\Delta G^{\circ}(T)\) without ad-hoc temperature laws.

• \(\eta_{\mathrm{ohm}}\): electrolyte resistance with Arrhenius ionic conductivity \(\sigma(T)\), plus interconnect bulk resistance.

• \(\eta_{\mathrm{act}}\): symmetric Butler-Volmer (\(\alpha = 0.5\)) with separate anode and cathode exchange-current densities \(i_{0,a}(T)\), \(i_{0,c}(T)\), solved in inverse-sinh form to avoid catastrophic cancellation near the open-circuit voltage (OCV).

• \(\eta_{\mathrm{conc}}\): selected per case – a Hafsi binary-diffusion model [1] (Wang, Zhao, Zaccaria, Hafsi-calibration), a stagnant-film variant (Costamagna), or a limiting-current form (Aguiar); anode diffusivity from molecular and Knudsen contributions.

• \(\ensuremath{R_{\text{contact}}}\cdot j\): the per-cell contact term supplied by the mechanical bridge (Section 2.2).

Internal direct reforming is implemented via Langmuir-Hinshelwood-Hougen-Watson (LHHW) kinetics following Xu and Froment [2], with the three reactions of Table 1 and a per-reaction effectiveness factor (\(\eta_{\mathrm{eff}} = 5\times10^{-4}\), the Ni/YSZ-anode value of Wang et al. [4] for pore-diffusion limitation). The equilibrium constants \(K_{\mathrm{eq},i}(T)\) come from the NASA polynomial Gibbs free energies, consistent with the Nernst calculation, except water-gas shift, which uses the empirical correlation of Zaccaria [5].

Time integration is implicit-Euler by default (semi-implicit, A-stable), with adaptive \(\Delta t\) bounded by minimum \(0.01\) s and maximum \(600\) s.

2.2 Mechanical submodel (contact-resistance closure)

The mechanical-to-electrical bridge is where most stack-level uncertainty enters. The closure adopted here is the mechanical-submodel work of Klein [6]: a geometric thresholding expression for contact area, a closure-plane construction from preload and effective stiffness, and a Holm-bridge resistance assembly. It gives the per-node contact-area fraction:

\[P_{\mathrm{contact}}(c) = \frac{\varphi_{\mathrm{rib}}}{2} \cdot \mathrm{erfc}\!\left[\frac{c - \mu_z}{\sigma_z \sqrt{2}}\right], \qquad \sigma_z^2 = \sigma_{\mathrm{cell}}^2 + \sigma_{\mathrm{ic}}^2 + \sigma_{\mathrm{assembly}}^2, \label{eq:mech}\]

where \(\varphi_{\mathrm{rib}}\) is the rib coverage of the interconnect (land area over land+gap), \(c\) is the closure plane (a function of the applied preload and an effective stiffness), and the three variance contributions enter in quadrature (cell topography, interconnect topography, assembly tolerance). The per-node contact resistance then follows from a Holm constriction model [3] with temperature-dependent interconnect resistivity; \(\rho_{\mathrm{ic}}(T)\) follows a Crofer-22-APU experimentally calibrated curve.

The closure is calibrated against finite-element contact-mechanics simulations (Klein [6]), not against measured stack data; it is smooth and tractable for state-observer use but remains phenomenological, and its validation against measured hardware is pending.

2.3 Stack solver

A multi-cell stack solver couples \(N_c\) cells via the interconnect plate. Inter-cell coupling enters through two physical channels: heat conduction through the IC plate and electrical contact resistance at each rib interface. A third, statistical channel supplies a per-cell total-ASR value – sampled at the input-generation stage and passed to the solver as a per-cell contact profile – that lumps three end-of-line contributors: cell-intrinsic ASR scatter, interconnect-plate ASR scatter, and assembly contact-resistance scatter. Boundary cells see ambient losses; interior cells see only their two IC-plate neighbors. Each cell solves the full 1D axial physics on its own three-layer PEN. Inter-cell heat conduction through the IC plate is applied as a post-solve operator-splitting correction: it uses the freshly computed current-step temperatures (not a one-step lag) and is energy-conservative by construction; all stack runs use \(\Delta t \le 0.5\) s. The stack solver assembles per-cell terminal voltages and computes the stack terminal voltage as \(\sum_i V_{\mathrm{cell},i}\).

Per-cell ASR is sampled from a clipped normal distribution \(\mathcal{N}(\ensuremath{\mu_{\text{ASR}}}, \ensuremath{\sigma_{\text{ASR}}}^2)\) over the production QC range \([15, 120]\) m\(\Omega\cdot\)cm². Rather than fit \(\ensuremath{\sigma_{\text{ASR}}}\) to our own results, we anchor it to published cell-to-cell voltage statistics: Guida et al. [7] report a start-of-life voltage coefficient of variation of \(0.6\)–\(1\,\%\) across two 5-cell stacks, and removing the flow/thermal share [8] leaves a manufacturing/contact dispersion of order \(10\) m\(\Omega\cdot\)cm² (at \(\ensuremath{\mu_{\text{ASR}}}= 50\) m\(\Omega\cdot\)cm², \(j = 0.5\) A/cm²). We adopt \(\ensuremath{\sigma_{\text{ASR}}}= 13\) m\(\Omega\cdot\)cm² as an illustrative baseline and sweep the conservative \([5, 40]\) m\(\Omega\cdot\)cm² envelope (Section 3.4).

2.4 Validation suite and parameter provenance

The same parameter set per case is run across three regimes: a steady-state current-voltage (I–V) suite of twelve experimental cases plus one model cross-check, from six published groups (Zhao and Virkar [9], Wang [10], Zaccaria [5], Aguiar [11], Costamagna [12], [13], Hafsi [1]); a dynamic load-step transient (Zaccaria 2016 Part I); and an internal-reforming spatial-profile case (Aguiar 2004 DIR, \(81.6\,\%\) CH\(_4\) conversion against the paper’s ca. \(80\,\%\)).

Experimental cases are scored against digitized or paper-tabulated laboratory I–V data, with recorded digitization provenance and a ca. 8–30 mV pixel-error band. A model cross-check is scored against a paper’s own analytical model curve where no experimental I–V is published (Aguiar 2004, H\(_2\)-only); it is reported separately and never counted toward the experimental pass rate. Kinetic parameters are taken from the source tables where available (e.g. Costamagna and Honegger 1998 Table III) or solved in closed form from reported measurements (e.g. a two-point Arrhenius fit of \(i_0\) at 700 and 800 °C). One caveat is stated plainly: for cells where the exchange-current prefactor \(\gamma\) is back-fit to a single I–V point, a per-case calibration residual of ca. 25 mV is carried; no parameter is swept to minimize the full-curve RMSE.

3 Results and Discussion

3.1 Steady-state validation

Across the steady-state suite, 8 of 12 experimental I–V cases pass within a uniform 50 mV root-mean-square error (RMSE) band against digitized or paper-tabulated laboratory data (Fig. 3). Four cases are documented xfails, each with a named physical cause (below), and one further case (Aguiar 2004, H\(_2\)-only) is scored only against the paper’s own model curve – a model cross-check, not experimental validation, and excluded from the experimental pass count. The strongest predictive result is the Hafsi 2024 pair: both pass with the anode effective diffusivity derived from the paper’s own reported microstructure (pore radius 57 nm, tortuosity 1.5, porosity 0.2), not a fitted constant.

The four xfails trace to understood causes that lie outside the core physics rather than challenging it: a cell geometry a 1D through-plane model cannot represent (an integrated-planar cell [13]), under-prediction at high current density, faithful reproduction of a published parameterization that is itself optimistic, and cathode-kinetic extrapolation below the calibration window. The last – the most informative – is examined next (Section 3.2).

As a cross-condition check (excluded from the pass count), the Zaccaria parameterization run at the Jülich Mark-F 18-cell operating point [14] – a different temperature and fuel – reproduces the body-of-stack cell voltage (\(0.84\)–\(0.85\) V) to within 18 mV; with the cathode parameterization inherited, this is a consistency check, not an independent validation.

3.2 Extrapolation limit: Zhao multi-temperature

Zhao and Virkar [9] provide the suite’s multi-temperature stress test. The lanthanum strontium manganite (LSM) cathode kinetics (\(\gamma_c = 1.64\times 10^{14}\) A/m\(^2\), \(E_{\mathrm{act},c} = 214\) kJ/mol) are anchored on the apparent \(i_0\) ratio between \(800\) and \(700~°C\), then applied across cells and temperatures. The model passes at \(700~°C\) and for the optimized cell at \(800~°C\) (Fig. 4). The two misses are not core-physics errors: at \(800~°C\) the standard cell tracks the paper’s own optimistic Table-5 Butler-Volmer formula rather than its measured Fig. 5 cell, and the \(600~°C\) case is a kinetic extrapolation limit.

The \(600~°C\) failure mechanism is the single-Arrhenius form, not a too-narrow fit window: the kinetics \(i_0(T) = \gamma_c \exp(-E_{\mathrm{act},c}/RT)\) are already anchored at both \(800~°C\) and \(700~°C\), but cannot capture the curvature of \(\log i_0\) vs \(1/T\) observed in LSM cathodes below ca. \(700~°C\) (multi-step oxygen reduction, oxygen-vacancy ordering).

3.3 Internal reforming: Aguiar 2004 DIR axial profile

The Aguiar 2004 direct-internal-reforming (DIR) SOFC case (S/C \(= 2\), 10 % pre-reformed CH\(_4\)) gives a total outlet CH\(_4\) conversion of 81.6 % against the paper’s ca. \(80\,\%\), and reaches 50 % conversion at \(z/L = 0.45\) (Fig. 5), within Aguiar’s reported \(z/L \approx 0.4\)–\(0.5\). Both the integrated conversion and the half-conversion location are reproduced even though the model uses Xu-Froment LHHW steam-methane-reforming (SMR) kinetics (\(E_{\mathrm{act}} = 240\) kJ/mol, with CO/H\(_2\)O denominator inhibition and a per-reaction effectiveness factor) rather than Aguiar’s Achenbach form (\(E_{\mathrm{act}} = 82\) kJ/mol, no adsorption denominator).

3.4 Assembly-twin scenario: 30-cell batch sensitivity

This subsection runs the full chain end to end – the paper’s contribution: turning a production-grade QC dispersion into a per-cell voltage spread across an assembled stack. The statistical results (Figs. 6a, 7, 8) use the linear ohmic mapping \(\ensuremath{\eta_{\text{contact}}}^{(i)} = \mathrm{ASR}^{(i)} \cdot j\) per Monte-Carlo realization, which follows from Eq. [eq:charge] when the non-contact overpotentials are common-mode across cells: absolute kinetic uncertainty then shifts the stack baseline but cancels in the per-cell spread, which depends only on the ASR dispersion and \(j\). The coupled stack solver (Sections 2.3 and 3.5) handles the dynamic load-step (Fig. 6b) and the HIL benchmarks.

Under this mapping the cell-to-cell variability is plate-dominated: body-of-stack \(\eta_{\text{contact}}\) at \(0.5\) A/cm² is ca. 25 mV, the highest-ASR outlier ca. 13 mV above it (Fig. 6a). A \(0.7\) A/cm² load-step through the coupled solver shows the per-cell spread tracking the load (Fig. 6b). The median peak-to-min \(\eta_{\text{contact}}\) at the data-anchored baseline (\(\ensuremath{\sigma_{\text{ASR}}}= 13\) m\(\Omega\cdot\)cm², derived in Section 2.3) is ca. 26 mV – the best-to-worst spread, distinct from the worst-cell absolute level treated below (Fig. 8). A spread of this size can produce a first-cell-to-fail signature over long-duration operation, while remaining below what the stack terminal voltage alone reveals.

Fig. 7 confirms the dominance ranking is invariant across the \([5,40]\) m\(\Omega\cdot\)cm² envelope, not just at the anchor: operating current density \(j\) is the largest lever, total-ASR dispersion \(\ensuremath{\sigma_{\text{ASR}}}\) second, plate mean \(\ensuremath{\mu_{\text{ASR}}}\) smallest (at the lowest dispersion \(j\) and \(\ensuremath{\sigma_{\text{ASR}}}\) become comparable). The baseline spread (median of the \(\ensuremath{\sigma_{\text{ASR}}}= 13\) panel) is \(26.1\) mV, bootstrap 90 % CI under \(\pm 0.5\) mV (\(N_{\mathrm{MC}} = 4000 \times 400\) samples).

Fig. 8 reports the worst-cell \(\eta_{\text{contact}}\) distribution at the anchor over four current densities. At \(0.5\) A/cm² the worst-cell median is 38 mV (95th percentile 44 mV); at \(1.0\) A/cm² it climbs to 76 mV (95th percentile 88 mV), approaching the individual-cell alarm regime while the stack terminal voltage may still be marginally acceptable. The worst-cell distribution is the production-relevant signal for first-cell-to-fail probability over the operating envelope.

Across the \([5, 40]\) m\(\Omega\cdot\)cm² sweep the spread median grows from 10 to 52 mV (26 mV at the \(\ensuremath{\sigma_{\text{ASR}}}= 13\) anchor), with bootstrap 90 % CIs \(\le \pm 0.5\) mV (\(N_{\mathrm{MC}} = 4000\)). The result is robust to distribution shape: lognormal and Weibull samplers matched on \((\ensuremath{\mu_{\text{ASR}}}, \ensuremath{\sigma_{\text{ASR}}})\) give worst-cell medians of 40 and 37 mV against 38 mV for the clipped normal.

3.5 Real-time performance

The 1D single cell runs at ca. 30\(\times\) real time at \(0.5\) A/cm², with a 99th-percentile step time of 139 ms – inside the 300 ms hardware-in-the-loop ceiling (only rare boundary-condition phase changes briefly exceed it). At stack scale the solver uses a representative-cell architecture: three cells (top, interior, bottom) carry the full 1D physics and are interpolated to the \(N_c\) interconnect-plate positions, so per-step cost is independent of \(N_c\). A coupled 30-cell transient then runs at \(4.7\times\) real time, while a full-cell variant (every cell solved independently) is two orders of magnitude above budget. The per-cell contact overpotential is \(\ensuremath{R_{\text{contact}}}\cdot j\) by definition (Eq. [eq:charge]), so the batch statistics in Section 3.4 apply this linear mapping directly; a uniform-contact 30-cell solve bounds the residual non-contact spread below ca. 1 mV, confirming cell-to-cell thermal coupling as a second-order effect (Section 1).

What per-cell resolution buys a stack builder. The worst-cell \(\eta_{\text{contact}}\) distribution (Fig. 8) is the production-relevant output of the chain. Before assembly it supports four decisions: predicting the batch voltage spread from measured production tolerances, flagging the worst cell early, setting tolerance budgets backward from end-use performance targets, and comparing cell or plate suppliers on a common operational footing.

Because the per-cell total-ASR dispersion at end-of-assembly is not published directly – that composite figure sits in proprietary stack-builder QC databases – the anchor is taken from published cell-to-cell voltage statistics and bracketed by the \([5, 40]\) m\(\Omega\cdot\)cm² sweep; a measured IC-plate QC time-series can replace the sampled distribution without changing the downstream chain.

4 Conclusions

We have presented a multi-scale SOFC stack model that maps as-built per-cell manufacturing tolerances to per-cell operational voltage at HIL-relevant rates. On the steady-state suite it passes 8 of 12 experimental I–V cases within a 50 mV RMSE band, with four documented xfails (each with a named physical cause) and one model cross-check; the Hafsi 2024 cases pass with the anode diffusivity taken from the reported microstructure, no fitted transport parameter. In the assembly-twin application, with \(\ensuremath{\sigma_{\text{ASR}}}\) anchored to published voltage statistics, the predicted peak-to-min \(\eta_{\text{contact}}\) spread across a 30-cell stack is ca. 26 mV – invisible at the stack terminal – with operating current density and total-ASR dispersion the dominant levers. The single-cell solver runs at ca. 30\(\times\) real time inside a 300 ms HIL budget, and the representative-cell stack solver holds that budget at \(4.7\times\). Future work: a parallelized full-cell coupled solver at HIL rates, an extended Kalman filter observer, and live IC-plate QC data replacing the sampled distribution.

Acknowledgements

The work has been conducted in the framework of the AMPS (Automated Mass Production of SOC Stack) project (www.amps-project.eu). The project is supported by the Clean Hydrogen Partnership and its members (GA: 101111882). The mechanical-submodel framework is F. Klein’s MSc work [6]. The graphical abstract was drafted with an AI image-generation tool (OpenAI DALL\(\cdot\)E) and edited by the authors; all technical content was specified and verified by the authors. Solver, inputs, and Monte-Carlo scripts will be deposited in the EFCF Zenodo community with a DOI on publication.

Keywords: EFCF 2026, SOFC, Digital Twin, Assembly Twin, Manufacturing Variability.