Multiscale modelling of fuel cells and related energy conversion systems

: Heat and power cogeneration plants based on fuel cells are interesting systems for energy-conversion at low environmental impact. Different fuel cells have been proposed, but experimental testing rigs are expensive and the development of commercial systems is time consuming if based on fully experimental activities. Furthermore, tight control of operation of fuel cells is compulsory to avoid damage, which must be based on accurate models, able to predict the cell behaviour and prevent stresses and shut-down. Some selected examples of steady state, dynamic and fluid dynamic modelling of different types of fuel cells are proposed, mainly PEMFC and SOFC type. The general ideas behind the thermodynamic, kinetic and transport description are recalled, with some examples of models derived for single cells, stacks and integrated power co-generation units.

PEMFC preferable for mobility applications, whereas high temperature FCs are best fitted for stationary applications. Examples of application of PEMFCs in vehicles are available from different suppliers, e.g. Toyota Mirai, Nikola One and Chevrolet Colorado ZH2. The size of PEMFC system is up to 250 kW, less than 10 kW for single houses, ca. 50 kW for buildings.
As for stationary or large scale applications, FCs may be used as stand-alone devices or connected to the grid. This option is somehow more complex, since the FC supplies DC, that be converted to AC and conditioned to the voltage of the grid. Furthermore, sufficient stability of the cell should be ensured for grid connection. The development of simulation tools and models to predicts the transient and steady state behaviour of the cell is a precious contribution to predict the safe operating modes. Fault prediction is one of the major issues to improve system reliability, considering permanent (membrane deterioration, poisoning, leakage and cell aging), transient (flooding or drying) and external (cooling system, feed supply) faults. Furthermore, some cells such as PEMFC and SOFC are particularly sensitive to operating conditions. Modelling offers a method for improve the control, optimisation and reliability of these FCs.
In addition, models to describe the poisoning rate by CO and the consequent voltage losses and regeneration conditions are useful for performance monitoring and correction [2].
The phenomena underpinning FC operation and connection with an upstream fuel processor or a downstream user occur at very different time and length scales, requiring proper modelling tools to achieve the desired information. Careful validation against the basic kinetic, thermodynamic and electrochemical principles, as well as against experimental data is compulsory, but also the choice of the proper modelling methods is fundamental. Avery detailed review makes a comprehensive description of the available equations to compute transport phenomena, electrochemical thermodynamic and kinetics and cell performance [3]. The scheme of the interconnections between these modelling activities and the different time and length scales of each of them are sketched in Figures  1 and 2.   [3] under the Creative Commons Licence. 4 of 26 In this work recent modelling and simulation papers are reviewed to offer an up-todate panorama of the opportunities for multiscale modelling of FCs. Examples will be discussed from microscopic phenomena, such as the retrieval of transport parameters and kinetics, to the models to size and rate single cells, stacks and integrated power systems, including fuel processors and connection to the grid. The main features of the models are here presented and discussed, leaving the interested reader to the original reference for detailed description.

Kinetics, thermodynamics and transport phenomena
Computationally fast electrochemical models are useful to predict FC behaviour and provide the control logic. A 0-D thermodynamically consistent electrochemical model for PEMFC has been developed, including the transport of gaseous species along the channel and through gas diffusion layer (GDL) [4]. An optimal set of calibration parameters has been also retrieved on a physical consistent basis. Very good extrapolation capabilities for operation points outside the calibrated parameters space has been also demonstrated.
A mathematical model was built to evaluate the thermodynamic performance of a SOFC coupled with a thermoelectric generator and an absorption heat pump to recover the waste heat of SOFC exhaust gas [5].
A single-anode SOFC was tested to retrieve the characteristic curves to build a 2D steady-state simulation code in Fortran. This was based on a robust electrochemical kinetic model, with semi-empirical approach where the parameters derived from data fitting were used in physically derived equations [6]. To size and rate FCs, reliable kinetic models are needed, but they are multivariable device. Simulation involves nonlinear parameter estimation issues, often involving more than one local minima imposing the choice of robust optimisation algorithms [7,8].
Gas-liquid mass transfer in the porous layers or in the current collectors has been modelled in the case of a PEMFC [9]. Liquid water, vapor, and non-condensable gasses in the pores were modelled to assess the best operating conditions and criticisms.
Mathematical modelling of a composite continuum-network formulation has been developed to describe species diffusion and convection in gas diffusion layers in PEMFC. This helps the fluid dynamic study of the system developing a control volume mesh to monitor anisotropic effective transport properties as diffusivity and permeability [10].

