## Abstract

The gas diffusion layer in the electrode of a proton exchange membrane fuel cell is a highly porous material which acts to distribute reactant gases uniformly to the active catalyst sites. We analyse the conservation laws governing the multiphase flow of liquid, gas and heat within the electrode. The model is comprised of five nonlinear-degenerate parabolic differential equations strongly coupled through liquid–gas phase change. We identify a scaling regime in which the model reduces to a free boundary problem for a moving two-phase interface. On each side of the moving boundary the nonlinear system is well approximated by its linearization whose relaxation times are much shorter than the front evolution. Using a quasi-steady reduction, we obtain an explicit leading-order evolution equation for the free surface in terms of the prescribed boundary conditions.

## 1. Introduction

Polymer electrolyte membrane (PEM) fuel cells are unique energy conversion devices, efficiently generating useful electric voltage from chemical reactants without combustion. They have recently captured public attention for automotive applications for which they promise high performance without the pollutants associated with combustion. The operation of PEM fuel cells involves the convective transport of gases in flow field channels, convective-diffusive transport through the gas diffusion layer (GDL), electro-chemical reaction in both the anode (fuel/negative) and cathode (air/positive) catalyst layers, conduction of protons through the PEM, conduction of electrons through the GDL and graphite plates, and removal of product heat through the coolant channels, see figure 1. The membrane's conductivity improves significantly with hydration level, which in turn is sensitive to the relative humidity of, and the presence of liquid water in, its environment. However, excessive production of liquid water can flood the porous catalyst layer leading to a drop in power. Water management is the balance between maintaining sufficient hydration levels in the membrane without flooding the gas conducting regions, and it requires an understanding of when and where liquid water is produced and how it is transported within the fuel cell.

There has been considerable attention given to the development of two- and three-dimensional numerical models to approximate the coupled two-phase mass and heat transport within unit cells, see He *et al*. (2000), Hsing & Futerko (2000), Um *et al*. (2000), Natarajan & Van Nguyen (2001), Rowe & Li (2001), Wang *et al*. (2001), Berning *et al*. (2002), Bradean *et al*. (2002), Djilali & Lu (2002), Pasaogullari & Wang (2004), Siegel *et al*. (2004) and Um & Wang (2004). These approaches have either treated the multiphase flow in a trivial way, advecting the liquid as suspended droplets in the gas flow, or used isothermal models in which the vapour and liquid are treated as a combined mixture, or have experienced convergence difficulties for physically reasonable values of the parameters. The majority of the treatments are at steady-state, and none have treated the transient evolution of the free surface between the two-phase and dry regions in the hydrophobic GDL which arises from the degenerate nature of the liquid water transport. A difficulty in treating the evolution of these surfaces is the dichotomy of time-scales present in the system. Even in the single phase regime, previous analysis (Promislow & Stockie 2001) has shown that the convective and diffusive gas transport relaxes on time-scales separated by 3–4 orders of magnitude. The inclusion of the slow liquid transients stretches the time-scale dichotomy from 10^{−6} s for pressure relaxation to 10^{2}–10^{3} s for motion of the free interface.

A second difficulty in numerical computations is the length-scale dichotomy. The large aspect ratio of the fuel cell, roughly 1000–1 along the channel versus through the membrane electrode assembly (MEA), renders the resolution of thin fronts within the through-plane direction of the MEA a significant obstacle to numerical stability. Several authors (Divisek *et al*. 2003; Freunberger *et al*. 2003; Kulikovsky 2003; Berg *et al*. 2004), have exploited this high aspect ratio, proposing so-called 1+1D models, which couple one-dimensional systems for the through-MEA transport (*z*-direction) to one-dimensional along-the-channel convective transport (*y*-direction) in the flow fields, see figure 1*b*. The through-MEA models provide source terms for the convective channel flow, consuming oxygen and producing water vapour as the gas flows from inlet to outlet. In turn, the gas flow provides the boundary conditions for the through-MEA codes. Such models have had success in determining the local current and voltage, and hence power generation, of the cell, however, they have not incorporated multiphase aspects of water management due to its significant computational complexity.

The analysis presented here affords an *analytic* description of the impact of multiphase flow on fuel cell performance, greatly extending the accuracy and applicability of previous models without significant increase in computational expense. Indeed, the inclusion of multiphase water management in models of stacks of 50–200 unit cells used for automotive applications, see Berg *et al*. (in press), would seem certain to require a reduction along these lines.

We present a model of gas, liquid and heat transport through a one-dimensional slice of the GDL of a PEM fuel cell. We prescribe the gas molar concentrations and temperature in the flow field channel, as well as the flux of product water and heat in the catalyst layer. We develop a class of quasi-steady solutions which describe the partition of product water into liquid and vapour fluxes in a boundary layer at the cathode GDL-catalyst layer interface, which resolve the slow time-scale evolution of the two-phase free-boundary within the GDL and the heat and gas flux into the flow field channels, and which predict the time-scales and operating conditions needed to fill or empty the GDL of liquid water.

