An Innovative Hybrid 3d Analytic-numerical Approach for System Level Modelling of Pem Fuel Cells

The PEM fuel cell model presented in this paper is based on modelling species transport and coupling electrochemical reactions to species transport in an innovative way. Species transport is modelled by obtaining a 2D analytic solution for species concentration distribution in the plane perpendicular to the gas-flow and coupling consecutive 2D solutions by means of a 1D numerical gas-flow model. The 2D solution is devised on a jigsaw puzzle of multiple coupled domains which enables the modelling of parallel straight channel fuel cells with realistic geometries. Electrochemical and other nonlinear phenomena are coupled to the species transport by a routine that uses derivative approximation with prediction-iteration. A hybrid 3D analytic-numerical fuel cell model of a laboratory test fuel cell is presented and evaluated against a professional 3D computational fluid dynamic (CFD) simulation tool. This comparative evaluation shows very good agreement between results of the presented model and those of the CFD simulation. Furthermore, high accuracy results are achieved at computational times short enough to be suitable for system level simulations. This computational efficiency is owed to the semi-analytic nature of its species transport modelling and to the efficient computational coupling of electrochemical kinetics and species transport. Symbols Channel cross-sectional area Velocity potential in channel domain concentration in channel domain water activity arg Argument of a nonlinear function Butler-Volmer function ℬ Velocity potential in GDL#1 domain concentration in GDL#1 domain , Water vapour concentration Velocity potential in GDL#2 domain ℂ concentration in GDL#2 domain • • • Membrane sulphonic group concentration Concentration of dissociated protons in membrane Total molar concentration of gas catalyst surface concentration of reaction product catalyst surface concentration of reactant reference concentration concentration of saturated water vapour Gaseous binary diffusion constant diffusion constant of protons in membrane diffusion constant of water in membrane , , coefficient of proportionality of proton diffusivity to water content , , coefficient of proportionality of water diffusivity to water content Faraday constant Linearised nonlinear function Nonlinear function Either or ℂ ℎ Sum of channel height and GDL thickness ℎ1 GDL thickness ℎ2 Channel height Electric/proton current density Cathode catalyst effective exchange current density , Molar flux of water Molar flux of protons Length of representative unit Depth of slice Index of modes along coordinate Coordinate normal to boundary plane/index of modes along coordinate Molar flux Energies 2013, 6 5428 Electroosmotic drag coefficient Pressure Gas constant/membrane ohmic resistance Source or sink term for species …