Fuel cell and stacks
A residential power generator (60 kW) based on a PEMFC stack was designed to operate at 75°C, 100% humidity (regulated through careful water management) and 3 bar pressure [11].
The behaviour of a PEMFC stack was investigated considering stability and efficiency under different conditions, to optimise the operating points during connection to the grid [12]. The equivalent circuit of the simulated PEMFC is reported in Figure 3.  The detailed mathematical model is as follows [12]: Block A in the detailed scheme reported in Figure 3, calculates the conversion (U) of hydrogen and oxygen, based on the volumetric flow rates (V), composition (x%, y%), pressure of both gases (P) and temperature (T), assuming ideal gas behaviour.
The partial pressure of the different species and Nernst voltage are expressed as: Block B calculates V0C and IE and block C measures M (Tafel slope).  (charge transfer coefficient), G (activation barrier) and Kc (voltage constant) depend on the polarisation curve at the rated operating conditions and on the variation of the gas pressures at varying temperature.
Voltage and O2 utilisation vs. time are represented in Figure 4 [12]. Overall the scheme of the FC stack, including 45 cells, and the details of the block for connection to the grid are exemplified in Figure 5 [12]. The simulation showed that variable pressure of air did not affect significantly fuel usage, but strongly the efficiency. Critical pressures have been determined to ensure stable operation as 0.8 bar for air and 1.4 for the fuel. The latter parameter affected both fuel consumption and efficiency [12].
The performance of a PEMFC depends on many structural and working parameters (flowrates, pressure and relative humidity, membrane hydration, T and P of the cell, etc.), so that the development of robust models to optimise its operating optimum saves precious time during the set up of the device [8,13,14]. Nevertheless, the complexity of the governing equations and the multiple dependence from so many variables makes the modelling activity very hard.
As usual, models can be classified as empirical, semi-empirical and mechanistic, with a growing degree of complexity and accuracy/robustness of representation of the phenomena. In the present case purely empirical models are very weak, while the mechanistic ones are usually too complex and computationally expensive due to the high number of factors and the algebraic-differential nature of the problem [8]. Overall, semi-empirical models may be a reasonable compromise, with regression of some unknown quantities and inclusion/solution of the fundamental equations. Among these, the zero-dimensional steady state models, such as the Generalised Steady State Electrochemical Model (GSSEM), have been mostly employed due to their relative simplicity and versatility [15]. However, the need to retrieve some of the parameters brings about an additional problem, the selection of the optimization algorithm. Usually the characteristic I-V curve of the cell is fitted through the model to accomplish the non-linear regression of the missing parameters. Therefore, the typical objective function is:

= −
The calculated cell potential includes all the hypotheses and formulations of the model proposed and the unknown parameters to be regressed. The search algorithms applied to FCs are different, e.g. genetic algorithm [16], artificial neural network [17], artificial bee colony [18], harmony search [19], whale optimization algorithm [20] and some hybrid formulations. Often, the physical meaning of the retrieved parameters is not demonstrated when the optimisation search algorithm is insufficiently robust [8]. Typically, they fail in description of the electrochemical phenomena, evaluation of cell resistance and reversible voltage, computation of the concentration losses and estimation of the effect of operating parameters.
Some of the features and criticisms of GSSEM models are discussed in the following, full details can be found in the references [8,15].
The cell voltage is calculated as: Where ENernst is the thermodynamic potential, while the  terms represent the losses due to activation, ohmic and concentration effects, T(K) is the cell temperature. A scheme of the polarisation curve of a PEMFC is exemplified in Figure 6. i are parameters that can be derived in principle from thermodynamic and kinetic mechanistic equations and have a precise physical meaning.
with R the gas constant (8.3143 J/mol K), F the Faraday constant (96487 As/mol), ne the number of electrons per mole and c the charge transfer coefficient of the cathode, which can be either calculated from the values of 3 or 4, as regressed from the model. A consistent set of parameters would return similar values when calculating c from both parameters, which does not happen in most cases, as well examined in [8,13].
Where pH2 and pO2 are the effective partial pressures of the two species, Pi are the inlet pressures at the anode and cathode, Pw sat is the vapor pressures, RHi are the relative humidities of the anode and cathode side and i is the current density of the cell.
The voltage drops are differently computed. The activation term is due to the slow kinetics which causes overpotentials and lower current densities. The current density can be calculated through the Bulter-Volmer equation (i0 is the exchange current density, i.e. the rate of the reactions at equilibrium when the net current is nil, nF the charge transferred, i the charge transfer coefficient for both reactions and E/Er the actual and equilibrium potentials, kB and h the Boltzmann and Planck constants, Ci the surface concentration of the reacting compounds and G the Gibbs free energy): The current densities at the cathode and anode can be calculated as follows and consequently the activation losses, higher for the cathode than for the anode: The reversible potentials for hydrogen oxidation at the anode and oxygen reduction at the cathode are respectively 0V and 1.23V at T=25°C and P=1 atm.
The ohmic losses include the resistance towards electron and proton transport: indicating with A the active area and with Ri the resistances towards electron flow (usually unknown and to be retrieved by regression) and the membrane one, l is the thickness of the membrane and M is its specific resistivity, reported in the above example for a Nafion membrane, whose resistivity is strongly dependent on its water content .
The concentration voltage drop is instead related to mass transfer limitations, that may limit the feed of reactants to the electrodes surface, causing a voltage loss.
iL being the limiting current density and B another unknown parameter of the cell to be retrieved via regression. In the optimisation it is better to set the upper and lower limits for the search of the different parameters. For instance I = 0.2-2, CH+ =10 -11 -10 -4 , Gch,i = 10-70 kJ/mol, RC = 1-8 10 -4 ,  = 10-24 and B = 0.0136-0.5, as suggested in [8].
Once the voltage of a single cell has been modelled, the voltage of the stack can be calculated by multiplying the one of the single cell by the number of cells of the stack. The total power output is computed from the product of the current and the voltage of the stack [13].
Artificial Neural Network was used as robust search algorithm to extract the FC parameters in [21], a slime mould optimisation algorithm was applied in [7], a flower pollination algorithm was applied in [22].
A similar model was applied to interpret the experimental data collected over a 30 cells stack by varying temperature, pressure and relative humidity [23]. The model correctly represented the stack operation within ca. 10% error. A decrease of the activation losses was evident at increasing temperature due to the increase of the exchange current density. Also mass transfer limitations decreased with increasing temperature. An increase of pressure improved the operation due to higher partial pressures of both the reactants. The effect of relative humidity was more sensitive to the cathode than the anode due to the higher flowrate. Overall, the polarisation curves were correctly represented by the model, though with a systematic voltage overestimation by ca. 10% due to simplified model assumptions.
A model implemented in Matlab © has been developed for a PEMFC and for an integrated electrolyser/fuel cell, with economic analysis [24]. The model develops through mass balances at the anode, cathode and membrane, as synthetically described below. The 1-D model assumes uniform distribution of reactants and current. At the anode, H2 is the only species considered, considering its inlet flow and the one consumed, quantified by the Faraday equation: where n is the number of cells in the stack, I the current intensity and SH2 the stoichiometric ratio.
At the cathode, besides the oxygen balance, water must be also considered. where SO2 the stoichiometric ratio, RHc is the relative humidity at the cathode and Psat is the saturation pressure for water, calculated as a function of temperature. The total flowrate at the cathode is the sum of that of oxygen and water [24].
As for the membrane, no crossover of oxygen and hydrogen is considered, but only protons and water transport. They are modelled considering different mechanisms: electro-osmotic drag (eod), diffusion (diff), pressure gradients and thermo-osmosis, of which only the first two are considered significant. where I is the current intensity, A the active area, nd the eod coefficient, membr the water content of the membrane and DH2O the water diffusion coefficient, quantified as described in [24].
The scheme of the model is reported in Figure 7. Costs of the order of 6 10 4 US$ have been estimated, which is comparably high with respect to rival technologies.
The effect of gas crossover is of course important, almost unavoidable in PEMFC and attempts to consider it should be done. N2 and O2 flow to the anode and H2 to the cathode leading to parasitic reactions, mixed potentials and an overall reduction of efficiency. Furthermore it leads to degradation due to the formation of pinholes. An example of modelling gas crossover in HT-PEMFC is reported by Chippar and Ju [25]. A very accurate steady state model of a SOFC tube, based on experimental data, has been developed in [26], returning in general a 1% error (max deviation 8%). The model included mass and heat balances for both the solid and fluid ensuring a satisfactory account also of the thermal management of the system.