We consider a one-dimensional slice in which part of the GDL is dry and part is two-phase, as depicted in figure 1*b*. We construct a class of solutions in which the liquid water in the two-phase region is above the immobile volume fraction, *β*_{*}, and the gas phase is nearly saturated. Phase change represents a singular perturbation to the system, as does the transition of the water volume fraction, *β*, from mobile (*β*>*β*_{*}) to immobile (*β*<*β*_{*}) volume fractions. We have two distinct phase-change induced boundary layers, a fixed one at the catalyst layer–GDL interface arising from the electrochemical reactions, and a free boundary separating the dry and two-phase regions of the GDL. We do not resolve the fixed layer, rather the fluxes of heat and total water (liquid+vapour) are prescribed on the catalyst side of the fixed boundary, and sufficient phase change is imposed within the fixed layer to bring the water vapour flux into balance with the vapour flux carried by the saturated two-phase gas, taking into account the impact of phase change on the heat flux.

The free boundary layer has a sub-layer structure. In the evaporative layer from , see figure 2, the liquid water volume fraction is greater than the immobile volume fraction, the gas phase is under-saturated, and the amount of phase change and the associated jump in fluxes are determined. In the infinitesimal sub-layer which has left and right limits and , respectively, the water volume fraction reaches the immobile volume fraction and Rankine–Hugoniot conditions, essentially conservation of mass across the interface, are enforced which ensure the dry-phase and two-phase solutions form a weak solution across the entire domain.

We non-dimensionalize the transport equations, which reveals three major factors affording simplification of the analysis. The current density represents a small molar flux, as a consequence the conservation laws governing heat and gas transport can be treated as linear on either side of the free boundary with little loss of accuracy. The capillary pressure dominates the gas pressure yielding a liquid water distribution which is a linear function of a liquid water potential to good approximation. The large separation of time-scales for the relaxation of the heat, gas and liquid transport from the time-scale of motion of the free boundary permits a quasi-steady approximation in which the PDE dynamics are taken at steady-state, adiabatically driven by the slowly moving free boundary.

The model we present applies to a wide range of PEM fuel cell operational regimes; there are some caveats to its application. We assume the flux of gas and heat out of the catalyst layer scale with the current density. This would not be the case, for example, with a dry anode feed and a wet cathode feed, which would generate water transport independent of the current level. We assume there is no lateral pressure gradient imposed in the *x*-direction across the GDL, such as would arise in interdigitated or serpentine flow fields: we consider straight flow fields. We assume the channel gas is undersaturated but nearly saturated as outlined in (5.13) and that the gas phase is undersaturated up to and into the free boundary. If the gas becomes oversaturated at the dry edge of the free boundary then the resultant condensation could form a region of immobile water, 0<*β*<*β*_{*}. We preclude this case by restricting our results to those situations where the gas phase undersaturation, given by (5.15), remains positive.

The main result is summarized in §5. We construct a leading order approximation to a family of sharp-interface, weak solutions of the governing system of conservation laws, along the lines of Fife (1988). The family describes the slow evolution of the free boundary at the two-phase interface, and analytically connects the prescribed channel values and heat and water production in the catalyst layer to both the production of liquid water in the catalyst layer and the fluxes of heat and mass into the channel flow stream. In addition, we show that to leading order the phase change is localized at the fixed boundary layer at the catalyst–GDL interface and the free boundary. The partitioning of product water into liquid and vapour is sensitive to the thermal conductivity of the solid phase of the GDL. Increasing this coefficient by factors of 2–10 has a significant impact on the water partitioning, and can *reverse* the direction of vapour phase flow, driving vapour from the channel to the cathode catalyst layer, see (4.49) and following discussion.

The paper is organized as follows. In §2, we present the notation and a dimensional version of the model, which is non-dimensionalized in §3 and put into system form in (3.19). The cathode boundary conditions which describe the product water partitioning are given in §2*c*. In §4, we rescale and linearize the governing system for the dry, two-phase and free boundary regimes. We present both fully time-dependent and quasi-steady reductions for the dry and two-phase regimes. Within the limits of the quasi-steady reductions, we derive explicit leading-order solutions. We determine the jumps across both the resolved portion and the infinitesimal sub-layer of the free boundary and obtain a differential equation for the motion of the free boundary in terms of the product water partition and the gas phase undersaturation. In §5, we derive the sharp-interface reduction for the overall problem, using the quasi-steady solutions and jumps across the free boundary to determine an expression for the gas-phase undersaturation. This yields a closed differential equation for the free boundary motion.