Introduction
Proton exchange membrane fuel cells (PEMFC) are emerging as an advantageous power source in the alternative power sources segment due to their high power density and zero tank-to-wheel emissions. However, as analysed in [1], commercialization of fuel cell technology is lagging behind owing to several disadvantageous characteristics of the fuel cell, summarised by [2] as: cost and complexity, immaturity, and its role as replacement technology. These facts call for increased application of advanced mathematical modelling and simulation tools to efficiently address and tackle the shortcomings of today's fuel cell technology [2]. PEMFC are characterized by fully interrelated transports of mass, charge and heat, which are governed by convection, diffusion, phase transition and electrochemical reactions. Many fuel cell models can be found in the literature or are available as packages in commercial software as systematically analysed in [3][4][5]. Among these models the most comprehensive are the ones based on 3D computational fluid dynamic (CFD) modelling with extended capabilities of solving electrochemical reactions (e.g., [6][7][8][9]). CFD-based models indeed provide insight into detailed phenomena along with high level of predictability and high accuracy; however, they are also characterized by very long computational times inherent to their comprehensive modelling framework. Long computational times are a major disadvantage whenever a large number of simulation results are required to be obtained in reasonable amount of time which is typical for tasks such as: (a) rapid testing of different design options [10,11] including the ultra large scale [9]; (b) modelling of larger systems powered by fuel cells, e.g., the entire fuel cell powered vehicle [12,13]; (c) modelling transient operation of the fuel cell powered systems [14,15] and (d) development of control functionalities [15,16]. Among these, fast computation is of most crucial importance in models for which real time capability is mandatory such as in [12,14,16] or any hardware-in-the-loop model application. Models featuring shorter computational times can be categorised as: -Fully empirical models such as equivalent circuits, neural-network models, support and relevance vector machines which are all highly accurate and require very short computational times. However, these, unlike the CFD models, lack any predictability since they are characteristic of a specific (i.e., non general) system part and are thus useful in simulations only as computer model representation of an existing or a non-modifiable component. Examples found in references: [17][18][19][20].
-Reduced dimensionality models which are predictable mechanistic models that feature reduced computational load yet retain higher accuracy of current density predictions. These are typically 2D models (e.g., [28], the 2D model in [29]) modelling straight channel fuel cell configurations (such as depicted in Figure 1) solving governing equations in the plane that is parallel to the gas flow and perpendicular to the membrane (e.g., green plane in Figure 1). Further simplification of 2D models leads to the so called 1D + 1D (e.g., [30]) models where mathematical treatment of the dimension that runs along direction of gas flow addresses only the most dominant physical phenomena (e.g., mass transport via bulk gas flow in channel). Applying these 2D models to a fuel cell geometry such as presented in Figure 1 leads to certain systematic discrepancies since the area through which mass is exchanged between the gas diffusion layer (GDL) and the channel is considerably smaller than the area where mass is exchanged between the membrane and the GDL. To compensate for this discrepancy some models employ additional correction parameters, obtained through fitting results of full 3D CFD models or experimental data, yielding pseudo 3D models (e.g., [31,32] and Sherwood number adjusted 2D model in [29]). An alternative to this pseudo 3D modelling is the so called 2D + 1D approach (as found in [29]) where the physical phenomena are fully addressed in the plane perpendicular to the gas flow and the reduced treatment of the 1D + 1D approach is used in the direction of gas flow.
These pseudo 3D and 2D + 1D models aim to take into account all three dimensions whilst being considerably faster than 3D CFD simulations. However with the exception of [31], they are still much slower than 1D models due to their inherent numerical nature. The model found in [31] features computational times comparable to 1D models owing to its analytic nature. However the approximate analytic solution for direct liquid fuel cells in [31] is only valid when the velocity along channels can be approximated as constant.
An article by Tavčar and Katrašnik [33] proposes a hybrid 3D analytic numerical approach to modelling species transport (HAN-ST) in a PEM fuel cell as a way to achieve full 3D resolution whilst featuring short computational times. The principles of HAN-ST modelling in [32] are presented on a theoretical straight channel co-flow fuel cell of very simple geometry and its core principle can be summarized as: The blue and the spotted translucent regions represent the membrane and GDL respectively, the green surfaces represent the plane of symmetry across the middle of a rib and the yellow surfaces the plane of symmetry between two ribs. (a) A five-rib parallel channel co-flow fuel cell geometry. The two symmetry planes that apply to all ribs are shown on one half of rib defining the representative unit; (b) Representative unit with indicated slice (red); (c) Slice as a sliced-out section of the representative unit that has sufficiently small depth to be treated as a 2D object and (d) Slice domains from top to bottom: cathode channel domain, cathode GDL domains (left cathode GDL#1, right cathode GDL#2), MEA domain, anode GDL domains (left anode GDL#1, right anode GDL#2), anode channel domain. A quasi 3D model constructed by taking a 1D numerical model for the gas-flow along the channel and superimposing onto it a 2D analytic solution for the plane perpendicular to the gas-flow. This proves to be computationally efficient due to the additional computational load introduced by the calculation of 2D analytic solution being of the same order of magnitude as the computational load of the base 1D calculation. It was shown in [33] that for the simple theoretical geometry of fuel cell (FC), HAN-ST enables resolving a 3D species transport with very high accuracy while preserving computational speed characteristic for 1D models owing to its partially analytic nature. Due to its high accuracy and computational efficiency the HAN-ST model offers a promising basis for an efficient standalone PEM fuel cell model suitable for system level simulations. HAN-ST as presented in article [33] is a species transport sub-model and needs to be coupled to an electrochemistry sub-model in order to yield a full standalone fuel cell model. Additionally the simple theoretical geometry of HAN-ST from [33] enables devising the 2D analytic solutions on a single domain which is practically impossible in cases of more complex realistic geometries.
As a full realistic standalone PEMFC model derived from the HAN-ST of [33], this paper presents and describes an innovative isothermal standalone hybrid 3D analytic numerical fuel cell model (HAN-FC) featuring: 1. The HAN-ST principle extended to more general realistic straight parallel channel FC geometries by devising the 2D analytic solution on a jigsaw puzzle of multiple domains and 2. An electrochemistry sub-model efficiently coupled to this extended HAN-ST sub-model. Thus, in terms of categorization introduced above, this is a 2D + 1D model resembling that found in [28], with the difference that it features analytically (instead of numerically) resolved species concentration in the 2D cross-sectional plane in both GDL and channel, making it computationally much faster. The nature of HAN-FC's analytic solution also fundamentally differs from the one in [31] where the 2D analytic solution is obtained for the plane that is parallel to the gas flow and perpendicular to the membrane. Differences of HAN-FC from other similar models found in literature are further discussed in Section 3.
The model assumptions and geometry of HAN-FC are outlined in Section 2 and the underlying governing equations of the species transport and electrochemical reactions are given in Section 3. For validation the simulation results of the HAN-FC is benchmarked with the results obtained by a professional 3D CFD tool: the simulation setup of HAN-FC and CFD is presented in Section 4, and the comparative results and their analysis is provided in Section 5. In all instances the accuracy of HAN-FC is to be understood as the accuracy relative to the CFD results, which is typically done for the 2D and pseudo 3D models [28,29,31,33]. Section 6 summarizes the main conclusions. Detailed mathematical derivation of the model is given in Appendix.

Model Assumptions and Geometry
The modelled fuel cell is a straight, co-flow, hydrogen-oxygen type PEM fuel cell (FC) with its geometry taken from a laboratory test fuel cell. The fuel cell topology and geometry are illustrated in Figure 1 where the following sub-elements can be distinguished: cathode channels, cathode-side gas diffusion layer (GDL), thin catalyst layer for oxygen reduction, hydrated proton exchange membrane, thin catalyst layer for hydrogen oxidation, anode-side GDL, anode channels. Due to the symmetric geometry of the fuel cell a half of one rib, as depicted in Figure 1, is taken as the representative unit of the fuel cell. The derivation of HAN-FC is based on three sets of assumptions: (1) regarding gaseous species transport, (2) regarding species transport in membrane and (3) regarding electrochemical reactions. The gaseous species transport assumptions are summarised as follows:

I.
A steady state solution of the problem is sought. II.
The problem is isothermal meaning a constant uniform temperature is assumed for the whole fuel cell and no energy equation is calculated. III.
The gases are treated as ideal.

IV.
There is no liquid water in GDLs and Channels. V.
Pressure variations are small enough that constant pressure is assumed in gas equations and the gas flow is assumed incompressible over the whole fuel cell. This is justifiable due to very small pressure variations in the fuel cell owing to the short straight channel geometry. VI.
The diffusion system is always bi-componential (either oxygen and water vapour or hydrogen and water vapour) with a constant binary diffusion coefficient. VII.
The effective (macroscopic) diffusion constant used in GDL is the diffusion constant for empty space (which applies in the channel) divided by the tortuosity of GDL. VIII.
Diffusion in gas in the direction of channel gas-flow is neglected due to the large aspect ratio of the relevant fuel cell dimensions.

IX.
In GDL there is no convective transport in the direction of channel gas flow. The very small pressure drop from inlet to outlet produces only negligible convective flow in the GDL in the direction along channel length. X.
The flow in channels is laminar.
HAN-FC extends the assumptions regarding membrane species transport made by HAN-ST amounting to the following: XI. There is no gas crossover in the membrane. Thus water is the only species transported through the membrane besides protons. XII. The mobility of both water molecules and protons (in form of hydronium ions) at a given point in membrane are linearly proportional to the membrane water content at that point. XIII. Species transport within the membrane only occurs in the direction perpendicular to the membrane sheet. XIV. The average number of water molecules dragged by the electro-osmotic pull of each proton is an independent constant. XV. The concentration of mobile protons is constant within the whole membrane. This is due to the humidification of membrane being always sufficiently high that full or close to full proton dissociation can be assumed. XVI. The membrane water content at the membrane/gas boundary is in equilibrium with water vapour concentration in the gas.
The assumptions regarding electrochemical reactions inherent to HAN-FC are: XVII. Only the electrochemical kinetics of oxygen reduction on the cathode catalyst is considered. The electrochemical kinetics of hydrogen oxidation on the anode catalyst is much faster and its contribution to the over-potential is several orders of magnitude smaller and thus neglected. XVIII. The catalyst layer is infinitely thin and characterised by an effective exchange current density as current per unit area [A/m 2 ]. XIX. The electrical resistance of the GDL is neglected and uniform constant electrical potential is assumed within the whole GDL.
All the above assumptions follow the assumptions made by other models of comparable modelling depth presented in the introduction.

Governing Equations
The modelled fuel cell is made of number of equal parallel symmetrical ribs where a half of one such rib is indicated between the green and the yellow plane in Figure 1a. The symmetry condition between ribs is defined in the GDL and membrane (yellow surfaces in Figure 1) requiring that there be no species flux across the symmetry plane. In the case of the two outermost ribs the wall boundary of the fuel cell also allows no species flux across that boundary. Thus, taking into account that zero velocity in the direction of channel gas flow is assumed in GDL (assumption IX), this wall boundary condition is the same as the symmetry condition allowing treating all ribs as identical. Since the two halves of one rib are mirror images and since all ribs are identical the so defined half of one rib is a representative unit of the fuel cell. This representative unit is sliced along the direction of gas flow into a number of thin slices as depicted in Figures 1b,c. Within each slice only variation of variables in the plane perpendicular to channel gas flow is addressed effectively making a slice a 2D object. The slice has three parts: the cathode gas part, the membrane electrode assembly (MEA) part and the anode gas part. The three parts are, as discernible form Figure 1d, further divided onto seven computational domains: each gas part is divided into one channel domain and two GDL domains (GDL#1, GDL#2 as depicted in Figure 2) whilst the whole MEA part is contained within its MEA domain. The so defined computational domains, depicted in Figure 2 as shallow rectangular cuboids, share the 2D nature of the slice and are thus simple rectangles. An analytic 2D solution for the distribution of species concentration in a slice is sought as a jigsaw puzzle of seven individual 2D solutions for each domain that are appropriately coupled to one another. The "Source" arrow represents the convective inflow of species from the channel domain of the previous upstream slice and the "Sink" arrow represents the outflow into the channel domain of the next downstream slice. Blue arrows represent the molar fluxes of water traversing the borders between domains that serve as boundary conditions for the individual solutions in domains. The ribbed purple surface represents the catalyst layer, i.e., the MEA/GDL interface. At the grey surfaces, i.e., the walls, and at the green and yellow surfaces, i.e., the two symmetry planes, the boundary condition requires zero flux through the boundary.
The analytic solutions for the consecutive slices are, due to assumptions VIII, IX and XIII, coupled only via the two bulk gas flows in channels. The cross-sectional velocity profile of the gas flow (i.e., profile in the 2D plane of a slice) is assumed to take on the dome shape of laminar flow (graphically shown in Section 5.2). This is similar to the 1D laminar velocity profile used in [34] and [28] however, unlike [28] where this profile does not change along the channel, in the HAN-FC model the profile scales with average cross-sectional gas velocity that varies along the channel. This variation of average gas velocity along the channel is calculated stepwise from slice to slice, i.e., as a numerical 1D pipe flow model. The net convective transport into and out of the 2D plane of a slice is then obtained by multiplying the 2D velocity profile by the 2D species concentration distribution profile. This 2D profile of convective transport thus differs from the uniform plug flow used in [29]. HAN-FC essentially treats the bulk gas flow as a 1D flow with a superimposed simple 2D velocity profile of laminar flow and a 2D species concentration profile which enables coupling of the 2D solutions of consecutive slices. The two 1D gas flows are calculated numerically meaning that the change in mean gas velocity is calculated stepwise from slice to slice.
Obtaining 2D analytical solutions for slices that are coupled to one another via the perpendicular numerically resolved 1D gas flow gives this approach the name "hybrid 3D analytic-numerical". The HAN-FC model thus gives full 3D information on species concentration distribution, differing from the approach of [31] where the model is essentially 2D and the effects in the third dimension are accounted for by means of calibration parameters.
The division of the gas part onto three computational domains is inherent to HAN-FC and differs from HAN-ST of [33] where the very simple geometry allows each gas part to be contained within a single computational domain. Sections 3.1 and 3.2 describe in detail the equations governing production, consumption and transport of species in the different types of calculation domains, Subsection 3.3 describes the coupling conditions between domains of one slice and Subsection 3.4 describes the procedure for obtaining a solution for the whole fuel cell by finding consecutive solutions for all the slices.

Gas Part
In the gas part only species transport takes place. The species transport in an ideal gas is modelled by Stefan-Maxwell diffusion equations which, for the case of binary gas mixtures (i.e., gas mixtures of two components) as assumed in VI, read [35]: where vectors and are the molar fluxes of the two components: the water vapour and the reactant (either oxygen for cathode gas part or hydrogen for anode gas part) respectively, and are the corresponding molar concentrations, is the net molar velocity of the gas, is the binary diffusion coefficient for the binary solution of the two components and = , , is the 3D nabla operator. Due to assumptions II and V, the temperature and pressure are constant and thus the total molar concentration = + = is also constant which means that and are complementary and the distribution of both concentrations can be obtained by solving only for and obtaining the other by: Taking the latter into account and adding together Equations (1) and (2) the expression for net molar velocity is obtained: All the subsequent derivation of governing equations for species transport in gas will be aimed at finding the distribution of only water vapour concentrations by means of solving Equations (1) and (4) since the reactant concentration is directly given by Equation (3) and thus the subscript shall be omitted.
The requirement for a steady state solution is expressed as: where is the water vapour molar concentration and the corresponding molar flux. The assumption of incompressible flow (V) requires that divergence of velocity be zero: At this point it is worth introducing a distinction that splits the general 3D notation into a 2D + 1D notation reflecting the different treatment of physical phenomena in the dimension along the channel gas-flow and in the dimensions perpendicular to the channel gas-flow. This distinction will be used in deriving the "2D + 1D type" equations in the following two subsections. Placing the gas part into the coordinate system defined in Figure 2 a 2D nabla operator∇, a 2D Laplace operator ∇ , a 2D velocity and a velocity component along the direction of channel gas flow are defined in the cross-sectional plane as: (the 3D nabla is not to be confused with its 2D version∇). The condition in Equation (6) can be thus rewritten as: Rewriting Equation (5) using the 2D + 1D notation and taking into account the assumption of no diffusive transport along the coordinate (VIII), leading to neglecting the term , reads: Equation (9) is solved for each computational domain. In favour of finding an analytic solution for species transport in the plane within a computational domain the term ∇ ⋅ ( ) in equation (9) is approximated as: where ̃ is the average concentration in the relevant computational domain. , being a 2D vector field, can be represented as a sum of gradient of a scalar field ( , ) and "2D curl" of a scalar field ( , ): where the "2D curl" is defined as ∇ ×= . Since divergence of curl is zero (i.e., ∇ ⋅ (∇ × ) = 0), Equation (9), together with approximation (10) reads: where the fact that and ̃ are constants in the relevant channel domain is taken into account. In the following manipulations of the species transport equation for the three calculation domains that constitute a gas part in a slice the following subscript notation is used to distinguish quantities in different calculation domains: subscript (as in e.g., , , ) pertains to the channel domain, subscript 1 pertains to the GDL#1 domain and subscript 2 to the GDL#2 domain.

GDL Domains
Assumption IX infers that the velocity of gas is in GDL has zero component along coordinate, i.e., = 0. Thus, Equation (12) in GDL simplifies to: where is the effective binary diffusion coefficient in GDL. Furthermore, gas-flow in GDL is governed by Darcy's law: where is the permeability of GDL and viscosity of the gas. Since Equation (14) is a form of Equation (11) with = 0 and = − the gas flow in GDL is fully defined by scalar field , in other words: gas-flow in GDL is potential. The definition region of Equation (13) is a " 1 by ℎ1" rectangle for GDL#1 domain and a " 2 by ℎ1" rectangle for GDL#2 domain as discernible from Figure 2. For devising analytic solutions in the GDL domains also boundary conditions at the four boundaries must be defined for each GDL domain. The requirement that there is no total gas-flow through walls and across the symmetry planes yields the boundary conditions of zero velocity component perpendicular to these boundaries: where = ℎ1 , = and = 0 denote the positions at: wall on top of GDL#1 domain, the symmetry plane on the left side of GDL#1 domain and the symmetry plane on the right side of GDL#2 domain respectively. Since there can be no species flux through these boundaries either there has to be, in addition to zero velocity component, also zero concentration gradient perpendicular to these boundaries: The molar flux of water vapour from MEA into GDL, schematically represented with the two narrower vertical blue arrows in Figure 2, defines two boundary conditions: where is volume fraction of the gaseous phase in the GDL, = 0 denotes the position at GDL/MEA interface and → and → are the molar fluxes of water exiting MEA and entering GDL#1 and GDL#2 domains respectively. According to Equation (4) the molar fluxes of water exchanged between MEA and GDL and the consumption of reactant at the MEA/GDL interface define the boundary condition for net molar gas flow at this interface: where and are the molar fluxes of reactant consumed at the interface between MEA and GDL#1 and GDL#2 domains respectively.
Species molar flux between GDL#1 and GDL#2, schematically represented with the blue horizontal arrow in Figure 2, defines the boundary condition at the border plane between the two GDL domains: and similarly, the continuity condition for the net molar gas flow: where = 1 denotes the position at the border between the two GDL domains, → is the molar flux of water vapour exiting GDL#2 and entering GDL#1 and → the gas velocity component perpendicular to the boundary between GDL#1 and GDL#2 evaluated at that boundary. GDL#2 exchanges species also with the channel domain, defining a boundary condition: where = ℎ1 denotes the position at the border between the GDL#1 and the channel domain and 1→ is the molar flux of water exiting GDL#1 domain and entering the channel domain at that channel/GDL boundary. The boundary condition for the velocity potential at that channel/GDL interface is obtained by assuming uniform pressure in the channel region and thus defining the value of at that boundary according to Equation (14): where = is the pressure in the adjacent channel domain. The 2D analytic general solution of species concentration distribution and the velocity potential distribution for every GDL domain is obtained as a linear combination of eigen functions of the ∇ operator. As defined in Subsection A.1.2 of Appendix A, there are two families of these eigen functions: -the first, defined by Equation (A31), deals with species transport induced by convective transport along the coordinate, and thus does not apply to GDL domains; -the second, defined by Equation (A41), deals with species transport induced by the boundary conditions at the boundaries of the 2D computational domains and applies also to GDL domains.
The general solution for each GDL domain (represented by Equations (A56) to (A60) in Appendix A) is obtained by treating the molar fluxes at its boundaries as input parameters. The full derivation is given in Subsection A.1 of Appendix A.

Channel Domain
Governing equations of gas flow in a channel are of Navier -Stokes type and thus generally the flow in channels is not potential i.e., the scalar field in Equation (11) is not zero. However, since in any case the non-potential part of the 2D velocity profile vanishes in Equation (12) the subsequent derivation of species transport equations in the channel domain neglects the non potential part of velocity field in the cross-sectional plane.
According to assumption VIII and IX the channel domain is the only domain type in which species transport along coordinate takes place thus the term ( ) in Equation (12) is nonzero only in the channel domain. The core of the numerical treatment of the gas flow in channel is the following finite-difference approximation for the derivative ( ): where = − is the depth of slice with and denoting the value of coordinate at the points where gas-flow enters and exits the channel domain respectively as indicated in Figure 2. Using approximation (25) in (12) yields: Equation (26) can be interpreted as a simple species transport equation in a 2D plane with source and sink terms where the source term is ( )| (i.e., the species inflow contribution) and the sink term is − ( )| (i.e., the species outflow contribution). This interpretation enables the species transport along the coordinate and in the plane to be treated within a channel domain with a single 2D differential equation: where the source term and the sink term , schematically represented by the yellow arrows in Figure 2, are expressed as: where ( , ) is the cross-sectional profile of the -component of velocity and superscript denotes the values from previous channel domain i.e., the channel domain in the upstream neighbour slice. The velocity profile always has, according to assumption X, the same dome-like shape of laminar flow (as found in Section 5.2) scaled by the average -component velocity in the channel domain in question: where is the average -component velocity (obtained according to 1D gas-flow Equation (85) or (86)) and Υ( , ) is the unit-less dome shape as given in [33]. The definition region of Equation (27) is, for the channel domain, a " 1 by ℎ2" rectangle whose four edges represent the boundaries of channel domain with: two walls (gray surfaces in Figure 2), one symmetry plane (translucent yellow surface in Figure 2) and one channel/GDL boundary (translucent dotted surface in Figure 2). There are two types of boundary conditions on the four edges: 1. The total molar gas flow and the molar flux of water vapour between the channel and the GDL, schematically represented with the wider vertical blue arrow in Figure 2, manifest as boundary conditions: and: where = ℎ1 denotes the position at the GDL/channel boundary and → is the gas velocity component perpendicular to the boundary between GDL#1 and channel evaluated at that boundary. (15) and (16) there is no flux through walls and across the symmetry plane, leading to the boundary condition for the other three edges:

Analogous to conditions
and: where = ℎ , = 2 and = 0 denote the position at the two walls and the symmetry plane respectively.
The 2D analytic general solutions of species concentration and velocity potential distribution for a channel domain are obtained as a linear combination of both families of eigen functions of the ∇ operator since species transport in the channel domain is induced both by the convective transport along coordinate and by the boundary condition at the boundary between channel domain and GDL#1 domain.
A general solution for a channel domain (represented by Equations (A53) to (A55) in Appendix A) is obtained by treating the average -component velocity of gas flow, the molar fluxes from GDL and the source term as input parameters. The full derivation is given in Subsection A.1 of Appendix A.

MEA Part
Contrary to the two gas parts that are each split onto three computational domains the whole MEA part is contained within one computational domain as discernible from Figure 1 and Figure 3. . Schematic representation of species production, consumption and transport in the MEA. The two ribbed purple surfaces represent the catalyst layers at the MEA/GDL interfaces. The brown arrows represent the species transport from/to cathode gas part, whilst cyan arrows indicate the species transport from/to the anode gas part and the grey arrow represents the proton flow across the membrane. Hydrogen and oxygen are consumed on the anode and cathode catalyst layer respectively. All the water is produced on the cathode catalyst layer and the molar flux of water traversing the membrane is equal to the flux entering the anode gas part.
In the MEA part, two types of physical processes are distinguished: • species transport across the membrane (water molecules and hydronium ions represented by the cyan "H 2 O" and the gray "H + " arrow in Figure 3); • species production and consumption by the electrochemical reactions on the catalysts that are on the top and bottom boundary of the MEA domain as observed in Figure 3.
The following two subsections thus deal with the species transport and electrochemical reactions respectively.

Species Transport across Membrane
The equations for membrane species transport are derived with membranes for low temperature PEMFCs in mind; however, as indicated in Section 3.5, they are also applicable to membranes for high temperature PEMFCs.
The following equations governing species transport across the membrane are expressed with respect to the orientation of the coordinate system indicated in Figure 3 which is such that = 0 represents the position at the anode catalyst layer and = the position at the cathode catalyst layer. According to assumption XIII the species transport in the membrane takes place only along coordinate, thus the dependency of variables with respect to coordinate is left out from notation in this subsection (recalling that within a slice the variables are also constant along coordinate).
Molar fluxes of water and protons are defined by equations taken from [34]: where membrane water content Λ is the concentration of water normalised with the sulphonic group concentration, is the concentration of dissociated protons, is the electroosmotic drag coefficient, is Faraday constant, is electric potential and and are the diffusion coefficients of water molecules and dissociated protons (in form of hydronium ions) respectively. According to assumption XII both these two diffusion coefficients are linear functions of the membrane water content with their coefficients of proportionality and depending only on temperature. The explicit notation of temperature dependence has been omitted due to the isothermal nature of the model. Following the assumption XV the concentration of dissociated protons is approximately equal to the membrane sulphonic group concentration • : meaning that any proton concentration gradient is negligibly small and thus the termin Equation (36) is neglected and is treated as constant along the coordinate, i.e.,: ≠ ( ). Following the derivation of [34] and taking into account this simplification of zero proton concentration gradient, Equations (35) and (36) can be written as: where the constants and are defined as: In steady state the continuity condition requires that: leading to the molar fluxes and being constant along the coordinate. In the previously defined coordinate system the membrane water content at the interface with the anode electrode and at the interface with the cathode electrode are expressed as Λ = Λ( = 0) and Λ = Λ( = ) respectively and the general solution for and Λ in differential Equations (39) and (40) can be expressed as: , the back diffusion flux of water, is also constant with respect to coordinate. Rearranging Equation (40), integrating it over the membrane thickness and taking into account that and are independent of gives: The integral on the left side of Equation (44) evaluates to: where Λ is the arithmetic mean of water content at the anode and cathode side, and the integral on the right side of Equation (44) represents the voltage drop due to the membrane's ohmic resistance Δ : Taking into account that the proton flux is proportional to the current density : Equation (44) can be formulated as the "Ohm's law" equation for charge transport across the membrane: where (Λ , Λ ) = ( ) is the membrane's specific ohmic resistance.

Electrochemical Reaction Kinetics
Following assumption XIX, the fuel cell operates at an uniform operational voltage ≠ ( , ).
is a sum of four contributions [36]: where Φ and Φ are anode and cathode Galvani potentials respectively and Δ and Δ are voltage losses due to ohmic resistance of proton transport in membrane and electron transport in GDL respectively. Generally all four depend on current density, however taking the definition of zero galvanic potential of a standard hydrogen electrode, following assumption XVII that neglects reaction kinetics effects on the anode and following assumption XIX that neglects electric resistance of GDL yields Φ ( , ) = Δ = 0. This simplifies Equation (49): where Φ = Φ (i.e., the subscript in Φ is from here on dropped for brevity). Current density, proportional to the rate of oxygen reduction reaction at the cathode catalyst, is given as a function of Φ by the Butler-Volmer equation: where Φ, , denotes a Butler-Volmer function, is the exchange current density defined as current per unit area of catalyst layer, and are the surface concentrations of reactant and product (in this case oxygen and water) respectively, is the stoichiometry ratio, is the faraday constant, is the electron transfer coefficient, is the reference concentration and Φ is the cathode open circuit Galvani potential at reference conditions (reference conditions mean that species participating in the cathode reaction assume reference concentration ). The variable represented by (Φ − Φ) is not to be confused with overvoltage which is the departure from the open circuit potential that would be exhibited at conditions equal to the actual conditions in the immediate surroundings of the electrode under operation.
Defining the current density directly by Butler-Volmer (or Tafel) Equation (51) is applicable only to very thin catalyst layers (CL). In [37] it is shown that when finite thickness of CL is considered the local current density additionally depends also on the local proton conductivity of the ionomer phase in the CL. According to Equation (40) the local proton conductivity is linearly proportional to the local ionomer water content, which means that the current density additionally depends also on Λ . In the limit of very small CL thickness this dependency vanishes and the function of current density reduces to Equation (51). The assumption of negligibly thin CL is typically made by the reduced dimensionality models, e.g., [28,31,32], it is made by the benchmarking CFD model and thus also used in the HAN-FC model. However, the function defining the current density in the HAN-FC model can easily be extended to account also for the effects of non-negligibly thick CL: instead of Equation (51) an equation for current density as a function of not only Φ, and but also Λ , as derived for example in [37], can be used. This can be done since the approach of derivative approximation and estimation iteration, described in Section 2.4.1., is applicable to any nonlinear function of any number of arguments.
Expressing Δ from Equation (50), substituting it into Equation (48) and assuming a given species concentration at boundaries (i.e., Λ , Λ , and are treated as parameters) the two unknowns and Φ are obtained by solving the set of the two transcendent equations: Current density in turn also defines the rates of species production and consumption: where 2 and are molar consumption/production per unit area of oxygen and water at cathode catalyst layer and 2 is that of hydrogen at the anode catalyst layer.
Since the species consumption/production takes place at the boundaries it is at this point worth defining boundary conditions for the MEA domain in the same fashion as was done for the gas part domains. As indicated by the cyan "H 2 O" arrow in Figure 3 the molar flux of water traversing the membrane is equal to the molar flux of water crossing the interface between the MEA and the anode gas part: where is molar flux of water across the interface between the MEA and the anode gas part and is, with respect to the orientation of coordinate system in Figure 3, defined as positive when pointing in the negative direction. Using Equations (43) and (47) in Equation (56) leads to the boundary condition at the interface between MEA and anode GDL: The steady state condition requires that all water produced at the cathode catalyst layer diffuses to the anode and the cathode gas parts thus: where is molar flux of water across the interface between the MEA and the cathode gas part. Equations (54), (57) and (58) lead to the boundary condition at the interface between the MEA and cathode the GDL:

Coupling of General Solutions in Domains
Coupling of general solutions in domains is done by fulfilling the conditions of continuity of species activity, of species molar fluxes, of net molar gas flows and of pressure at those boundaries that are shared by two domains. In the gas part the continuity of species activity simply requires continuity of species concentration. The two couplings between gas part domains are thus straight forward: At the boundary between the GDL#1 and the GDL#2, domain continuity of activity requires that: continuity of molar flux, using Equation (21), requires that: and continuity of net molar gas flow, using Equation (14), requires that: continuity of pressure, using Equation (22), requires that: Similarly, using Equations (23), (31) and (32), the coupling between GDL#1 and channel domain requires that: where the continuity of pressure is already comprehended in the boundary condition (24).
Coupling the MEA domain to the GDL domains has to take into account the fact that water in membrane is in absorbed liquid phase while in GDL it is in gaseous phase. Activity of water vapour is simply its relative humidity: where is concentration of saturated water vapor, while activity of absorbed water in membrane follows the relationship from [23]: Assumption XVI requires continuity of water activity also at the membrane/GDL interface meaning that the activity of liquid water in membrane at that interface is equal to the activity of water vapour in GDL at that interface. This equality together with Equation (67) leads to the activity continuity condition at the boundary between the MEA domain and the domains of anode and cathode GDL: where ( ) and ( ) are water activity at the anode and cathode MEA/GDL interface and and in superscript or subscript denote values pertaining to the anode and the cathode side respectively. The continuity of molar fluxes of water vapour from MEA to the anode and the cathode side, defined as: and continuity of molar fluxes of reactant consumed at the anode and cathode defined as: lead, according to Equations (57), (59), (17) and (18) The main governing equations and corresponding boundary conditions are summarised in Table 1.

Electrochemical kinetics
Similarly, the membrane ohmic resistance can, according to Equations (69) and (70), be expressed as a function of and i.e., (Λ , Λ ) = Λ( ), Λ( ) = ( , ). The system of Equations (52) thus reads: Equations (75) to (80) therefore couple the species transport in membrane and electrochemical reactions on the cathode catalyst to species transport in both gas parts.

Solution for the Whole Cell
A full 3D solution for the whole representative unit is obtained as series of consecutive 2D analytic solutions for slices. Each solution for a slice is coupled to the solution for the neighbouring upstream slice by taking the anode and cathode sink terms of the upstream slice as the source terms for the slice in question. Obtaining the analytic 2D solution for a slice requires knowing the value of the mean gas velocity component in the sink term defined in Equation (29). The mean gas velocity component is calculated by equations of 1D gas-flow given in the Subsection 3.4.2. Description of obtaining the analytic 2D solution for a slice is given in the following Subsection 3.4.1.

2D Analytic Solution for a Slice
A general solution of the 2D diffusion problem in each of the seven computational domains of a given slice (Figure 1d) is obtained by finding the domain specific eigen functions (also called harmonics) of the ∇ operator. These harmonics are 2D trigonometric functions characterised by mode numbers and : numbering modes along coordinate and the modes along coordinate. (Fully detailed construction of the eigen functions is given in Appendix A). For each domain the full solution comes in the form of a specific linear combination of its eigen functions, i.e., a type of 2D Fourier series (Equations (A53) to (A60) in Appendix A). Finding the full solution for water concentration and velocity potential in a slice thus means finding all Fourier coefficients in the linear combinations of eigen functions of the seven computational domains. These Fourier coefficients are defined by the couplings between domains. Two types of coupling are distinguished: coupling between domains within a slice and coupling between channel domains in consecutive slices.
• Within a slice the coupling of two domains at their common boundary via the coupling conditions described in Section 3.3 is, on the level of eigen functions, manifested in algebraic coupling of the two sets of Fourier coefficients pertaining to the two linear combinations of harmonics in the two domains as described in Section A.3.1 of Appendix A.
• Coupling to the neighbour upstream slice via the two source terms in cathode and anode channel domain fully defines the solution for a slice (at a given operational voltage). Coupling to the downstream slice has no influence on the solution in the slice in question since sink terms are fully defined by concentration distribution in the two channel domains and are, due to assumption VIII, not influenced by the conditions in the neighbour downstream slice. Source terms, defined in Equation (28), also come in the form of a linear combination of eigen functions [Equation (A27)] thus the seven linear combinations that constitute the full solution for a slice are an algebraic function of the two linear combinations of the two source terms.
Equations (27) and (13) that govern species transport in gas parts are linear differential equations and, within a slice, the couplings between linear combinations of eigen functions of the coupled domains of gas parts are handled by linear algebra as described in Section A.4 of Appendix A. Contrary to this, the couplings between MEA and the two gas parts contain the electrochemical equations and the relationship of activity to water content dependence which are nonlinear. In order to handle also these nonlinear phenomena with linear algebra the approach of derivative approximation and estimation-iteration, described in the following paragraph, is used.
This approach assumes that within a slice each argument of these nonlinear functions (i.e., Φ( ), ( ) or ( )) has only small deviations from its average along the coordinate. Let stand for any of the aforementioned nonlinear functions and arg , arg , … stand for its arguments: Each argument is split into two parts: where arg stands for any of the arguments, Π is a predictor of the average of arg over (not necessarily the average itself but a close enough estimate of it) and Δarg( ) is its deviation from the predictor Π along .
The assumption of small deviations means that Δarg( ) is sufficiently small to justify the following derivative approximation: Additionally, the average species concentrations in the corresponding domains, required by Equations (13), (27) and others, are calculated in this step. This procedure of linearising the problem using predictors and re-evaluating the predictors using the solution of the linearised problem is represented in the flowchart in Figure 4 as the repetition loop that follows a negative (red box) test of convergence criteria (yellow box). This loop is iterated until sufficient convergence i.e., until positive outcome (green box) of convergence criteria test. It has to be noted that the derivative approximation is done only on the level of a slice requiring that the variation of arguments of nonlinear functions is small only along the coordinate. Variation of these arguments along the length of the representative unit ( direction) need not be small since the nonlinear variation from slice to slice is fully captured through the estimation-iteration procedure.

1D Gas-Flow Equations
Within the iterative procedure described in previous subsection and depicted in Figure 4 the channel domain mean gas velocity , needed in Equation (30), is also repeatedly re-evaluated (purple box) according to the following equations: where is the gas -component velocity of gas in channel domain in previous upstream slice, = 1 × ℎ2 , as discernible from Figure 2, is the area of channel cross-section ( = 1 × ℎ2 for anode channel and = 1 × ℎ2 for cathode) and is the net molar flux added or taken from the gas flow within the slice depth defined as: where = 1 + 2 = 1 + 2 is the width of a slice, and are the molar fluxes of water from the membrane to the anode and cathode gas parts as defined by Equations (57) and (59) and and are the changes in molar flux of hydrogen and oxygen due to consumption at the catalyst as defined by Equations (53) and (55). Once a converged solution for a slice is obtained, the sink terms of the anode and cathode channel domains are passed on to the next downstream slice to serve as the two source terms.

HAN Concept Extendibility
Assuming no liquid water in the channels and GDL the model is readily applicable to high temperature PEM fuel cells (HT-PEM) that typically feature only single phase flow. This is owed to the fact that PBI membranes that are used in HT-PEMs feature linear relationship between conductivity and relative humidity [38] and no electroosmotic drag [39], making the equations for membrane species transport as given in Section 3.2 directly applicable also to the PBI membranes.
The presented framework of solving the diffusion equation on rectangular definition regions can be applied to treat also other phenomena governed by diffusion-type mechanisms such as heat transport, multi-component diffusion and Darcy-type flow of liquid water induced by gradients of capillary pressure, leading to a comprehensive treatment of also low temperature PEM fuel cells.
The principles of HAN model have been demonstrated on a simple straight channel fuel cell, however, the HAN model can be further extended to model also fuel cells with more complex geometries. In serpentine channel geometries the pressure difference between neighbor channel menders leads to a bulk gas cross-flow between the meanders, i.e., nonzero flow across the yellow plane in Figures 1 and 2 (which is thus no longer a symmetry plane). Having a routine for calculating the pressure drop along the channel as given in [33] and noting that the channel pressure directly translates to the velocity potential boundary condition, as defined by Equation (24), the velocity field of this cross-flow is obtained using the formalism in Section 2.1. This leads to coupling each slice not only to its upstream neighbor but also to its two side neighbor slices.

Simulation Setup
To make the simulation setup of HAN-FC and CFD models as identical as possible all material and operational parameters, with the exception of exchange current density, were set to the same values in both HAN-FC and CFD model settings. The parameter defining the effective exchange current density, for which the user manual of the CFD software gives insufficient information, was obtained through calibration procedure performed on a case of higher operational voltage (low current densities) while the comparative analysis was performed for a case of lower operational voltage (high current densities). Calibration of HAN-FC for obtaining the parameter of effective exchange current density is described in Section 5.1 while other parameters used in HAN-FC and CFD models are summarised in Table 2. * In CFD this is the value of outlet pressure, maximal deviations from this value anywhere in the representative unit are +0.53% and −0.28% justifying assumption V of isobaric HAN-FC; ** In CFD these are the mean values since the CFD simulation takes into account the dependence of binary diffusion coefficient on the ratio of two the components. However, the relative standard deviation is below 0.09% with maximal deviations from the mean value anywhere in the representative unit being +0.45%and −0.26% justifying the assumption VI, i.e., the use of component ratio independent diffusion coefficient in HAN-FC.
The gas composition was, according to assumption VI, taken as bi-componential, i.e., H 2 with H 2 O (g) for the anode and O 2 with H 2 O (g) for the cathode. The parameters for membrane assumed no gas crossover.
CFD simulation was done using AVL FIRE ® v2011 with the "Fuel cell" module (validation of which is found in references [40] and [41]). Two meshes differing in meshing density shown in Figure 5 were used for simulations. The denser mesh was used for obtaining results with high spatial resolution needed for plotting smooth 3D graphs in Section 5. The coarser mesh was used for estimating the shortest computational times of a CFD simulation that still gives accurate results on current density. Meshing of the representative unit for CFD simulation shown in Figure 5 was done as follows (numbers apply to the denser/coarser mesh): the mesh of the representative unit is made of anode and cathode side meshes containing equal number of subdivisions. Both GDLs are thus sectioned into 12/6 sections along GDL thickness ℎ1, 12/6 along the width of representative unit and both channels are sectioned into 12/6 sections along channel height ℎ2, 6/3 along the channel width 1 and the whole model is sectioned into 80/40 sections along the direction of gas flow (i.e., along the length of the representative unit ). Meshing of the membrane is done internally by the AVL FIRE fuel cell module as described in [42] and [33]. AVL FIRE fuel cell module treats the catalyst layers as zero-thickness surfaces. The isothermal conditions were obtained by omitting calculation of energy equation. • In HAN 6 first 6 modes per each coordinate are taken into account in harmonics in each gas part computational domain.
• In HAN 9 first 9 modes per each coordinate are taken into account in harmonics in each gas part computational domain.
Taking into account a greater number of harmonics brings higher accuracy but also requires longer computational times as presented at the end of Subsection 5.2 and in Subsection 5.3.

Calibration and Results
HAN-FC was calibrated by finding the value of effective exchange current density that leads to the best agreement between HAN-FC and CFD results on current density at operational voltage of 0.85 V. The HAN-FC model with the so obtained exchange current density parameter was then evaluated globally by a comparative plot of the polarisation curve and locally by a detailed comparative analysis of spatial distribution of key variables at the operating point of highest current.
To aid the interpretation of the graphs in the following figures. Figure 6 shows the modelled representative unit placed in the coordinate system with respect to which all graphs of CFD and HAN-FC results are plotted. It should be noted that the choice of directions and origins of coordinates in Figure 6 is such that it enables easy interpretation of graphs in Section 5.2 and differs from the choice in Figure 2 or Figure 3 where the coordinate directions and origins are chosen in such a way to best suit expressions of governing equations in Section 3. Additionally it should be noted that the gap between anode and cathode GDL in Figure 6 does not reflect the thickness of the membrane but only the positioning of the meshes of the anode and cathode side in the CFD model as seen in Figure 5. Figure 6. Representative unit in coordinate system. Cathode channel is coloured blue, cathode GDL magenta, anode GDL green and anode channel red. The central symmetry plane of a rib and the symmetry plane between ribs (respectively depicted as green and yellowish surfaces in Figure 1) that define representative unit are at y = 0 and at y = 0.75 mm respectively. Anode and cathode inlets are at, z = −13.5 mm outlets at z = 13.5 mm and the cross-sectional plane midway between inlet and outlet is at z = 0.

Calibration
At operational voltages of 0.85 V and lower the Butler-Volmer Equation (51)  Finding the best fit of HAN-FC to CFD is done by adjusting the parameter 0 * in Equation (91) which avoids the need to define the exact value of Φ . The best fit of current density results of HAN-FC to those of CFD, graphically depicted in Figure 7, is found at * = 0.1195 / .

Plots of Comparative Results
In this subsection a number of comparative results are presented where in all graphs the CFD results are coloured green and are obtained using the denser mesh and HAN-FC results are coloured brown and are obtained using HAN 6 model. At the end of this subsection also a plot comparing results of HAN 3 and HAN 6 is shown. Figure 8 shows the polarisation curve as calculated by the CFD and by the HAN-FC where the points representing the CFD results are obtained by calculating current density at five predefined operational voltages (with all other operational parameters fixed) whereas the points representing the HAN-FC results are obtained by finding the values of voltages that give the same current densities as the CFD results. The lowest operational voltage chosen for CFD simulation was 0.75 V, since at lower voltages liquid water, which cannot be comprehended by the HAN-FC model, starts occurring at the outlets. Comparison of Figure 7 and Figure 9 shows that the operational voltage of 0.75 V, i.e., the operational point with highest current density, leads to considerably different operational conditions compared to the calibration case of 0.85 V. All following analyses are done for the operational point with the highest current density that fully exposes all governing phenomena that are comprehended in governing equations of Section 3. These analyses are thus aimed at validating the capability of HAN-FC to accurately treat these governing phenomena in a general straight parallel channel fuel cell.
As the most important simulation result Figure 9 shows graph of current density distribution. The primary objective of a fuel cell simulation is obtaining the net eclectic current at given operational conditions and in this respect the close agreement between HAN-FC and CFD results in Figure 9 speaks of high level of fidelity of the HAN-FC model. However it is instructive to analyse the underlying physical phenomena responsible for the shape of plots in Figure 9 and validate HAN-FC model under higher level of scrutiny.
The trend of current density rising with increasing values of (i.e., from inlet to outlet), as discernible from Figures 9a,b, Figure 10) and rises with decreasing cathode Galvani potential Φ (plotted in Figure 11). Since the dependence of current density on is linear and dependence on Φ is exponential the influence of decrease of Galvani potential measuring a few hundredths of a volt, as observed in Figure 11, has a much larger effect than the decrease in oxygen concentration of around ten percent as observed in Figure 10. The decrease of Galvani potential is, according to Equation (50), equal to decrease of membrane overpotential. A decrease in membrane overpotential comes as a consequence of decreased membrane ohmic resistance that is inversely proportional to the mean membrane water content [Equation (48)]. Figure 12 shows that the two gas-flows get increasingly more saturated with water vapour on their way from inlet to outlet and consequentially also the membrane water content rises from inlet to outlet as discernible from Figure 13. Altogether this means that the increase in membrane hydration increases current density and indeed the plots of mean water content in Figure 12 and current density in Figure 9a are very similar. The trend of the increasing current density on the way from inlet to outlet thus reflects the influence of membrane hydration on fuel cell performance.    Additionally, Figure 13 reveals that the membrane water content rises slightly also with increasing values of y. This is explained with the help of Figure 2: Water diffusing from MEA into GDL at lower values enters GDL#1 domain (as indicated by the left narrower vertical blue arrow) and from there it diffuses into the channel domain (as indicated by the wider vertical blue arrow) where the stream of gas-flow caries it away. Water diffusing into GDL at higher values enters GDL#2 domain (as indicated by the right narrower vertical blue arrow) and from there it first has to diffuse into GDL#1 domain (as indicated by the horizontal blue arrow) before reaching the gas stream in channel domain. This leads to a higher water vapour build-up in the GDL#2 domains compared to GDL#1 domains and consequently regions of membrane adjacent to GDL#2 domains feature higher water content than regions adjacent to GDL#1 domains. This in turn also explains the higher current densities observed at higher values of in Figure 9c.
All comparative graphs of the MEA variables in the plane presented in Figure 9, 10, 11 and 13 show close agreement between HAN-FC and CFD results and thus show important agreement of results on distribution of cathode Galvani potential and mean membrane water content, shown in Figures 11 and 13, that directly reflect the nonlinear phenomena of membrane ohmic resistance, electrochemical reaction kinetics and membrane water activity [reflected for example in Equations (48), (51) and (68)]. This validates the capability of HAN-FC to accurately treat the coupled linear and nonlinear phenomena through the estimation-iteration loop described in Section 3.4.1. Furthermore Figure 9b and especially Figure 9c show that HAN-FC with its derivative approximation is also capable of accurately capturing variation of variables along the y coordinate.
It is also valuable to comparatively validate HAN-FC with results on species concentration distribution presented in Figures 12 and 14. The brown plot in Figure 14 represents HAN-FC's full 2D analytic solution for water vapour concentration distribution in the slice that is midway between inlet and outlet. Figure 15 shows the component of net molar gas velocity (i.e., component perpendicular to membrane) at the GDL/MEA interface. The discontinuity and the notable ondulations in the HAN-FC plot in Figure 15 are inherent to the domain approach and the finite number of harmonics used in the HAN-FC model. The comparative graphs in Figure 14 and 15 validate the approach of splitting diffusion problem into computational domains and finding solutions as linear combinations of eigen functions. Figure 12 shows the species concentration distribution on the rib's central symmetry plane and thus shows the variations along the length of the representative unit. The variation of HAN-FC's results along the length of representative unit reflects the variation of the 2D analytic solution from slice to slice. Since 2D solutions in consecutive slices depend on the coupling between slices via sink and source terms the close agreement between HAN-FC and CFD in Figure 12 validates the approach to treating the convective gas-flow with source and sink terms in conjunction with the 1D bulk gas-flow model.  As a CFD simulation result Figure 16 shows a profile of z-component velocity of cathode gas on the cross section that is midway between inlet and outlet. It is discernible from this plot that gas z-component velocity in GDL is negligible and that the velocity profile in the channel assumes the dome shape of laminar flow. This shape of the velocity profile is the same for any cross section meaning that only the height of the dome shape varies and that z-component velocity in GDL is always negligibly small. Since the same also holds true for anode gas, this justifies assumption IX.
Overall Figures 9-15 show a very good agreement between CFD and HAN-FC in terms of trends in analysed variables and in terms of their absolute values. It has been shown that by calibrating the model on a case where the resulting current density is relatively uniform (Figure 7), very good agreement with CFD results were achieved for a case where the variation of current density is substantial (Figure 9), confirming the capability of HAN-FC to model coupled 3D species transport and electrochemical phenomena with high fidelity. Figure 16. Cross-sectional profile of absolute velocity of cathode gas at cross-section midway between inlet and outlet as obtained by CFD simulation. Figure 17 shows a comparison between two analytic HAN-FC solutions for a slice where one takes into account three modes along each coordinate within a domain and the other six modes. The two analytic solutions are hardly distinguishable from one another with the difference in accuracy being noticeable only at the interfaces between computational domains where the analytic solutions for domains are sewed together. This means that as little as three modes per coordinate taken into account in harmonics of each computational domain already capture the essential features of the concentration distribution. Furthermore plots of HAN-FC results in Figures 9-15 stay practically identical when plotted using results obtained with HAN 3 instead of HAN 6 indicating that HAN-FC can give accurate results at low computational demands. Figure 17. Comparison between HAN 6 (brown) and HAN 3 (yellow). Close-up of the concentration distribution in cathode GDL domains and part of cathode channel domain of the slice that is midway between inlet and outlet is shown. Table 3 shows computational times for models with different spatial resolution. The spatial resolution in case of CFD models is defined by the number of volume elements of the mesh whereas the spatial resolution of HAN-FC model is defined by the number of harmonic modes per coordinate taken into account in each computational domain.  3 3 modes 1 1 0.85 0.85 HAN 6 6 modes 1 1 3.7 3.7 HAN 9 9 modes 1 1 6.0 6.0 * Iterations of obtaining solution for whole representative unit.

Computational Times
The denser CFD mesh features twice as many divisions along each coordinate as coarser one thus featuring eight times more volume elements. Larger number of volume elements also requires more iterations to reach adequate convergence. The large difference between computational times of HAN-FC and those of CFD originates from the fact that CFD requires from a few hundred to a few thousand iterations of the solution for the whole mesh to reach convergent results, whereas the isothermal and isobaric HAN-FC calculates everything in one go. This single iteration of the solution for the whole representative unit is not to be confused with iterations performed within each slice that are performed by the estimation-iteration routine for treating nonlinear functions. Typically 5 to 6 of these iterations were required in each slice to reach relative accuracy of 10 −4 (convergence criteria in yellow box in Figure 4).
Recalling from Figure 17 the high accuracy of the lower resolution HAN-FC together with Table 3 reveals that HAN-FC gives accurate results at computational times that are two to three orders of magnitude shorter than that of CFD. The core of computational operations of HAN-FC is the linear algebraic solving for eigen function coefficients, described in Appendix A, which is in the case of simulations presented in this paper done using vector algebra as introduced in HAN-ST of [33]. Furthermore, the computational times of HAN-FC model evaluated in this paper have a considerable margin for additional improvement by programming the model in a programming language such as C or Fortran that are more suitable for numerical computation than Mathematica and by executing the simulation on multiple CPU cores.

Conclusions
The paper presented a Hybrid Analytic-Numerical approach to fuel cell modelling (HAN-FC) efficiently applied to modelling a realistic straight parallel channel fuel cell. The HAN-FC approach efficiently combines three main features: • The 1D numerical treatment of gas flow; • The analytic 2D solution for the species concentration profile in the plane perpendicular to the gas flow and • The electrochemical sub-model.
The key to the hybrid 3D solution is constructing it as a series of consecutive 2D analytic solution that are obtained by splitting the 2D diffusion problem onto separate computational domains which on one hand enables modelling of realistic fuel cell geometries and on the other hand gives computationally efficient results with high accuracy already at very moderate number of harmonic modes taken into account. The efficiency of incorporation of the electrochemical sub-model lies in the derivative approximation and estimation-iteration approach and is reflected in the fact that electrochemical calculations do not need a separate routine but are performed alongside calculations of the nonlinear coupling of water activity in gas to the water activity in membrane. Overall, HAN-FC proves to be very accurate and computationally efficient and as such a very promising standalone fuel cell model for system level simulations. Further challenges in developing HAN-FC for system level simulations remain the treatment of liquid water in GDLs and channels.

Conflicts of Interest
The authors declare no conflict of interest.

Appendix A
A steady state solution for the whole representative unit is made of consecutive steady state solutions for each slice. In this Appendix a derivation of the full solution for a slice is presented. First, general solutions for the seven computational domains (indicated in Figure 1d) shall be constructed and then the full solution shall be obtained by applying appropriate coupling among the seven computational domains. The derivation makes the assumptions described in Section 2.

A.1. Gas Part
This subsection deals with obtaining the three general solutions for the three computational domains of the gas part, schematically depicted in Figure 2, and with coupling of these general solutions. First the species transport equation in each computational domain shall be split onto sub-equations, following with the construction of the two types of eigen functions: eigen functions dealing with species transport induced by the convective transport along the z coordinate (source and sink terms) and eigen functions dealing with species transport induced by the boundary conditions at the boundaries of the 2D computational domains. Finally general solutions shall be expressed in form of Fourier series as linear combinations of eigen functions.

A.1.1. Diffusion Differential Equation
2D species molar flux in the xy plane can be, analogously to Equation (1), defined as: where approximation (10) and neglect of any non-potential part of velocity field have been taken into account. By defining a new quantity ℤ: both species transport equations and total gas flow equations are translated into the same diffusion equation, namely: Gradients of both and ℤ define vector flows: which are the gas velocity that measures the net gas flow and the species molar flow respectively. Furthermore the zero velocity divergence requirement (8) in the GDL where the z-component velocity is zero ( = 0) reads: Equation (13) for species transport in the GDL transforms into the same form when expressed with ℤ: where the diffusion constant assumes value and steady state conditions are taken into account (i.e., = = 0 ). Similarly the expression for species transport in channel (27) becomes: Using Equation (8) and finite-difference approximation (26) (i.e., ≅ − ) the exact same form is obtained for velocity potential in the channel: where the corresponding source and sink terms are defined as = and = . Since both and ℤ obey the exact same equations, the following derivation of general solutions for gas part shall be shown only for ℤ. Furthermore, because Equations (A4), (A6) and (A7) are classical equations of the steady state diffusion problem with source and sink terms the quantity ℤ shall in the following be referred to as concentration and the corresponding species flux as diffusion.
In the following derivation the subscript position will be used for indexing eigen functions and their corresponding Fourier coefficients, thus the distributions of concentration ℤ and velocity potential in the channel domain,,GDL#1 domain and GDL#2 domain are renamed: Each computational domain in the gas part is coupled to the neighbour domains via the molar flux across the boundary reflected in the boundary conditions (17), (18), (21), (23) and (31). Additionally the channel domain is also coupled to the two channel domains in the upstream and downstream neighbour slices via the source and sink terms [Equations (28) and (29)]. Each of these couplings has its own contribution to the inflow (or outflow) of species into (out of) the domain. Thus the and similarly Equation (A6) that holds for GDL#1 and GDL#2 domains: where stands for either or ℂ (GDL#1 or GDL#2) and stands for any of the , , and ℎ representing the four boundaries. Equations (A14) and (A15) already reveal that and do not interfere with each other's influx coupling condition since assumes no outflow and thus does not add anything to the term and vice-versa for the . The three contributions to the concentration distribution in channel domain satisfy the following boundary conditions: where the latter expression acknowledges that , ( , ) = 1. The second last term in (A34) contains explicit time dependence which may seem to work against finding a steady state solution, however this equation applies to one concentration distribution contribution not the total concentration distribution and since a steady state solution of the total concentration distribution in a domain is sought this leads to the requirement that the time dependant terms of all the contributions must cancel each other out. Before defining this condition more concretely the other type of eigen functions needs to be constructed.

A.1.2.2. Eigen Functions for Flow from Boundary
The flux from boundary gives rise to six contributions (pertaining to the three different gas part domains): ( ) , ( ) , ( ) , ( ) , ℂ ( ) , ℂ ( ) . The six families of eigen functions behind this six contributions are practically the same differing only in their orientation or dimensions of their definition regions. Thus, as an example, only the family of eigen functions for ( ) will be derived in this subsection with the other five given in Appendix B. The with the definition region 0 < < ℎ1 ∧ 0 < < 1 as discernible from Figure 2. Satisfying conditions (A64), (A65) and (A66) is crucial for obtaining a physically meaningful steady state solution presented in Section A.3 that deals with algebraic manipulation. Identical derivation applies also to the total gas velocity field : equations for in the three gas part domains are obtained by substituting , , and ℂ with , ℬ and in the derivation.

A.2. MEA Domain
A general solution for water concentration and flux in membrane is obtained using Equations (42) and (43) taking into account that Λ , Λ are functions of ( ) and ( ): From Equations (A68) and (68) it follows that is linear in ( ) and nonlinear in ( ) and ( ). In order to couple to other variables using linear algebra the nonlinear dependence is linearised with the derivative approximation described in Section 3.4.1. The following paragraph gives derivative approximation specifically for . for > 0:

A.3. Coupling of Linear Combinations of Eigen Functions
Computational feasibility requires a finite number of harmonics taken into account in the Fourier series and series of , and . Greater number of harmonics taken into account improves accuracy but lengthens computational times. Let N +1 be the maximal number of modes in coordinate and also the maximal number of modes in coordinate, i.e., both index and index running from 0 to N (Generally the maximal number of and indices need not be the same, however this was chosen in this case for simplicity). This means the total number of , harmonics is (N + 1) 2 and the total number of and harmonics is N + 1. Thus, in sums of the following expressions the infinity as the upper summation limit shall be replaced by N (or 2N + 1) in case of MEA domain as explained in A.3.1).
There are two types of couplings: • Couplings between domains of a slice that require continuity of water activity and molar fluxes and • Couplings between channel domains in consecutive slices that require computing sink terms that serve as source terms in neighbour downstream slices.