Failure modelling
P. Costamagna et al. widely studied in the last decade fuel cell systems, often based on SOFCs, developing specific models for the description of stacks and electrodes, integrated fuel processors and more recently their failure modes. In particular, two classical failure detection approaches, based on fault signature matrix or on data with statistical classification, were applied independently or combined. The random forest classification was used to select regular and faulty modes. A hybrid approach coupling a model-based scheme with a statistical classifier showed the most effective [27][28][29][30].
The degradation modes for SOFC have been extensively analysed by Eichhorn Colombo and Kharton [31], including coking, sulfur poisoning, formation of local hot spots due to excessive current densities, frequent power cycling. In order to model the probability of failure, or in other terms the reliability of the system, a hierarchy of analysis is proposed as in Figure 8 and the failure probability is correlated with the stresses of the cell or stack, which determines the correlation between current density and failure probability. The model can be integrated with a non linear equation solver to build the expected cell behaviour at a given time, building for instance the polarisation curve including the degradation of the cell (Figure 9). Different formulations can be given for the degradation rate, either derived from experimental data, or more easily visualised in terms of virtual age, which quantifies the extent of loss of each element as a function of current density, voltage or operating temperature [31]. Based on statistics, if a system is constituted of components connected in series, the system fails if any of the component fails. If they are connected in parallel, failure occurs when all the components fail. If a SOFC system is based on n independent identical stacks, k stacks being in operation occurs with a probability of: Where Pf is the failure probability of a stack. The reliability of the system is calculated as: Overall, at constant current density SOFC degradation rate increases with operating temperature, while at constant temperature increases with the current density. The failure probability can be plotted vs. the virtual age of the system, as in Figure 10.  [32]. Reproduced from [31] by kind permission of Elsevier.

Dynamic modelling
In addition to steady state modelling, dynamic simulation is important to highlight the transient behaviour of the system. Different models have been briefly reviewed in [9] and Matlab © and Simulink © tools are often used to this scope, as detailed in [13].
Besides the computation of all the variables described for steady state models, dynamic modelling adds a dynamic computation of the activation losses and of cell thermodynamics due to the formation of a charge double layer (CDL). The following time-dependent equations add: Where Rc sums the steady state losses for activation and concentration and C is the double layer charge [13]. Furthermore, the temperature of the cell does not remain constant and its value affects most of the variables of the model. Its variation with time is: Where CT is the thermal capacitance, H the global heat transfer coefficient, Tf the reference temperature and Tcell the lumped temperature of the FC [13].
For instance, a dynamic model of a Horizon H-500 PEMFC has been described [33], where the double layer formations has been described as a capacitor. A dynamic model of a PEMFC has been developed and validated against experimental data in [34] and in [35], focusing on the temperature profile through the cell, produced by the reaction and due to ohmic losses and removed through convection in the channels. Steady state assessment of voltage and ohmic losses was accompanied by dynamic modelling of concentration and activation losses in [36] and the dynamic response under step current variations in [37].
The integration of the heat exchange apparatus is also important, especially for FCs operating at high temperature. Heat can be used to sustain reforming to produce hydrogen, to feed a parallel turbine cycle based on steam or an Organic Rankine Cycle (ORC) or for residential heat cogeneration in the case of low enthalpy sources. An example of integration with an ORC is reported in [13] and with a water absorption cycle is described in [38]. In addition, the dynamic simulation of hot water supply simultaneously with power generation in a house has been reported in [39].