## 2. Presentation of the model

### (a) Notation

The free boundary layer is demarcated by the points , , , see figure 2. The evaporative boundary layer is the region , while the infinitesimal sub-layer has left and right limits and . The jump of a quantity across the resolved boundary layer is denoted by(2.1)the jump across the sub-layer is denoted by(2.2)while the jump across the whole free surface is denoted by(2.3)The *j*th component of a vector * V* is denoted

**V**_{j}. The subscripts and

*z*denote differentiation with respect to that variable. Except for the scaling parameters

*t*

_{*},

*I*

_{*}and

*T*

_{*}, the sub- or superscript ‘*’ denotes a quantity evaluated at the two-phase boundary point

*z*

_{*}. The superscript ‘t’ after a row-vector indicates transpose. The projections are given by and . The positive

*z*-direction is from channel to membrane, positive fluxes are channel to membrane.

### (b) The GDL model

The dimensional unknowns are the molar concentrations of oxygen, water vapour and nitrogen, denoted by , and , the liquid water volume fraction *β*, and the temperature . The total (convective plus diffusive) flux of species *α* is denoted by *N*_{α}. The variables represent averages over the pore-scale structure of the liquid–gas–solid domains. The gas phase occupies the fraction 1−*β* of the void space, which is *ϵ*=0.74 of the total volume fraction. Conservation of mass and energy is expressed through the following PDEs which are coupled not only through the constitutive relations for the fluxes, but also through the condensation, , which appears as a source or sink term:(2.4)(2.5)(2.6)(2.7)(2.8)The total flux of gas species *α* is the sum of the convective and diffusive fluxes(2.9)As a convenience, we have assumed a Fickian formulation for the diffusive flux, with a common diffusivity, *D*, for all gas species. It is straightforward to incorporate a Maxwell–Stefan formulation, see Taylor & Krishna (1993) and Stockie *et al*. (2003), for the diffusive gas fluxes into the analysis presented here.

The gas pressure is given by the ideal gas law,(2.10)where the total gas concentration satisfies . The gas velocity is given by Darcy's law,(2.11)where *K* is the permeability, *μ*_{g} is the gas viscosity and the gas relative permeability, *k*_{rg} is given by(2.12)

An essential feature of the multiphase model are the constitutive relations, which specify the liquid relative permeability and the capillary pressure in terms of the liquid volume fractions. These relations incorporate the considerable role played by the hydrophobicity of the GDL into the model. Our approach is not new to the fuel cell literature, see Natarajan & Van Nguyen (2001) and Divisek *et al*. (2003), and is borrowed from the hydrology literature on hydrophobic soils, see Nieber *et al*. (2000) and references within. The liquid velocity is defined by Darcy's law,(2.13)where the liquid relative permeability *k*_{rl} is given by(2.14)where (*z*)_{+}≡Max{*z*, 0}. This constitutive form for the relative permeability supposes that for low liquid volume fractions, *β*<*β*_{*}, the water is contained in disconnected droplets which are immobile. For volume fractions larger than the percolation threshold *β*_{*}, the droplets merge to form a locally connected domain which supports transport. A second impact of hydrophobicity is on the capillary pressure, which describes the pressure jump across the pore-scale liquid–gas interfaces,(2.15)where, following Leverett (1941), we take the capillary pressure to have the form(2.16)Here, *γ*_{s} is the liquid surface tension, *S*_{p}>0 is a scaling constant which includes both contact angle effects and the substantially larger scaling for capillary pressure in hydrophobic media. For the dependence of capillary pressure on liquid volume fraction, we use a simplified van Genuchten function, see van Genuchten (1980),(2.17)which describes the increasing pressure jump required at high liquid volume fractions to force the water into increasingly smaller pores. The liquid flux, , is defined in terms of the liquid velocity, the liquid volume fraction, and the constant liquid molar density, *C*_{l}, as(2.18)

The heat flux is given by(2.19)where *c* is the specific heat, *ρ* is the density, and is the effective thermal conductivity (table 1). The effective thermal conductivity is taken to be that of the solid phase, since it is dominant,(2.20)while the heat capacity and convective heat flux are averaged over the phases(2.21)(2.22)

The condensation rate is given by(2.23)where *L*_{p}=1 μm is the typical pore size in the GDL, the coefficient of evaporation is approximated by Langmuir's formula, see Zemansky & Dittman (1981),(2.24)where *m*_{v}=0.018 kg mol^{−1} is the molecular weight of water, and the liquid–gas interface density is given by(2.25)Here, the exponent 2/3 arises from the assumption of spherical water droplets and gas bubbles, with the corresponding relationship between volume and surface area. The vapour saturation concentration, is given as(2.26)where the saturation pressure of water vapour is given from known empirical formulae, see Springer *et al*. (1991).

