Source-linked AI summary
Free open reference implementation of a two-phase PEM fuel cell model
Roman Vetter, Jürgen O. Schumacher
TL;DR
PEM fuel-cell modeling lacks transparent, accessible reference implementations, making common development and comparison difficult. This paper provides a compact open MATLAB implementation of a steady-state, non-isothermal two-phase five-layer MEA model, with simulations showing strong moisture dependence and material-property sensitivity. The resulting tool supports rapid testing, comparison, validation, and extension within its one-dimensional scope.
Problem
Transparent and easy-to-understand PEM fuel-cell implementations are scarce, forcing repeated redevelopment and limiting accessible model comparison.
Method
The paper implements a steady-state, non-isothermal two-phase macro-homogeneous five-layer MEA model as a compact standalone MATLAB program focused on through-plane transport.
Results
The simulations show performance is strongly moisture-limited, with dry/hot scenario T2 producing the worst prediction and water-transport parameterizations materially affecting predicted performance.
Takeaways & Limitations
The open implementation enables rapid material-parameter studies, model comparison, validation, and extensions using an accessible reference code.
Takeaways & Limitations
The model uses a simplified Dirichlet liquid-water boundary condition and a one-dimensional setup that makes some stoichiometry tests inapplicable.
Abstract
from arXiv · showhide
In almost 30 years of PEM fuel cell modeling, countless numerical models have been developed in science and industrial applications, almost none of which have been fully disclosed to the public. There is a large need for standardization and establishing a common ground not only in experimental characterization of fuel cells, but also in the development of simulation codes, to prevent each research group from having to start anew from scratch. Here, we publish the first open standalone implementation of a full-blown, steady-state, non-isothermal two-phase model for low-temperature PEM fuel cells. It is based on macro-homogeneous modeling approaches and implements the most essential through-plane transport processes in a five-layer MEA. The focus is on code simplicity and compactness with only a few hundred lines of clearly readable code, providing a starting point for more complex model development. The model is implemented as a standalone MATLAB function, based on MATLAB's standard boundary value problem solver. The default simulation setup reflects wide-spread commercially available MEA materials. Operating conditions recommended for automotive applications by the European Commission are used to establish new fuel cell simulation base data, making our program a valuable candidate for model comparison, validation and benchmarking.
1 Introduction
The paper addresses the lack of transparent, accessible PEM fuel-cell implementations by presenting a lightweight open reference model for five-layer MEAs. It targets through-plane transport modeling with a compact implementation suitable for extension and comparison.
- The paper targets the shortage of fully transparent, easy-to-understand PEM fuel-cell implementations, which has forced repeated redevelopment and slowed modeling adoption.
- It presents a free, standalone macro-homogeneous five-layer MEA model for low-temperature PEM fuel cells.
- The model represents three-dimensional physics through a one-dimensional formulation focused on dominating through-plane transport processes.
- The implementation accounts for non-isothermal and two-phase behavior while using a compact, readable codebase.
- The five-layer geometry divides physical through-plane transport processes among distinct subdomains.
2 Mathematical model
The model represents a five-layer MEA with coupled electrochemical, charge, heat, gas-species, vapor, ionomer-water, and liquid-water transport processes. It combines conservation laws and constitutive transport relations with reaction, phase-transition, and heat-source terms, while noting unresolved challenges in liquid-water interface treatment.
- Geometry: The one-dimensional MEA comprises a PEM between two catalyst layers and two gas diffusion layers, with gas channels and plates represented as boundaries.Microporous layers are not explicitly modeled but can be added.
- Governing equations: Eight coupled steady-state second-order PDEs conserve charge, heat, gas species, dissolved water, and liquid water through the MEA.Nonlinearity arises because transport coefficients and source terms depend on the model’s dependent variables.
- Transport processes: Charge and energy transport use Ohm’s law for electronic and protonic fluxes and Fourier’s law for heat conduction, while gas species follow Fick’s law under uniform pressure.The model neglects thermal diffusion and uses the ideal gas law for interstitial gas concentration.
- Water transport: Water transport combines ionomer back diffusion and electro-osmotic drag with unsaturated-flow liquid transport driven by saturation gradients.Liquid water saturation denotes the fraction of pore space filled with liquid water.
- Sources and reactions: Electrochemical reactions and phase transitions couple the equations through source terms for reactant consumption, water production, sorption, evaporation, condensation, and associated heat effects.Faraday’s law determines hydrogen and oxygen consumption, while latent, Joule, ionic, activation, and Peltier contributions enter the heat balance.
- Boundary conditions: Liquid-water boundary conditions remain delicate because dynamic droplet formation and detachment are difficult to represent with steady-state area-averaged conditions.The text identifies physically accurate liquid-water interface formulations as an ongoing research topic.
3 Parameterization
The model uses literature-based parameterizations for commercial MEA materials and for electrochemical, ionomer, water, and porous-media transport properties. These choices provide a reference setup for PEMFC comparison and benchmarking.
- Material selection: Nafion NR-211 membranes and Toray TGP-H-060 carbon-paper GDLs parameterize a commercially representative MEA for reference simulations.Standard literature data supply the remaining material-independent electrochemical and physical properties.
- Electrochemical parameters: The cathode uses temperature- and concentration-dependent ORR exchange-current parameterization, while HOR and ORR are represented with Butler–Volmer kinetics.The ORR is treated as the rate-limiting half-reaction, and the cathode platinum loading is set three times higher than the anode.
- Ionomer-related parameters: Nafion ionomer properties include a percolation-based, Arrhenius-corrected conductivity law and a refitted water diffusivity parameterization based on experimental data.The diffusivity fit is plotted at 80 °C and corrected for partial ionomer content in catalyst layers using a Bruggeman factor.
- Water and porous-media transport: Water and gas transport use literature correlations for saturation pressure, electro-osmotic drag, sorption, vapor mass transfer, gas diffusion, capillary pressure, and hydraulic permeability.Partially flooded porous-media gas diffusivities include porosity, tortuosity, and a (1 − s)^3 saturation correction.
- Water and porous-media transport: The CGDL/GC liquid-water boundary uses a capillary estimate of sC ≈ 0.12 derived from pc ≈ 2 kPa for Toray TGP-H-060.The pressure boundary does not constrain the interfacial liquid-water flux, which adjusts to satisfy the imposed saturation condition.
4 Numerical implementation
The model is distributed as a standalone MATLAB function using built-in numerical tools and a free BSD-licensed source release. Moderate solver tolerances provide rapid voltage sweeps with a quantified discretization error.
- Implementation: MMM1D solves the coupled model with MATLAB’s standard bvp4c boundary-value solver and requires only built-in functionality.The solver uses finite differences, a 3-stage Lobatto IIIa implicit Runge–Kutta method, and automated residual-based mesh selection.
- Implementation: The complete source code is freely available under the 3-clause BSD license for commercial and noncommercial use with attribution.The standalone implementation is intended to maximize compatibility and portability.
- Numerical performance: At relative and absolute tolerances of 10^-4 and 10^-6, the reference implementation averages 54 mesh nodes and completes a full 50 mV voltage sweep in a few seconds.The reported maximum absolute (relative) discretization error is 4.4 mA/cm^2 (0.23%).
5 Simulation results
The model produces through-plane polarization, potential, flux, and base-case summary outputs, then applies EU-harmonized JRC stress tests to establish comparison data. Results show strong moisture dependence, with temperature, humidification, and pressure effects shaping predicted performance.
- Base case: The model automatically generates polarization curves, through-plane potential profiles, and flux profiles across the MEA for the base case.The base-case outputs include a polarization curve over cell voltages and profiles evaluated at multiple voltages.
- Base case: Membrane proton-potential curvature reflects spatially varying proton conductivity caused by declining ionomer water content toward the anode.Strong electro-osmotic drag toward the cathode drives the water-content decline; higher diffusivity, lower drag, or higher low-water conductivity would flatten profiles and increase current density.
- Base case: The base-case output reports key figures including ohmic membrane resistance, average MEA temperature, and mean ionomer water content, evaluated at U = 0.6 V.The integration covers the catalyst-coated membrane and can be repeated at other operating points.
- EU harmonized stress tests: The JRC stress tests vary temperature, gas humidification, and pressure to provide normalized data for PEMFC model comparison, validation, and benchmarking.Tests T1–T7 are applicable to the 1D model; T8 and T9 are excluded because they vary stoichiometry.
- EU harmonized stress tests: Stress-test performance is strongly moisture-limited: humid T1 performs best at moderate current densities, while dry/hot T2 performs worst because its membrane is extremely dry.At low cell voltages, hotter but drier T3–T5 can outperform T1 because the 45 °C temperature in T1 limits ionic conductivity; T4 slightly exceeds T3 through stronger dissolved-water back diffusion.
- EU harmonized stress tests: Dry T6 and T7 show little pressure sensitivity because the differential fuel cell does not reach the high-current regime where fuel or oxidant starvation becomes important.The pressure response changes with thicker, less permeable diffusion media and more humid gases.
6 Conclusion
The paper provides a compact, open MATLAB implementation of a two-phase PEMFC model designed for accessible simulation, modification, comparison, and validation. Its usefulness is bounded by deliberately simplified physics and substantial uncertainty in transport parameters, while the framework remains extensible.
- The open MATLAB implementation enables readily accessible PEMFC simulation, parameter substitution, feature addition, and parameter studies.It contains less than 400 lines of commented code, requires no third-party software or compilation, and permits commercial and noncommercial use.
- The model intentionally omits numerous physical and chemical effects, including convection, pressure-driven permeation, degradation, unsteady phenomena, and ice formation.These omissions may matter under certain operating conditions.
- Literature variation in model parameters and transport coefficients significantly affects results, with liquid-water flux identified as a particularly unresolved source of uncertainty.The authors note that satisfactory universal solutions are not yet available for these modeling aspects.
- The model can run a single simulation in about one second on an ordinary laptop, supporting time-critical applications and rapid parameterization tests.The authors also describe simplifications that can reduce resource requirements further.
- The framework is designed for extension through added subdomains, alternative transport equations, additional liquid-water treatment, and refined parameterizations.Examples include MPLs, Brinkman or Maxwell–Stefan transport, Nernst–Planck equations, and temperature-dependent material properties.
Nomenclature
The nomenclature defines variables for thermodynamics, geometry, temperature, voltage, gas composition, and transport properties used throughout the PEMFC model. It includes quantities for both anode and cathode gas channels and for liquid-water behavior.
- The nomenclature includes thermodynamic quantities such as temperature, pressure, Gibbs free energy, enthalpy, and molar volumes.
- Transport and electrochemical nomenclature covers gas diffusivity, dissolved-water diffusion, reaction symmetry, porosity, overpotential, permeability, and surface tension.
- It defines geometric and state variables including through-plane coordinate x, cell voltage U, mean temperature T, and reference temperature Tref.
- Gas-composition variables distinguish anode water vapor and hydrogen from cathode water vapor and oxygen, alongside saturation vapor and dry-gas fractions.