Fluid dynamic modelling
The detailed description of the flow field at the anode and cathode is an important piece of information. It affects the distribution of reactants (and thus cell performance), but also heat and water management. Accurate description is specifically needed at the cathode side, characterised by smaller gas diffusion coefficients. If flooding should occur, gas diffusion through the gas diffusion layer and at the catalyst layer would be inhibited, with overheating and dying of the membrane. Modelling the flow pattern implied the knowledge of the geometry of the channels, their discretization into a proper mesh and the solution of the governing equations of such problem through home developed algorithms, often aided by commercial or open packages, such as Ansys Fluent © or Open Foam  .
A button SOFC has been described with a 3D Computational Fluid Dynamic (CFD) model to retrieve the temperature, current and concentration patterns in a cross section of the cell, allowing better control of the operating conditions and preventing thermal stress and cell damage [40]. The model was based on electrode kinetics, mass transport, and elaborated in the ANSYS © Fluent environment, using a specific add-on for SOFC. An experimental investigation supported the activity, referring to a button SOFC with a La0.8Sr0.2MnO3 cathode and NiO-YSZ anode. The nodes of the model implemented were based on modelling fluid flow, including heat and mass transport into the channels and pores of the anode and cathode, current and potential in the solid material (including porous layers) and electrochemical reactions at the interface between the electrode and electrolyte [40]. The key steps of the model are the following.
The conservation of momentum is: In the equation  is the density, v the velocity, yi the composition in mass fraction and Si the sink or source of the i-th species.
The charge conservation is calculated as: with  as electric conductivity and  the Nernst potential calculated as: However, the actual cell potential is calculated taking into account the jump potential and the activation over-potential and the current density is calculated through the Butler-Volmer equation.
where act is the activation overpotential for the anode (a) or cathode (c). Equations for heat transfer computation are also available, such as: where Sh is the heat source and E the volumetric source. Details are reported in the original reference [40]. The model was implemented and loaded on a mesh of 3 million volume cells. The results showed increasing current density and power density for a given voltage and flowrates, and the parameters also increased with increasing flowrates ( Figure 11). The concentration profiles of hydrogen and oxygen at the electrodes are also shown. A very detailed 3-D fluid dynamic model has been developed by Atyabi et al. [41] to compare the flow pattern of oxygen between straight parallel and sinusoidal cathodic channels. CFD via finite volume and SIMPLE algorithm has been applied to model cathodic flowrate. The results showed that higher velocity was attained in the sinusoidal configuration, with increased diffusion towards the GDL, though at the expenses of higher pressure drop. Through the quantification of a uniformity index of different variables, the sinusoidal channels were characterised by more uniform oxygen mass fraction and temperature distribution, reducing the thermal stress [41].
In the same line, Kuo et al. [42,43] demonstrated that a wave flow field improves mass transfer and convection, with more uniform temperature profile. As well, bio-inspired flow patterns improved uniformity of reactants and consequently power density, with respect to parallel-spiral flow [44]. Counterflow configuration led to higher membrane conductivity than co-current flow [45].
A CFD framework has been used also to model diffusion and convection in porous transport layers, such as GDL. A macroscopic control volume was developed and implemented in a finite volume code with Ansys Fluent, embedding the microscopic porous network into the mesh [10].
Transient effects are significant for mobile applications and were the focus of a 3-D CFD study of a PEMFC with a serpentine channel. A step-change in the mass flow rates ends in current and power densities increase with time at low cell voltage due to concentration losses, but negligible at high voltages [46].