### (c) Boundary conditions

We prescribe fluxes at the catalyst layer–GDL boundary. At the channel–GDL boundary, we prescribe Robin type conditions to account for boundary layers introduced in the gas channel flow by the fluxes from the GDL (table 2). The fluxes out of the catalyst layer are described by three interrelated parameters: the current density, ; the heat flux out of the catalyst layer, ; and the water (liquid plus vapour) flux out of the catalyst layer, . These last two quantities are thought of as the heat and water production of the catalyst layers, which are clearly related to the current density , adding or subtracting heat and water flux between the membrane and catalyst layers. The heat and water fluxes between the membrane and catalyst layer depend upon the state of the membrane and the anode channel, as such we take and as given quantities. In the fixed boundary layer at the catalyst layer, the phase change resolves the catalyst water flux into a vapour flux and a liquid water flux(2.27)The vapour flux will be determined by the temperature profile in the vapour-saturated two-phase region. In the event that the vapour flux is toward the membrane, the liquid water flux into the GDL can exceed the water flux out of the catalyst layer, . The phase change within the catalyst–GDL boundary layer also impacts the heat flux, , into the GDL,(2.28)

We take the same mass transport parameter, , for each of the three gas species. The heat transport is primarily through the channel landing and we take a value of the heat transport parameter, , which is consistent with a solid–solid interface, see Fowler (1997). The cathode channel total gas concentration, , satisfies(2.29)and will be prescribed in place of . The value given in table 1 corresponds to a cathode pressure of three atmospheres at 353 K.

## 3. Non-dimensionalization