A.3.1. Couplings between Domains of a Slice
The coupling of two domains sharing a boundary is done by equalising the expressions for water activity and molar flux in each of the two domains at their common boundary. The boundary values are expressed as Fourier series and thus, the coupling is done by equalising the corresponding Fourier coefficients.
First let us deal with coupling of the molar fluxes which is due to construction of very straightforward. At the boundary between channel domain and GDL#1 domain the continuity condition, according Equations (A23) and (A24), is expressed as: The corresponding boundary condition for the velocity field is, according to Equations (19), (20) and (71) In order for Equations (A88) to have a unique nontrivial solution the number of terms must be the sum of the number of , terms and the number of ℂ, terms, i.e., index must run over twice as many values as index .
Coupling the water activity is less straightforward. First let us deal with continuity of water activity. Water vapour concentration in the domains of the gas parts is, according to Equation (A2) obatined as: At the boundary between channel domain and GDL#1 domain the water vapour concentration (which is proportional to water activity) is from the two sides expressed as: The pressure continuity conditions that apply to the net gas velocity field are somewhat simpler. Following Equation (A63) continuity of pressure at the boundary between the two GDL domains requires that: The fact that zero can be put on the right hand side of Equation (A107) comes from the fact that the velocity field potential is defined to an arbitrary constant, i.e. = ∇ = ∇( + ). Having uniquely defined velocity potential at the GDL/channel boundary no other additional coupling for velocity potential needs to be defined.

A.3.2. Coupling between Slices
Coupling between slices is done by taking the sink terms of upstream slice to serve as the source terms of the slice in question. This is done by having a general expression for the sink term for anode and cathode side and evaluating it once the full 2D solution for the slice in question is obtained and passing it to the next downstream slice.