Water management
Water management is a fundamental isse to ensure stable performance for fuel cells operating at temperature lower than 100°C (e.g. in PEMFC), where liquid water is expected to form and accumulate in the cell [47][48][49]. On one hand the use of humid inlet gases is compulsory to avoid the drying of the membrane. The conductivity of the membrane is tightly related to its water content. If the cathodic gas nearby the membrane is not saturated with vapor, the membrane may release water, decreasing its proton conductivity and adding ohmic losses. The water transport in the membrane is also further complicated by the tranport of protons across it, which occurs as hydrated ions. The existence of a water concentration gradient from the cathode to the anode can stimulate the back diffusion of water. Based on the local temperature and pressure in the pores of the cathode, water can be either in vapor or liquid state. On the other hand the water forming during the reaction must be effectively discharged from the channels to avoid accumulation and thus flooding of the cell [50,51]. In such a case unstable performance may be observed, with abrupt voltage change when water is expelled. If both conditions coexist periodically, the expected lifetime of the cell is shortened [52]. Furthermore, the excess water prevents the correct transport and contact between the gases and the electrodes.
The former issue can be accomplshed by including an external humidifier, in form of a heat exchanger that provides some heat and vapor to the incoming gases. This is of course an additional item that adds costs and management issues to the assembly. On the contrary, a solution has been proposed to spatially resolve the heat removal from the cell so to achieve an appropriate humidity of the gases [52]. For instance, is the temperature of the cathode is allowed to increase only gradually prodiding accurate cooling, the gas can remain saturated along the channel, without the need of external humidity supply. A model has been developed to design a proper thermal profile, neglecting pressure drop along the channel, the thickness of the gas diffusion layers and hydraulic permeation acros the membrane. The model demonstrated the possibility to use the product water to keep the gases saturated through the channel, conclusion confirmed by experimental validation [52].
Detailed models have been also developed to investigate droplets formation and discharge in the channels of a PEMFC [51]. Multiphase flow field of oxygen, hydrogen and water must be studied to optimise operation and ensure the long term stability of the device. Microchannels are designed to narrow the depth and width of the flow field and modelled usually through finite element simulation. The pore size distribution of the carbon paper affects the droplets size and shape. Smaller channels increase the droplet pressure difference and improve their discarge rate [51].

Integrated fuel processors
A dynamic model based on simulation and experimental data was developed for a stand-alone PEMFC, including a reformer, a FC and the power conditioning system [53].
A similar system was also modelled for a residential power generation unit based on methanol reforming [54][55][56].
Hydrogen production via steam reforming of biomass-derived compounds, such as bioethanol is getting increasing attention [57][58][59][60][61][62][63], including the design of integrated processes for fuel processing and conversion into FCs. Ethanol in particular has the advantage of being a non-toxic liquid, being flexibly reformable and its material yield and energetic input can be tuned varying the water/ethanol ratio in the reacting mixture [64]. It can also be used at low purity levels [58,65], thus limiting the reactant cost. The catalysts for ethanol reforming are a substantially commercial reality [66,67] and can meet wide market options, such as residential power co-generation [68,69] or centralised production in biorefineries. An integrated plant for the joint electrical and thermal power production from diluted bioethanol, with a size of 8 -11 kW (of which Pelectical ≥ 4 kW) has been designed and rated, based on a PEMFC, with special emphasis on thermal integration [70][71][72]. This unit exploits diluted bioethanol as inexpensive feed with respect to the anhydrous pure ethanol. This point is the key for the economy of the generator and opens possibilities for thermal energy valorization.

Conclusions
Models for the description of different FC types have been derived over the years. They allow computation of the steady state or transient behaviour of the cell and can be used for the faster and inexpensive simulation of cell operation under different conditions. Simpler systems include only single cells or stacks, while more complex layouts include the fuel processor and, possibly, the connection with the grid. Failure models are also available to predict the reliability of the cells.
The thermodynamic description of the electrochemical system is available, together with kinetics for different FCs, mainly PEMFC and SOFC devices, which are the most investigated at the moment. Different overpotentials are also discussed, with the relative modelling tools. Finally, transport phenomena are accounted for through proper modelling, to estimate the efficiency of the system, but also to cope with the removal of water through the channels.
Overall, relatively easy models can be implemented with reasonably accurate predictions of the cell behaviour, allowing sizing and rating, but also optimisation in connection with fuel processors and lifetime prediction.
Funding: This research received no external funding.
Data Availability Statement: All the data are available in this paper and in the cited references.