We scale the gas species concentrations by the channel total molar concentration. We choose the temperature scale to balance the heat flux produced by evaporation against the molar fluxes of liquid and vapour water, setting the coefficient of *Γ* in the heat equation to one. Distance is scaled by the thickness *L* of the diffusion layer, and time is scaled by *τ*_{*}=*ϵL*^{2}/*D*, the characteristic diffusive time. The dimensionless dependent and independent variables are(3.1)and the dimensionless fluxes arewhere *α* runs over the gas species oxygen, nitrogen and water vapour. The dimensionless current, or ion flux, is scaled by(3.2)

A key observation is that the current, while a large charge flux, represents a weak molar flux due to the size of Faraday's constant *F*=96 485 C mol^{−1} and the thinness of the GDL. At typical current densities , the dimensionless ion flux *I*=*O*(10^{−3}–10^{−2}).

We eliminate the nitrogen concentration *C*_{n} in favour of the dimensionless total gas concentration and rescale the conservation laws for the heat and mass transport in the GDL as(3.3)(3.4)(3.5)(3.6)(3.7)The total flux of gas species *α* is(3.8)while the gas flux is purely convective,(3.9)Here, the convective gas parameter, , is a Péclet number and the dimensionless heat capacity *σ* is given by(3.10)where liquid–solid and liquid–gas heat capacity ratios are and , where the ratio of liquid to gas molar concentrations, , has been factored out of the liquid–gas ratio.

The dimensionless liquid water flux *N*_{l} is given by(3.11)where measures the contribution of gas convection to the liquid motion and , balances capillary pressure against diffusivity. The effective water diffusivity *ν* is given by(3.12)

The dimensionless condensation rate *Γ* is given by(3.13)where the condensation rate *H* is given by(3.14)The temperature in this formula has been approximated by the dimensionless channel temperature . The dimensionless vapour saturation concentration is given by(3.15)

The dimensionless heat flux is given by(3.16)where the Lewis number , represents the balance of thermal conductivity against molecular diffusivity.

Collecting the fluxes and unknowns , and , the constitutive relations for the fluxes may be written as(3.17)with the matrix given by(3.18)where describes the impact of convective transport on the heat flux. The system of conservation laws takes the form(3.19)where describes the impact of phase change on each species, and the matrix is given by(3.20)

## 4. Scalings and reductions for the three regimes

The asymptotics for the gas, liquid and heat transport are predicated upon the small parameters, of which there are five. Chief among these is the dimensionless flux density *I*, which is the driving force for the flow. The inverse evaporation rate 1/*H* drives a singular perturbation which defines the boundary layer asymptotics. Other small parameters which are tangentially involved are the liquid–gas molar density ratio, *δ*_{l}, the Lewis number, *R*_{T}, which balances thermal conductivity against molecular diffusivity, and the inverse scaled channel temperature 1/*T*_{c} (table 3). We scale the boundary conditions at the membrane by the current density(4.1)where the scaled membrane fluxes(4.2)are written in terms of the prescribed catalyst layer fluxes , , and the undetermined GDL vapour flux . On the two-phase side of the interface, we prescribe the unknowns(4.3)where **V**_{*} is to be determined. On the dry side of the interface we prescribe the fluxes(4.4)where again **F**^{d} is to be determined. Using the projection *π*_{d} defined in §2*a*, we write the Robin conditions at the channel as(4.5)and the Dirichlet condition for the liquid water as(4.6)The channel values , are given prescribed, and the 4×4 matrix *M*_{r} is diagonal,(4.7)where and .

### (a) Dry regime

After subtracting-off the channel values, we scale the dry regime solution with the dimensionless current density *I*, after subtracting-off the channel values. In the dry region there is no liquid water present (*β*=0) and the gas is undersaturated, so there is no phase change. The scaled unknowns are , where(4.8)The relative permeabilities simplify to *k*_{rg}=1 and *k*_{rl}=0, and the conservation laws (3.3)–(3.6) and boundary conditions become(4.9)(4.10)(4.11)where ^{d} and ^{d} are 4×4 matrices. The matrix ^{d} is independent of **V**^{d}, and is given by(4.12)We expand the *z*-derivatives to find(4.13)where is given by(4.14)We neglect the *O*(*I*) terms in (4.13), and find a leading order constant-coefficient parabolic system,(4.15)where(4.16)The characteristic polynomial of ^{d} is(4.17)Using the asymptotic relations *δ*_{l}, 1/*T*_{c}≪1 we approximate the eigenvalues to find(4.18)(4.19)(4.20)The fastest time-scale is associated with the eigenvalue *λ*_{3}≈−1.7×10^{3} and has eigenvector , which describes the linearized dimensionless gas pressureThe smallest eigenvalue *λ*_{4}≈−0.25 is associated with relaxation of the temperature profile.

#### (i) The quasi-steady dry regime

In the quasi-steady regime, the two-phase interface point *z*_{*} and the dry fluxes **F**^{d} are time dependent but on a time-scale several orders of magnitude longer than the relaxation times 1/*λ*_{1}, …, 1/*λ*_{4} of the governing systems of conservation laws (4.15). Given the uniformly parabolic nature of governing equations (4.15), the corresponding solutions are quasi-stationary, adiabatically driven by the slow evolution of the front and fluxes, up to the order of the time-scale ratio, which is seen in (4.72) to be *O*(*δ*_{l}). In this regime, the fluxes are spatially constant to leading order and the governing equation reduces to(4.21)(4.22)(4.23)Since has an *O*(1) inverse, the system has the simple linear solution(4.24)

### (b) The two-phase regime

In the two-phase regime, the liquid water is everywhere greater than the immobile volume fraction, *β*_{*}, and the gas is saturated to leading order. As in the dry regime, the scaled variables describe variations from the channel values, except for the scaled water vapour and water volume fractions, which describe variations from the saturation pressure (3.15) and the immobile volume fraction respectively,(4.25)Change of variables is chosen to render the diffusivity matrix, ^{w} see (4.35), *O*(1) with an *O*(1) inverse. In particular, the (*Iδ*_{l})^{2/3} scaling for *β*^{1} is determined from the liquid flux boundary conditions at the membrane, where the effective liquid diffusivity, given by (3.12) scales as(4.26)where . Near the membrane, neglecting the convective terms, the liquid flux (3.11) has the scaling(4.27)which yields an *O*(1) value for *β*^{1} from the liquid flux boundary condition at the membrane .

In linearizing the conservation laws, it is convenient to introduce , *σ*_{*}=*σ*(*β*_{*}), and . In the new variables in the two-phase region, the liquid relative permeability simplifies,(4.28)The dimensionless vapour saturation concentration, is approximated to *O*(*I*^{2}) by linearizing about the channel temperature, *T*_{c},(4.29)The coefficient *X*_{s}=*C*_{sat}(0) is the saturation mole fraction at the channel, while gives the linearized response of the saturation mole fraction to changes in temperature, and calculated from the particular form of the dimensional saturation pressure. We collect the leading order terms in the conservation laws, finding(4.30)where(4.31)and(4.32)

Having neglected the *O*(*H*) terms associated with the vapour, the matrices and are non-invertible and the leading order equation for is degenerate. Inclusion of the *O*(*H*) terms forces a boundary-layer type solution for which is incompatible with the boundary conditions. Instead, we separate this equation out,(4.33)and eliminate from the conservation laws, obtaining a reduced 4×4 system. We drop the terms in the last row of , and replace (4.30) with a system for ,(4.34)where(4.35)and(4.36)We rewrite the conservation laws as(4.37)where(4.38)has a spread of eigenvalues from *O*(10^{3}) for pressure relaxation to *O*(10^{−2}) for the relaxation of the liquid water profile, for *β*^{1}=*O*(1).

#### (i) Quasi-steady two-phase regime

As in the dry regime, we consider the two-phase regime to be driven adiabatically by the slow variation of the two-phase front location *z*_{*}(*τ*) and the values at the two-phase interface **V**_{*}(*τ*). As a convenience, we replace *β*^{1} with the liquid ‘potential’(4.39)so that the scaled unknowns are(4.40)

With this substitution, the quasi-steady version of the conservation law (4.30) takes the form(4.41)where(4.42)is a constant matrix. The system (4.41) appears to be degenerate since det ^{wq}=0, but is well-posed. Indeed, the solvability condition(4.43)requires that , since * S* is not orthogonal to the kernel. This implies that, away from the two-phase boundary layer, the phase change is a higher-order effect arising from nonlinear interactions. The solvability condition also gives uniqueness, the only solution to(4.44)subject to is . In particular, we observe that the fluxes are constant to leading order within the two-phase region .

We drop , and introduce the reduced variables and the reduced diffusivity matrix ^{wr} which is the 4, 4 minor of ^{wq}. With a two-phase front located at *z*=*z*_{*}, the solution to (4.44) subject to (4.1) and (4.3) is a simple linear function of position(4.45)

The vapour flux carried by the saturated gas, *N*_{v}, is constant to leading order, and is given in terms of the temperature and total gas gradients in (4.33). This yields a leading order expression for the vapour flux into the GDL–catalyst boundary layer,(4.46)The analytic solution (4.45) gives the slope of the scaled molar gas concentration and temperature in terms of the prescribed fluxes,(4.47)(4.48)Substituting these expressions into (4.46) leads to a linear equation for , which has solution(4.49)The denominator is positive, while the term on the left in the numerator is negative, since the heat flux is away from the membrane ; while the total water flux is dominated by the positive term . The term on the right in the numerator is positive, since . We have two competing effects, the heat and water flux push the vapour flux negative, while the thermal conductivity, *R*_{T}, pushes it positive. The heat flux and the thermal conductivity are only partially linked, at quasi-steady state the heat produced in the MEA must leave through either the cathode or anode side of the membrane, the total heat flux is essentially prescribed. The thermal conductivity ascribes a temperature gradient necessary to accommodate the heat flux, and has only a secondary affect on the heat flux itself. For small values of *R*_{T}, as we have taken here for carbon fibre paper, the fluxes dominate and the vapour flux is negative. Indeed in the limit *R*_{T}*T*_{C}≪1, the vapour flux has the simple form(4.50)For a more thermally conductive GDL, intermediate between Toray carbon fibre paper and pure graphite, for example *κ*_{s}≈30 J (m s K)^{−1} with *R*_{T}≈1, the positive terms dominate and the vapour flux is *towards the membrane*. In the limiting case *R*_{T}*T*_{C}≫1, the vapour flux depends only upon the channel saturation mole fraction(4.51)The high thermal conductivity lowers the temperature difference between membrane and channel, forcing condensation at the GDL–catalyst boundary layer, where the liquid flux(4.52)is more negative than the total water output of the membrane. This condensation is a sink for water vapour, and the vapour flux towards the membrane gives the oxygen a convective boost, enhancing the transport of oxygen towards the membrane. More specifically, the convective boost increases the oxygen concentration at the membrane for a prescribed oxygen flux (current density) reducing oxygen reduction reaction overpotential losses. The benefit comes at the cost of increased liquid water production, and hence increased heat removal load for the coolant. Assessing the impact of the increased water production on flooding in the catalyst layer is beyond the scope of this analysis.

### (c) The two-phase interface

The free boundary layer has a double-layer structure. The first, which we resolve, is the evaporative layer in which the leading order phase change occurs. The infinitesimal layer is the point at which the liquid volume fraction reaches the immobile volume fraction, and the governing equation loses parabolicity. We address each layer in turn.

#### (i) The evaporative layer

The key to the free boundary evolution is to determine the total amount of phase change which occurs within the evaporative layer (see figure 2). We will express this net phase change and the flux jumps in terms of *Θ*_{*}, the degree of vapour undersaturation at the two-phase point ,(4.53)We assume that the liquid water is everywhere greater than the immobile volume fraction, *β*_{*}. We rescale the variables as in the two-phase region, except for the water vapour whose scaling *H*^{−α}*I* contains a parameter *α* determined below,(4.54)The local phase change takes the form(4.55)where .

We rescale the spatial variable within the evaporative layer, introducing . The natural parabolic rescaling of time *τ*_{b}=*Hτ* suggests relaxation times which are yet 5–6 magnitudes of order faster than those of the two-phase regime. For this reason we consider only the quasi-steady version of the conservation laws(4.56)

The fluxes are related to the scaled variables through the constitutive laws as(4.57)where the 5×5 matrix ^{b} is given by(4.58)

Combining the conservation laws (4.56) with the constitutive relations (4.57), we obtain a 5×5 system(4.59)We observe that(4.60)where and(4.61)The parameter *ξ*_{v} is independent of *H* and strictly positive. The equation for uncouples at leading order,(4.62)(4.63)(4.64)We assume *Θ*_{*}>0, so the constant . This equation has a classic boundary layer solution,(4.65)Returning to the conservation laws (4.59), we find that(4.66)To balance the jumps in the fluxes with the magnitude of the catalyst layer fluxes suggests the scaling(4.67)and the associated choice of *α*=1/2 in the scaling of . From (4.57), we find the leading order jump in the fluxes across the boundary layer,(4.68)

#### (ii) The sub-layer and the front motion

The free boundary has an infinitesimal sub-layer with left and right limits and , where the liquid volume fraction reaches the threshold, *β*_{*}, for motion. For piece-wise smooth solutions of the full problem (3.3)–(3.7) to form a weak solution when glued together across the interface, Rankine–Hugoniot conditions must be enforced. The non-degenerate variables *C*, *T*, *C*_{o}, and *C*_{v}, must be continuous across the interface and the corresponding flux jumps must be zero,(4.69)(4.70)The degenerate variable, *β* jumps from *β*_{*} down to 0 across the inner layer, and the liquid flux is related to the interface motion and the liquid jump through a degenerate Rankine–Hugoniot relation (Smoller 1983; Fife 1988),(4.71)The jump in liquid water flux across the evaporative layer is given by (4.68), while the flux is constant to leading order across the two-phase region . We express the front motion in terms of the water membrane flux and the undersaturation(4.72)where the vapour flux is given by (4.49). The front motion is on a time-scale *O*(10^{−5}), which is slower than the next slowest time-scale, relaxation of the *β* profile, by the factor *δ*_{l}≪1, and is 8–9 orders of magnitude slower than the fastest time-scale, gas pressure relaxation in the dry regime.

## 5. The quasi-steady sharp-interface reduction

We define a two-phase sharp-interface solution of the governing system of conservation laws (3.19) to be one which is a piecewise solution on either side of a front *z*=*z*_{*}(*τ*), which is dry (*β*=0) on one side, with a mobile two-phase regime on the other (*β*>*β*_{*}) and which satisfies the jump conditions(5.1)(5.2)and the boundary condition on the two-phase side of the interface(5.3)Here, . For such a solution, §4 shows that one can add an evaporative layer which adjusts the fluxes * N* according to (4.68) so that the remaining jump satisfies the Rankine–Hugoniot conditions (4.69)–(4.71), and forms a weak solution across the interface.

The dichotomy between relaxation times of the conservation laws and evolution of the front scale like *δ*_{l}, the density ratio between liquid water and gas. It is consistent to take the solutions in the dry and two-phase regions to be at steady state, with representing small forcing terms which can be neglected at the order we consider, except for their cumulative effects, on time-scales of which are manifest through the evolution of the free boundary and the slow change of the dry fluxes.

We construct leading-order sharp-interface solutions from the quasi-steady solutions constructed in §4*a*,*b* whose fluxes are piecewise spatially constant and given by(5.4)The jump conditions (5.2) applied to **F**^{m} and **F**^{d} yield the equation(5.5)We rewrite the dry region solution (4.24) as(5.6)Defining(5.7)(5.8)(5.9)(5.10)we obtain the expression(5.11)which yields expressions for the vapour mole fraction and dimensionless vapour saturation at the interface(5.12)We introduce the scaled channel undersaturation,(5.13)and combine the dry-region solutions (5.12) with the definition (4.53) and (4.67) of to obtain a single, linear, equation for , which has the solution(5.14)Substituting for **G**^{m}, **G**^{s}, **H**^{m}, **H**^{s} yields a leading order expression for ,(5.15)where the uniformly positive denominator *Δ* is given by,(5.16)

We summarize our result below,

*Main Result 1*. Given channel values **V**_{c}, for which the scaled channel undersaturation, , given by (5.13) is *O*(1) and given catalyst layer heat, , and water, fluxes which are *O*(*I*), we determine the catalyst layer vapour flux by (4.49) and the scaled membrane fluxes by (4.2). So long as the dry phase undersaturation , determined by (5.15), is positive, then the function(5.17)with the free boundary *z*=*z*_{*}(*τ*) given by (4.72) is, to leading-order, a sharp-interface solution of the system of conservation laws (3.19). Here, **V**_{−} is determined from (5.6) through the change of variables (4.8), and **V**_{+} is determined from (4.45) through the change of variables (4.25) and (4.39).

## 6. Results and discussion

We have presented an asymptotic reduction of the degenerate system (3.19) of coupled nonlinear partial differential equations describing multi-phase flow in a PEM fuel cell GDL. We obtain a reduced system which is comprised of a pair of linear equations on either side of a sharp interface whose evolution is given in terms of liquid flux balance. In the adiabatic limit, we obtain explicit solutions to the free interface motion of the reduced problem. This eliminates the primary source of stiffness that is the bane of numerical simulations, affording a computational speed-up of 3–4 *orders* of magnitude over the full system. Building this model into a unit cell simulation code promises huge reductions in computational cost, and admits the possibility of performing either full stack-based calculations, or doing inverse calculations and parameter estimation.

We have derived analytic expressions for the partition of product water into liquid and vapour at the GDL–catalyst layer interface and for the evolution of the free boundary. These analytical expressions provide insight into the influence of material parameters on water management and can be used in the design of fuel cell components and in the choice of materials to better optimize performance under a wide range of operating conditions.

In applications, the reduced system is embedded in a 1+1D computational scheme for the overall fuel cell. This includes a model of the membrane's water content and temperature, the anode GDL, and the variation of the oxygen and water vapour contents in the flow-field channels in the along-the-channel direction. To present numerical results from the reduced system, we simulate this coupling by providing along-the-channel data for the oxygen and water vapour concentrations, temperature, current density and catalyst layer production of heat and total water from a previous 1+1D computation reported in Bradean *et al*. (2002). These values vary in the *y* direction, but are constant in time and do not couple back to the reduced simulations.

The computations for two different runs are presented in figures 3 and 4. In figure 3, the GDL thermal conductivity is taken as *κ*_{s}=1 J (m s K)^{−1}, and the inlet and outlet ends of the channel remain dry, while the two-phase interface, plotted at 10 s intervals, grows from zero at *t*=0 to occupy the majority of the mid-channel region by *t*=80 s. The vapour flux from the catalyst layer is largest (most negative) at the dry cathode inlet and jumps as the channel flow saturates slightly downstream. In the high thermal conductivity run, figure 4, the catalyst layer total water flux, , goes positive (towards anode) at the cathode outlet *y*=0.67 m where all the product water crosses over to the anode side to hydrate the dry inlet stream. The amount of evaporation at the two-phase boundary is an order of magnitude lower than either the liquid water or vapour fluxes from the catalyst layer. In the bottom run, with *κ*_{s}=10 J (m s K)^{−1} the front velocity is an order of magnitude greater, and the two-phase region occupies all but the very end section of the cathode outlet. The two-phase front is plotted at intervals of 5 s up to *t*=15 s. The dimensionless catalyst vapour flux has dropped from roughly −7×10^{−3} to −1.6×10^{−3} while the liquid water flux has grown by approximately an order of magnitude. For the high thermal conductivity case, the evaporation at the front is large at the cathode inlet, and significantly slows the front evolution there.

The results presented have several extensions. The formal result presented in §5 can be made rigorous: the reduction of the full system to the sharp-interface solutions fits into the renormalization group approach introduced by Promislow (2002). By rescaling the spatial domain the front can be kept at a fixed position, with slowly varying coefficients introduced to the differential operators on either side of the front. The full system (3.19) can be written in the form(6.1)and the exact solution can be decomposed as(6.2)where * V* is the formal solution. The role of the renormalization group approach is to find an evolution for the free parameters, here the front position

*z*

_{0}, which guarantees that the remainder

*remains small. The requires a detailed understanding of the linear operators obtained by linearizing the flow*

**W***F*about

*(*

**V***z*;

*z*

_{0}). The application of the renormalization approach to this class of zeroth order phase transitions is a promising avenue of future work.

The model can be extended to two space dimensions by inclusion of the transverse or *x* variation (see figure 1) in the GDL which arises due to the landing area obstructing the gas, liquid and heat flow. The chief obstacle is the considerable difficulty in tracking the interface which will take the form of a curve, or family of curves, in the two-dimensional *x–y* plane. These curves can self-intersect, or have non-trivial interaction with the boundary. Adaptation to a phase-field model may be appropriate.

For applications of the reduced system to design and engineering of PEM fuel cells, the reduced system must be integrated into a 1+1D unit-cell code for the full PEM fuel cell and validated both against the full system and against experimental data. This reflects an extension of an on-going program with Ballard Power Systems, see Berg *et al*. (2004). Maintaining computational stability while coupling the reduced model to the hyperbolic flow in the gas channels and the elliptic equation for the voltage balance poses numerical challenges.

## Acknowledgements

KP thanks the NSF for their support of this research through NSF DMS 04055965. All authors thank Ballard Power Systems and Mathematics of Information Technology and Complex Systems (MITACS) for their financial support of this project.

## Footnotes

- Received August 31, 2004.
- Accepted September 12, 2005.

- © 2006 The Royal Society