Simulation of Mass and Heat Transfer in an Evaporatively Cooled PEM Fuel Cell

Abstract

Evaporative cooling is a promising concept to improve proton exchange membrane fuel cells. While the particular concept based on gas diffusion layers (GDLs) modified with hydrophilic lines (HPILs) has recently been demonstrated, there is a lack in the understanding of the mass and heat transport processes. We have developed a 3-D, non-isothermal, macro-homogeneous numerical model focusing on one interface between a HPIL and an anode gas flow channel (AGFC). In the base case model, water evaporates within a thin film adjacent to the interfaces of the HPIL with the AGFC and with the hydrophobic anode GDL. The largest part of the generated water vapor leaves the cell via the AGFC. The transport to the cathode side is shown to be partly limited by the ab-/desorption into/from the membrane. The cooling due to the latent heat has a strong effect on the local evaporation rate. An increase of the mass transfer coefficient for evaporation leads to a transport limited regime inside the MEA while the transport via the AGFC is limited by evaporation kinetics.

1. Introduction

The combination of renewable energies and hydrogen as an energy carrier is a promising alternative to fossil fuels for the energy mix of the future [1,2]. However, hydrogen conversion systems such as proton exchange membrane fuel cells (PEMFCs) still require improvements to be competitive and to ensure large-scale deployment. The management of water and heat is one of the key performance bottlenecks of PEMFCs. It needs to ensure adequate levels of membrane humidification to realize proper protonic conductivities. Typically, this is achieved by external humidifiers that load the reactant gas flows with water vapor upstream from the cells. These devices, which are bulky and expensive, are often considered as one of the sub-systems that diminish the overall efficiency of the fuel cell system and should be removed if possible. Another goal of the water management is to keep the reactant gases access to the catalyst layers (CLs) unblocked by liquid water accumulations. Liquid water inside the cell can accumulate from condensation of the vapor coming from the humidified gas streams in addition to the chemical production of water at the cathode CL at high current densities. Hence, fully humidified gases should be avoided on the cathode side. In terms of heat management, fuel cell systems typically include a water recirculation to cool the cell during high loads, which is important to avoid degradation. This, however, complicates the structures used for flow field plates.

A promising concept to optimize water and heat management is evaporative cooling [3,4,5]. It is based on the evaporation of water directly inside the cell and the simultaneous dissipation of latent heat. This scheme can operate with dry feed gases, which avoids the need for the external humidifiers mentioned above. One of the evaporation cooling concepts was developed at the Paul Scherrer Institute, which is solely based on the modification of the wettability in a gas diffusion layer (GDL) [6,7,8,9]. The wetting properties are locally turned from hydrophobic to hydrophilic within separated parallel lines, which are placed perpendicular to the gas flow channels (GFCs) (Figure 1a). These hydrophilic lines (HPILs) can wick water from a liquid water channel (LWC) within the anode flow field. The water in the HPILs then evaporates and, hence, provides simultaneous humidification and cooling.

This concept has been investigated in experiments. Cochet et al. [10] focused on the evaporation in a complete but non-operating thermal test cell (i.e., without electric current). They measured the heat flux on the anode and cathode side at different anode and cathode gas speeds. The total heat flux for gas speeds representative of a long cell behavior was higher than 1 W cm${}^{-2}$, which has the same order of magnitude as the waste heat released by an operating PEMFC. An analytical model was developed under several simplified assumptions to relate the measured heat fluxes to the unmeasured water vapor flow rates leaving the anode and cathode channels as a function of temperature and gas channel speeds. This analytical model is described in Appendix A. By changing the mass flow rate in the LWC, they varied the contact surface between water in the HPILs and the anode gas flow, which was extracted from Neutron radiography images (Figure 1a). The measured heat flux as a function of this contact surface agreed well with the prediction from the analytical model. In a subsequent study, Cochet et al. [11] extended this analysis of the non-operating cell by analyzing the relative humidity at the anode and cathode outlets as a function of the contact surface using the analytical model of [10]. Furthermore, they investigated evaporative cooling and humidification during operation. They demonstrated that the evaporative cooling cell with dry gases performed in a similar way as a PEMFC that operates in a classical fashion with pre-humidified gases.

These experiments and the analytical model provided a lot of insights into the evaporation cooling concept. However, only the heat flux was measured at the ends of the anode and cathode flow channel plate (FCP). No measurements were done or are possible for the distribution of temperature and relative humidity inside the cell, evaporation along the HPILs, water vapor transport paths and level of membrane humidification. In this regard, a numerical model with a thorough characterization of the material properties of each layer could provide additional insights and help further improve the understanding of the mass and heat transport processes in an evaporative cooling cell. The goal of this work is to develop such a model.

The macro-homogeneous modeling approach has been frequently used as a method to investigate the different coupled transport phenomena in PEMFCs (e.g., [12,13]). However, the applicability of this approach to simulate transport processes in GDLs of PEMFCs has been questioned in the past (e.g., [14,15]) due to the absence of a robust representative elementary volume (REV) considering the large pore sizes (e.g., pore diameter between 20 and 100 $\mathsf{\mu }$m for a Toray GDL) compared to the thickness (170–400 $\mathsf{\mu }$m) of uncompressed GDLs. For this reason, direct simulations are used to study evaporation inside the GDLs [16,17]. Most recently, van Rooij et al. [18] developed a direct pore-level numerical modeling framework to analyze the heat and mass transport due to evaporation. Evaporation was modeled as a mass source term using the kinetic gas theory and a $\delta$ function to represent the location of the predefined interfacial area. They investigated the effect of morphology by approximating the geometry of the GDL by an artificial lattice. Their setup was adapted to the ex-situ experiments of [16,17], and, hence, not to an evaporative cooling cell containing HPILs inside a GDL. These direct pore level simulations have the disadvantage to be computationally expensive. Therefore, they typically cover only the GDL and exclude other parts of the PEMFC such as the membrane.

For this reason, the macro-homogeneous modeling approach has continuously been adopted in simulations of PEMFCs due to the advantage of being able to efficiently couple the two-phase transport with transport of charge, heat and electro-chemical reactions (e.g., [19]). Dujc et al. [20] presented a macro-homogeneous two-phase model of a PEMFC to simulate the effects of a patterned GDL on the cathode side. The results were in good agreement with experimental ex-situ imbibition data obtained by neutron radiography imaging. Their model confirmed that a patterned GDL, applied to the cathode side of PEMFCs, promotes diffusion of oxygen from the gas channels to the catalyst sites within the hydrophobic portions by focusing the liquid water in the hydrophilic portions. This model, however, is rather simplified. First, it is isothermal and therefore cannot take the heat transport into account. Second, it includes only the cathode GDL and, hence, cannot assess the impact on the water management of the whole cell. Third, it does not account for the transport in the GFCs and the transition to the porous media domains. These aspects are crucial to realize a model of the experiments conducted by Cochet et al. [10] who measured the heat flux and used their analytical model to calculate the amount of water vapor released through the gas channels.

In this study, we present for the first time a non-isothermal 3-D macro-homogeneous model that contains all layers of the evaporative cooling test cell of Cochet et al. [10] (Section 2). In Section 3, we provide a detailed description of the water and heat transport processes induced by evaporation in a base case model. By varying the mass transport coefficient for evaporation, we analyze the role of evaporation kinetics in comparison to the transport limitations in the MEA and AGFC. Subsequently, we identify the main model limitations and propose improvements for future studies. Finally, we compare the model with the analytical model of Cochet et al. [10]. Our work sheds light on the water and heat transport paths and the humidity distribution inside a non-operating evaporatively cooled PEMFC, which serves as a basis to further improve the evaporative cooling concept and PEMFCs in general.

2. Methods

2.1. Model Setup

Figure 1b shows the numerical model setup of an evaporatively cooled PEMFC that encompasses all layers of the experimental setup of Cochet et al. [10] (Figure 1a). It focuses on a single interface between a HPIL and an AGFC. In addition to the central channel (LWC on the anode side, GFC on the cathode side), one GFC is taken into account on both the anode and cathode side. The MEA contains a central PEM sandwiched between the GDLs, MPLs, and CLs on the anode and cathode side. A vertical plane of symmetry is placed through the central flow channel. Counterflow is applied in the AGFC and CGFCs.

2.2. Basic Assumptions

The model is set up under the following assumptions for a non-operating evaporatively cooled PEMFC:

• All processes have reached a steady state;
• The porous medium of the GDL and other layers of the MEA can be described as a macro-homogeneous medium with effective transport properties;
• Single-phase flow of gas and liquid phase is assumed in the gas and liquid water flow channel, respectively. The flow is laminar under all conditions investigated, that is, the Reynold number is significantly smaller than 2000, which is considered to be the onset of sustained turbulence in pipe flow [21];
• Two-phase flow of gas and liquid phase exists in the porous medium of GDLs, MPLs and CLs. Gas and liquid water saturations are obtained by using water retention curves rather than resolving the phase interface;
• Inertial forces for the transport of gas and liquid water can be neglected;
• Gas phase behaves as an ideal gas;
• No electrochemical reactions occur;
• Local thermal equilibrium holds, that is, the temperature of the solid and fluid (and its phases and species) is the same;
• There are no thermal contact resistances;
• Gravitational forces can be neglected.

2.3. Conservation Equations

In the following, we describe the transport equations for liquid water, gas, heat, and dissolved water. The source term S, the baseline parameterization, the numerical implementation in COMSOL Multiphysics and the boundary conditions are described in Section 2.4, Section 2.5, Section 2.6 and Section 2.7.

2.3.1. Transport of Liquid Water and Gas

For the transport of liquid water and gas, we distinguish between the flow in the channels (AGFC, CGFC, LWC) and in the porous media. In the GFCs and LWC, we solve for the conservation equation of mass and momentum of a single phase (gas or liquid water), respectively, as

$$\nabla ·\left({\rho }_{\beta }{u}_{\beta }\right)=0$$

(1)

and

$$\nabla ·K_{\beta }-\nabla P_{\beta }=0$$

(2)

$P_{\beta }$, ${\rho }_{\beta }$ and $u_{\beta }$ are the pressure, density and intrinsic average velocity of the gas mixture ($\beta =$ gas) or liquid ($\beta =$ liq) phase, respectively. The medium is assumed to be viscous with the constitutive relationship defined as

$$K_{\beta }={\mu }_{\beta }\left(\nabla u_{\beta }+{(\nabla u_{\beta })}^T-\frac{2}{3}(\nabla ·u_{\beta })I\right),$$

(3)

where $K_{\beta }$ is the stress tensor, ${\mu }_{\beta }$ is the dynamic viscosity and I the identity matrix. In the porous medium of the MEA, we solve for a two-phase flow so that the entire pore space is saturated with liquid and the gas phase, satisfying the constraint

$$s_{\mathrm{gas}}+s_{\mathrm{liq}}=1,$$

(4)

where $s_{\mathrm{gas}}$ and $s_{\mathrm{liq}}$ are the saturation of the gas and liquid phase, respectively. $s_{\mathrm{liq}}$ is determined via the capillary pressure using a water retention curve, as described in Section 2.5.4. The equation for conservation of mass is defined as

$$\nabla ·\left({\rho }_{\beta }{u}_{\beta }\right)=S_{\beta },$$

(5)

where $S_{\beta }$ is the source term of the phase. The momentum conservation is described by the Brinkman equation as

$$\nabla ·K_{\beta }-\nabla P_{\beta }-u_{\beta }\left(\frac{{\mu }_{\beta }}{{\kappa }_{\mathrm{abs}}{\kappa }_{\mathrm{rel},\beta }}+\frac{S_{\beta }}{{\left({ϵ}_{\mathrm{p}}{s}_{\beta }\right)}^2}\right)=0,$$

(6)

where the porous medium is characterized by the porosity ${ϵ}_{\mathrm{p}}$, the intrinsic material permeability ${\kappa }_{\mathrm{abs}}$ and the relative permeabilities ${\kappa }_{\mathrm{rel},\beta }$. The effective constitutive relationship for the porous medium is

$$K_{\beta }=\frac{{\mu }_{\beta }}{{ϵ}_{\mathrm{p}}{s}_{\beta }}\left(\nabla u_{\beta }+{(\nabla u_{\beta })}^T-\frac{2}{3}(\nabla ·u_{\beta })I\right).$$

(7)

2.3.2. Transport of Gas Species

The gas species transport equation for gas species X is defined as

$$\nabla ·\left({\rho }_{\mathrm{gas}}{\omega }_X{u}_{\mathrm{gas}}+j_{\mathrm{gas},X}\right)=S_X,\mathrm{anode}:X={\mathrm{H}}_2,{\mathrm{H}}_2\mathrm{O},\mathrm{cathode}:X={\mathrm{H}}_2\mathrm{O},{\mathrm{N}}_2.$$

(8)

As in the experimental study of Cochet et al. [10], dry hydrogen and nitrogen enter the anode and cathode gas channel, respectively, so that we do not consider oxygen on the cathode side in this study. We solve for the mass fraction of water vapor on both sides, while the mass fractions of hydrogen and nitrogen on the anode and cathode sides, respectively, are calculated from the constraint ${\displaystyle \sum _{X=1}^N}{\omega }_X=1$. The diffusive gas flux $j_{\mathrm{gas},X}$ is defined using the Maxwell-Stefan equation following Curtiss and Bird [22] as

$$j_{\mathrm{gas},X}=-{\rho }_{\mathrm{gas}}{\omega }_X\sum _{Y=1}^N{\mathbb{D}}_{X,Y}{d}_Y,$$

(9)

where $\mathbb{D}$ is the multicomponent Fick diffusivity. d is the diffusive driving flux defined for gas species X as

$$d_X=\nabla y_X+\frac{1}{P_{\mathrm{gas}}}(y_X-{\omega }_X)\nabla P_{\mathrm{gas}},$$

(10)

where $y_X=\frac{{\omega }_X}{M_X}{M}_{\mathrm{n}}$ and $M_X$ are the molar fraction and molar mass of each species, respectively. $M_{\mathrm{n}}=\frac{1}{{\sum }_{X=1}^N\frac{{\omega }_X}{M_X}}$ is the mean molar mass.

2.3.3. Transport of Heat

The heat transport equation is governed by advective and diffusive transport defined as

$$\begin{array}{c}\hfill {\rho }_{\mathrm{f}}{C}_{\mathrm{f}}{u}_{\mathrm{f}}·\nabla T+\nabla ·j_T=S_T,j_T=-k_{\mathrm{eff}}\nabla T,\end{array}$$

(11)

where $S_T$ is the heat source term and $k_{\mathrm{eff}}$ the effective thermal conductivity. In the porous domains, the fluid density ${\rho }_{\mathrm{f}}$, heat capacity $C_{\mathrm{f}}$ and velocity $u_{\mathrm{f}}$ are calculated as an average over gas species and liquid water depending on the local saturation (Section 2.5.6).

2.3.4. Transport of Dissolved Water

Inside the ionomer, water is assumed to be in a dissolved phase [23]. The dissolved water content in the ionomer $\lambda$ is defined as the ratio of the number of water molecules to the number of charge sites in the CLs and membrane. The transport of $\lambda$ is governed by diffusion as

$$\nabla ·j_{\lambda }=S_{\lambda },j_{\lambda }=-M_{{\mathrm{H}}_2\mathrm{O}}\frac{{\mathcal{D}}_{\lambda }}{V_{\mathrm{m},\mathrm{i}}}\nabla \lambda ,$$

(12)

where ${\mathcal{D}}_{\lambda }$ is the effective water diffusivity and $S_{\lambda }$ is the source term. $V_{\mathrm{m},\mathrm{i},\mathrm{dry}}=\frac{\mathrm{EW}}{{\rho }_{\mathrm{i},\mathrm{dry}}}$ is the equivalent molar volume of the ionomer, which describes the ratio of polymer volume per moles of sulfonic acid groups. EW is the equivalent weight and ${\rho }_{\mathrm{i},\mathrm{dry}}$ is the mass density of the completely dry ionomer, respectively (

Table A3. Dry ionomer properties (for Nafion-NR-211 type).

Symbol	Explanation	Value
EW	equivalent weight	1.02 kg mol${}^{-1}$ [42]
${\rho }_{\mathrm{i},\mathrm{dry}}$	mass density	1970 kg m${}^{-3}$ [43]

).

2.4. Source Terms

The conservation equations are coupled via the source terms related to evaporation/condensation and ab-/desorption in the MEA summarized in

Table 1. Source terms.

Symbol	AGDL	AMPL	ACL	PEM	CCL	CMPL	CGDL
$S_{\mathrm{liq}}$	$-S_{\mathrm{ec}}$	$-S_{\mathrm{ec}}$	$-S_{\mathrm{ec}}$	-	$-S_{\mathrm{ec}}$	$-S_{\mathrm{ec}}$	$-S_{\mathrm{ec}}$
$S_{\mathrm{gas}}=S_{{\mathrm{H}}_2\mathrm{O}}$	$S_{\mathrm{ec}}$	$S_{\mathrm{ec}}$	$S_{\mathrm{ec}}-S_{\mathrm{ad}}$	-	$S_{\mathrm{ec}}-S_{\mathrm{ad}}$	$S_{\mathrm{ec}}$	$S_{\mathrm{ec}}$
$S_T$	$S_{T,\mathrm{ec}}$	$S_{T,\mathrm{ec}}$	$S_{T,\mathrm{ec}}+S_{T,\mathrm{ad}}$	-	$S_{T,\mathrm{ec}}+S_{T,\mathrm{ad}}$	$S_{T,\mathrm{ec}}$	$S_{T,\mathrm{ec}}$
$S_{\lambda }$	-	-	$S_{\mathrm{ad}}$	-	$S_{\mathrm{ad}}$	-	-

.

Source terms for liquid water, water vapor and dissolved water include the phase transformation evaporation/condensation $S_{\mathrm{ec}}$ (Section 2.5.5) and membrane ab-/desorption $S_{\mathrm{ad}}$ (Section 2.5.8). The heat source term $S_T$ includes latent heat of evaporation/condensation $S_{\mathrm{T},\mathrm{ec}}=-H_{\mathrm{ec}}{S}_{\mathrm{ec}}$ and latent heat of ab-/desorption $S_{\mathrm{T},\mathrm{ad}}=H_{\mathrm{ad}}{S}_{\mathrm{ad}}$, with $H_{\mathrm{ec}}$ (Equation (31)) and $H_{\mathrm{ad}}$ (Equations (53) and (54)) being the respective molar enthalpies.

2.5. Base Case Operating Conditions and Parameterization

The properties of the different model domains are given in

Table A1. Properties of the model domains.

Property	Unit	FCP	AGDL	AMPL	ACL	PEM	CCL	CMPL	CGDL
Thickness L	$\mathsf{\mu }$m	1000	150	65	4	10	4	45	150
Porosity ${ϵ}_{\mathrm{p}}$	1	1	0.72	0.5	0.4	–	0.4	0.5	0.72
Pore tortuosity ${\tau }_{\mathrm{p}}$	–	–	2.89	1.59	2.25	–	2.25	1.59	2.89
Pore surface area density $a_{\mathrm{p}}$	m${}^2$/cm${}^2$	–	27.6	27.6	27.6	210	27.6	27.6	27.6
Ionomer tortuosity ${\tau }_{\mathrm{i}}$	–	–	–	–	1.4	1	1.4	–	–
Ionomer volume fraction ${ϵ}_{\mathrm{i}}$	–	–	–	–	0.3	1	0.3	–	–
Absolute permeability ${\kappa }_{\mathrm{abs}}$	$\mathsf{\mu }$m${}^2$	–	2.14	0.041	0.1	–	0.1	0.041	2.14
Through-plane dry thermal conductivity $k_{\mathrm{dry}}$	W/(mK)	16.3	0.64 ${}^{\ast 1}$	0.3	22	0.11 ${}^{\ast 2}$	22	0.3	0.64 ${}^{\ast 1}$

FCP: Stainless steel 316L. AGDL: Modified Toray TGP-H-060 carbon paper with HPILs. AMPL: Carbel MP30Z (W. L. Gore and Associates) with additional microporous layer. ACL, PEM and CCL: catalyst coated membrane Primea 5710 (W.L. Gore & Associates, Inc., Newark, DE, USA). CMPL and CGDL: SGL 24BC, with an integrated MPL containing 23% in weight of PTFE. ${*}^1$: evaluated for T = 80 °C and P_cl= 6 MPa in Equation Equation (39). ${*}^2$: evaluated for T = 80 °C in Equation Equation (34).

.

2.5.1. Operating Conditions

Table 2. Base case operating conditions.

Symbol	Explanation	Value	Unit
$T^{\mathrm{B}}$	Temperature at anode and cathode FCP ends	80	°C
$P_{\mathrm{gas}}^{\mathrm{B}}$	Gas pressure at outlets of GFCs	2.0	bar
RH${}^{\mathrm{AGFC},\mathrm{in}}$	Relative humidity at inlet of AGFC	0	–
RH${}^{\mathrm{CGFC},\mathrm{in}}$	Relative humidity at inlet of CGFC	0	–
${\dot{V}}_{\mathrm{gas}}^{\mathrm{AGFC},\mathrm{in}}$	Volumetric flow rate at AGFC inlet	24	mL min${}^{-1}$
${\dot{V}}_{\mathrm{gas}}^{\mathrm{CGFC},\mathrm{in}}$	Volumetric flow rate at CGFC inlet	96	mL min${}^{-1}$
$P_{\mathrm{c}}^{\mathrm{B}}$	Capillary pressure at outlet of LWC	0.5	mbar
${\dot{m}}_{\mathrm{liq}}^{\mathrm{B}}$	Mass flow rate at LWC inlet	8	g h${}^{-1}$

lists the operating conditions for the baseline parameterization, which determine the boundary and initial conditions discussed in Section 2.7.

2.5.2. Liquid Water Properties

The water density ${\rho }_{\mathrm{liq}}$ is given at standard atmospheric pressure as [24]

$$\begin{array}{c}\hfill {\rho }_{\mathrm{liq}}={\rho }_{\mathrm{crit}}\left(1+1.9927{\widehat{T}}^{1/3}+1.0997{\widehat{T}}^{2/3}-0.51084{\widehat{T}}^{5/3}-1.7549{\widehat{T}}^{16/3}-45.517{\widehat{T}}^{43/3}-6.74694{10}^5{\widehat{T}}^{110/3}\right),\end{array}$$

(13)

where $\widehat{T}=1-T/T_{\mathrm{crit}}$ and $T_{\mathrm{crit}}$ and ${\rho }_{\mathrm{crit}}$ are the critical temperature and density, respectively (

Table A2. Liquid water properties.

Symbol	Explanation	Value
${\rho }_{\mathrm{crit}}$	Critical mass density of water	$322\mathrm{kg}/{\mathrm{m}}^3$
$T_{\mathrm{crit}}$	Critical temperature of water	647.1 K
$P_{\mathrm{crit}}$	Critical pressure of water	22.064 MPa

). The dynamic viscosity of the liquid water phase is given as [25]

$${\mu }_{\mathrm{liq}}=\sum _{i=1}^4{c}_i{\left(\frac{T}{300\left[\mathrm{K}\right]}\right)}^{b_i}\left[\mathsf{\mu }\mathrm{Pa}\mathrm{s}\right],$$

(14)

where coefficients $c_i$ and $b_i$ are given in Ref. [25]. The thermal conductivity is following the internationally recommended correlation at $1\mathrm{bar}$ up to $110 °\mathrm{C}$ [26] as

$$k_{\mathrm{liq}}=\sum _{i=1}^4{c}_i{\left(\frac{T}{300\left[\mathrm{K}\right]}\right)}^{b_i}\left[\frac{\mathrm{W}}{\mathrm{m}\mathrm{K}}\right],$$

(15)

where coefficients $c_i$ and $b_i$ are given in Ref. [26]. The specific heat capacity of liquid water is set to the value provided in reference [27] for $T=350$ [K] and $P=1$ [bar].

$$C_{\mathrm{liq}}=4194.5\left[\frac{\mathrm{J}}{\mathrm{kgK}}\right].$$

(16)

2.5.3. Gas Properties

The gas density is defined assuming the ideal gas law as

$${\rho }_{\mathrm{gas}}=M_{\mathrm{n}}\frac{P_{\mathrm{gas}}}{RT}.$$

(17)

The gas viscosity ${\mu }_{\mathrm{gas}}$ is calculated as an weighted average from the viscosities of the individual gas species as

$${\mu }_{\mathrm{gas}}=\sum _{X=1}^N{y}_X{\mu }_X,$$

(18)

where the viscosities of the gas species ${\mu }_X$ is calculated as

$${\mu }_X=\sum _{i=0}^6{b}_i{\left(\frac{T}{1000\left[\mathrm{K}\right]}\right)}^i\left[\mu \mathrm{P}\right]$$

(19)

with coefficients $b_i$ given in Table 2 in reference [28].

The multicomponent Fick diffusivities ${\mathbb{D}}_{X,Y}$ are defined via the effective binary Maxwell-Stefan diffusivity ${\mathcal{D}}_{X,Y}$ following [22] as

$${\mathcal{D}}_{X,Y}=-\frac{y_X{y}_Y}{{\omega }_X{\omega }_Y}\frac{{\displaystyle \sum _{\gamma \ne X}}{\mathbb{D}}_{X\gamma }{\left(\mathrm{adj}{B}_X\right)}_{\gamma Y}}{{\displaystyle \sum _{\gamma \ne X}}{\left(\mathrm{adj}{B}_X\right)}_{\gamma Y}},$$

(20)

where ${\left(B_X\right)}_{\gamma Y}=-{\mathbb{D}}_{Y\gamma }+{\mathbb{D}}_{X\gamma }$. $\mathbb{D}$ is symmetric and obeys the relationship ${\sum }_X{\omega }_X{\mathbb{D}}_{X,Y}=0$. See Appendix B.5.2 for explicit representation of ${\mathbb{D}}_{X,Y}$ for the anode and cathode gas species. The effective binary gas diffusivity ${\mathcal{D}}_{X,Y}$ are computed as in [19,29] following the Chapman–Enskog kinetic gas theory for non-polar gases (and assuming that it can also be applied for water vapor) as

$${\mathcal{D}}_{\mathrm{X}},\mathrm{Y}={\zeta }_{\mathrm{p}}\frac{3}{8}\sqrt{\frac{RT}{2\pi }\left(\frac{1}{M_X}+\frac{1}{M_Y}\right)}\frac{KT}{P_{\mathrm{gas}}{\left(\frac{{\sigma }_X+{\sigma }_Y}{2}\right)}^2{\Omega }_{X,Y}},$$

(21)

where ${\zeta }_{\mathrm{p}}$ is the microstructure factor of the porous medium, K is the Boltzmann constant and ${\sigma }_X$ is the Lennard–Jones collision diameter (

Table A4. Gas properties: values marked by *${}^1$ are from Ref. [32].

Symbol	Explanation	Unit	H2	N2	H2O
M	Molar mass	g mol${}^{-1}$	2.016 *1	28.013 *1	18.015
${\sigma }_X$	Lennard-Jones collision diameter	Å	2.915 *1	3.667 *1	3.166
${\epsilon }_X$	Lennard-Jones potential depth (divided by Boltzmann constant)	K	38.0 *1	99.8 *1	78.2

). The Lennard–Jones collision integral is defined as [30,31,32]

$${\Omega }_{X,Y}=\frac{1.06036}{{\left(T^{\ast }\right)}^{0.15610}}+\frac{0.19300}{\exp \left(0.47635T^{\ast }\right)}+\frac{1.03587}{\exp \left(1.52996T^{\ast }\right)}+\frac{1.76474}{\exp \left(3.89411T^{\ast }\right)},$$

(22)

where the dimensionless temperature $T^{\ast }$ is related to the Lennard-Jones potential depth ${\epsilon }_X$ (Table A4) by $T^{\ast }=\frac{T}{\sqrt{{\epsilon }_X·{\epsilon }_Y}}$.

2.5.4. Porous Media Properties, Wettability and Capillary Pressure-Saturation Relationship

$s_{\mathrm{liq}}$ is determined as a function of the capillary pressure $P_{\mathrm{c}}$ as

$$s_{\mathrm{liq}}=f\left(P_{\mathrm{c}}\right),$$

(23)

where the capillary pressure $P_{\mathrm{c}}$ is defined as [8,20]

$$P_{\mathrm{c}}=P_{\mathrm{liq}}-P_{\mathrm{gas}}.$$

(24)

The capillary pressure-liquid water saturation curve f is obtained for the hydrophobic anode GDL and HPIL by applying a piecewise cubic fit to data given in [8] (Figure A1).

We define the immobile liquid water saturation as $s_{\mathrm{liq},\mathrm{im}}=0.057$. We assume that the wettability of the other domains of the MEA is the same as of the hydrophobic part of the AGDL. Note that under the operating conditions listed in Table 2, this results in $s_{\mathrm{liq},\mathrm{red}}=0$ everywhere except of the HPIL. The resulting $s_{\mathrm{liq}}$ determines the microstructure factor ${\zeta }_{\mathrm{p}}$ as

$${\zeta }_{\mathrm{p}}=\frac{{ϵ}_{\mathrm{p}}}{{\tau }_{\mathrm{p}}}{(1-s_{\mathrm{liq}})}^3$$

(25)

and the relative permeabilities ${\kappa }_{\mathrm{rel},\mathrm{wp}}$, which are defined via the saturation of the wetting (wp) and non-wetting (nwp) phase, respectively, as

$${\kappa }_{\mathrm{rel},\mathrm{wp}}=s_{\mathrm{wp}}^3,{\kappa }_{\mathrm{rel},\mathrm{nwp}}=s_{\mathrm{nwp}}^3$$

(26)

with an exponent of 3 following through-plane experimental data of hydrophobic GDLs [33]. In the hydrophilic material, the wetting and non-wetting phase are the liquid and gas phase, respectively. In the hydrophobic material, the wetting and non-wetting phase are the gas and liquid phase, respectively.

2.5.5. Evaporation

Evaporation and condensation is modeled as a volumetric source term $S_{\mathrm{ec}}$ for the gas and liquid phase using continuum expression of the Hertz-Knudsen equation [19,20,34] as

$$S_{\mathrm{ec}}={\alpha }_{\mathrm{ec}}\left({\rho }_{\mathrm{gas},\mathrm{sat}}-{\rho }_{\mathrm{gas},{\mathrm{H}}_2\mathrm{O}}\right)a_{\mathrm{p}}\{\begin{array}{cc}{s}_{\mathrm{liq},\mathrm{red}}\hfill & {\rho }_{\mathrm{gas},{\mathrm{H}}_2\mathrm{O}}<{\rho }_{\mathrm{gas},\mathrm{sat}}\left(\mathrm{evaporation}\right),\hfill \\ (1-s_{\mathrm{liq},\mathrm{red}})\hfill & {\rho }_{\mathrm{gas},{\mathrm{H}}_2\mathrm{O}}>{\rho }_{\mathrm{gas},\mathrm{sat}}\left(\mathrm{condensation}\right),\hfill \end{array}$$

(27)

where ${\rho }_{\mathrm{gas},{\mathrm{H}}_2\mathrm{O}}$ and ${\rho }_{\mathrm{gas},\mathrm{sat}}$ are the densities of water vapor at $P_{\mathrm{gas}}$ and at saturation pressure $P_{\mathrm{sat}}$, respectively. $a_{\mathrm{p}}$ is the pore area density. $a_{\mathrm{p}}{s}_{\mathrm{liq},\mathrm{red}}$ can be considered as the up-scaled liquid water interfacial area density $a_{\mathrm{liq},\mathrm{gas}}$, where $s_{\mathrm{liq},\mathrm{red}}$ is the reduced liquid water saturation, which is defined as

$$s_{\mathrm{liq},\mathrm{red}}=\frac{s_{\mathrm{liq}}-s_{\mathrm{liq},\mathrm{im}}}{1-s_{\mathrm{liq},\mathrm{im}}}.$$

(28)

${\alpha }_{\mathrm{ec}}$ is the mass transfer coefficient for evaporation/condensation that is defined as

$${\alpha }_{\mathrm{ec}}=\chi {\mathrm{\Gamma }}_{\mathrm{s}}{\mathrm{\Gamma }}_{\mathrm{m}}\sqrt{\frac{RT}{2\pi M_{{\mathrm{H}}_2\mathrm{O}}}},$$

(29)

where ${\mathrm{\Gamma }}_{\mathrm{s}}=0.1$ is an interfacial accommodation factor, ${\mathrm{\Gamma }}_{\mathrm{m}}=5\times {10}^{-4}$ is an uptake coefficient that accounts for the combined effects of heat and mass transport in the vicinity of the liquid-vapor interface. $\chi$ is a factor that is used in this study to investigate the role of the combined role of ${\mathrm{\Gamma }}_{\mathrm{s}}{\mathrm{\Gamma }}_{\mathrm{m}}$ in defining ${\alpha }_{\mathrm{ec}}$. We use $\chi =10$ in base case model presented in Section 3.1 to approach a fully saturated gas at the borders of the HPIL. The choice of $\chi$ and the impact on the results is discussed in Section 3.2.

$P_{\mathrm{sat}}$ is defined as [24]

$$P_{\mathrm{sat}}=P_{\mathrm{crit}}exp\left[\frac{T_{\mathrm{crit}}}{T}\left(-7.8595\widehat{T}+1.8441{\widehat{T}}^{1.5}-11.787{\widehat{T}}^3+22.681{\widehat{T}}^{3.5}-15.962{\widehat{T}}^4+1.8012{\widehat{T}}^{7.5}\right)\right],$$

(30)

where $P_{\mathrm{crit}}$ is the critical pressure of water. The temperature dependence of latent heat for evaporation/condensation is parameterized as [19]

$$H_{\mathrm{ec}}=2914.79\left[\mathrm{kJ}/\mathrm{kg}\right]\times \exp \left(0.261\ln \left(\widehat{T}\right)-0.044\ln {\left(\widehat{T}\right)}^2-0.0044\ln {\left(\widehat{T}\right)}^3\right).$$

(31)

2.5.6. Effective Thermal Properties

We use the parameterization of [19] to define the thermal conductivity of the humidified membranes as

$$k_{\mathrm{eff}}^{\mathrm{PEM}}={ϵ}_{\mathrm{i},\mathrm{liq}}{k}_{\mathrm{liq}}+\left(1-{ϵ}_{\mathrm{i},\mathrm{liq}}\right)k_{\mathrm{dry}}^{\mathrm{PEM}}.$$

(32)

${ϵ}_{\mathrm{i},\mathrm{liq}}$ denotes the volume fraction of water in the ionomer

$${ϵ}_{\mathrm{i},\mathrm{liq}}=\frac{\lambda V_{\mathrm{m}}}{\lambda V_{\mathrm{m}}+V_{\mathrm{i}}},$$

(33)

where $V_{\mathrm{m}}=\frac{M_{{\mathrm{H}}_2\mathrm{O}}}{\rho \mathrm{liq}}$ is the molar volume of liquid water. $k_{\mathrm{dry}}^{\mathrm{PEM}}$ [19,35] is the dry thermal conductivity of the ionomer

$$k_{\mathrm{dry}}^{\mathrm{PEM}}=\left(0.451-0.286\frac{T}{300\left[\mathrm{K}\right]}\right)\frac{\mathrm{W}}{\mathrm{m}\mathrm{K}}.$$

(34)

The gas phase conductivity $k_{\mathrm{gas}}$ and specific heat capacity $C_{\mathrm{gas}}$ depend on the gas composition. We model it as a linear combination of the conductivities of the individual gas components, with the mole fraction as weights and species conductivities $k_X$ from Green and Perry [36]:

$$k_{\mathrm{gas}}=\sum _X{y}_X{k}_X$$

(35)

and

$$C_{\mathrm{gas}}=\sum _X{y}_X{C}_{\mathrm{X}}.$$

(36)

In the GFCs, $k_{\mathrm{eff}}^{\mathrm{AGFC},\mathrm{CGFC}}=k_{\mathrm{gas}}$, ${\rho }_{\mathrm{f}}^{\mathrm{AGFC},\mathrm{CGFC}}={\rho }_{\mathrm{gas}}$, $u_{\mathrm{f}}^{\mathrm{AGFC},\mathrm{CGFC}}=u_{\mathrm{gas}}$ and $C_{\mathrm{f}}^{\mathrm{AGFC},\mathrm{CGFC}}=C_{\mathrm{gas}}$.

In the LWC, $k_{\mathrm{eff}}^{\mathrm{LWC}}=k_{\mathrm{liq}}$, ${\rho }_{\mathrm{f}}^{\mathrm{LWC}}={\rho }_{\mathrm{liq}}$, ${u_{\mathrm{f}}}^{\mathrm{LWC}}=u_{\mathrm{liq}}$ and $C_{\mathrm{f}}^{\mathrm{LWC}}=C_{\mathrm{liq}}$. For the porous media, the effective thermal conductivity is

$$k_{\mathrm{eff}}^{\mathrm{GDL},\mathrm{MPL},\mathrm{CL}}=(1-{ϵ}_{\mathrm{p}})k_{\mathrm{solid}}+{ϵ}_{\mathrm{p}}{k}_{\mathrm{f}}=k_{\mathrm{dry}}+{ϵ}_{\mathrm{p}}{k}_{\mathrm{f}}.$$

(37)

The theoretical thermal conductivity of the solid bulk material is defined as

$$k_{\mathrm{solid}}=\frac{k_{\mathrm{dry}}}{1-{ϵ}_{\mathrm{p}}},$$

(38)

where $k_{\mathrm{dry}}$ is the effectively measured thermal conductivity of the dry porous layer. For GDLs, Vetter and Schumacher [19] fitted the following function to experimental data of SGL 10 series [37]:

$$k_{\mathrm{dry}}^{\mathrm{GDL}}=\left(0.776-0.430\frac{T}{300\left[\mathrm{K}\right]}\right){\left(\frac{P_{\mathrm{cl}}}{1\left[\mathrm{bar}\right]}\right)}^{0.21}\frac{\mathrm{W}}{\mathrm{m}\mathrm{K}},$$

(39)

where $P_{\mathrm{cl}}$ is the clamping pressure. For the CLs and MPLs, a constant value is used for the dry thermal conductivity (Table A1) and humidity dependence of the effective thermal conductivity is added through Equations (37) and (40) in the CLs just like in the GDLs. For the FCPs, $k_{\mathrm{eff}}^{\mathrm{FCP}}=k_{\mathrm{dry}}^{\mathrm{FCP}}$ (Table A1) and $u_{\mathrm{f}}=0$.

$k_{\mathrm{f}}$ the conductivity of the fluid filling the pore space. For the fluid conductivity $k_{\mathrm{f}}$, specific heat capacity $C_{\mathrm{f}}$ and density ${\rho }_{\mathrm{f}}$, we assume that liquid water and the gas mixture form transport channels in through-plane direction along which heat is transported in parallel:

$$k_{\mathrm{f}}=s_{\mathrm{liq}}{k}_{\mathrm{liq}}+(1-s_{\mathrm{liq}})k_{\mathrm{gas}}.$$

(40)

$$C_{\mathrm{f}}=s_{\mathrm{liq}}{C}_{\mathrm{liq}}+(1-s_{\mathrm{liq}})C_{\mathrm{gas}}.$$

(41)

$${\rho }_{\mathrm{f}}=s_{\mathrm{liq}}{\rho }_{\mathrm{liq}}+(1-s_{\mathrm{liq}}){\rho }_{\mathrm{gas}}.$$

(42)

The fluid velocity is defined as

$$u_{\mathrm{f}}=\frac{s_{\mathrm{liq}}{\rho }_{\mathrm{liq}}{u}_{liq}+(1-s_{\mathrm{liq}}){\rho }_{\mathrm{gas}}{u}_{\mathrm{gas}}}{{\rho }_{\mathrm{f}}}.$$

(43)

2.5.7. Dissolved Water Diffusivity

The temperature dependent dissolved water diffusivity in the ionomer is defined as [13]

$${\mathcal{D}}_{\lambda }={\zeta }_{\mathrm{i}}{D}_{\mathrm{F}}\exp \left[\frac{E_{\mathrm{d}}}{R}\left(\frac{1}{353.15\mathrm{K}}-\frac{1}{T}\right)\right].$$

(44)

The parameterization of the Fickian diffusivity ${\mathcal{D}}_{\mathrm{F}}$ as a function of $\lambda$ is a subject of controversy [19]. Especially the existence of a local maximum is unclear. Here we use the parameterization of Vetter and Schumacher [13] with a peak at $\lambda \approx 4$:

$${\mathcal{D}}_{\mathrm{F}}=\frac{3.842{\lambda }^3-32.03{\lambda }^2+67.74\lambda }{{\lambda }^3-2.115{\lambda }^2-33.013\lambda +103.3}{10}^{-6}{\mathrm{cm}}^2/\mathrm{s}.$$

(45)

The microstructure factor of the ionomer phase ${\zeta }_{\mathrm{i}}$ is given by Vetter and Schumacher [13] following the idea of the microstructure factor ${\zeta }_{\mathrm{p}}$ as

$${\zeta }_{\mathrm{i}}={ϵ}_{\mathrm{i}}/{\tau }_{\mathrm{i}}^2,$$

(46)

where ${ϵ}_{\mathrm{i}}$ and ${\tau }_{\mathrm{i}}$ are the ionomer volume fraction and ionomer tortuosity, respectively. The activation energy $E_{\mathrm{d}}$ is defined as

$$E_{\mathrm{d}}=\left(38.0{ϵ}_{\mathrm{i},\mathrm{liq}}^2-47.9{ϵ}_{\mathrm{i},\mathrm{liq}}+29.2\right)\frac{\mathrm{kJ}}{\mathrm{mol}}.$$

(47)

2.5.8. Absorption and Desorption

The absorption and desorption of water vapor into and from the ionomer phase of the CL’s is simulated via the sorption source term $S_{\mathrm{ad}}$ defined as

$$S_{\mathrm{ad}}=M_{{\mathrm{H}}_2\mathrm{O}}\frac{{\lambda }_{\mathrm{eq}}-\lambda }{V_{\mathrm{i}}}\times \{\begin{array}{ll}\frac{{\alpha }_{\mathrm{a}}}{L^{\mathrm{ACL}}}\hfill & \mathrm{if} \lambda <{\lambda }_{\mathrm{eq}}\left(\mathrm{absorption}\right)\\ \hfill \frac{{\alpha }_{\mathrm{d}}}{L^{\mathrm{CCL}}}& \mathrm{if} \lambda >{\lambda }_{\mathrm{eq}}\left(\mathrm{desorption}\right),\end{array}$$

(48)

where ${\alpha }_{\mathrm{a},\mathrm{d}}$ are the ab-/desorption mass transfer coefficients defined as [19,38]

$${\alpha }_{\mathrm{a},\mathrm{d}}={ϵ}_{\mathrm{i},\mathrm{liq}}\exp \left[\frac{20[\mathrm{kJ}/\mathrm{mol}]}{R}\left(\frac{1}{303.15\left[\mathrm{K}\right]}-\frac{1}{T}\right)\right]\times \{\begin{array}{cc}\hfill 1.14\times {10}^{-5}\left[\mathrm{m}{\mathrm{s}}^{-1}\right]& \mathrm{if} \lambda <{\lambda }_{\mathrm{eq}}\left(\mathrm{absorption}\right)\\ \hfill 4.59\times {10}^{-5}\left[\mathrm{m}{\mathrm{s}}^{-1}\right]& \mathrm{if} \lambda >{\lambda }_{\mathrm{eq}}\left(\mathrm{desorption}\right)\end{array}$$

(49)

As in previous studies, we take into account that the membrane is at the same time in contact with liquid water and water vapor. Furthermore, Schroeder’s paradox is taken into account. Thus, the equilibrium hydration number ${\lambda }_{\mathrm{eq}}$ of the membrane is defined as

$${\lambda }_{\mathrm{eq}}=s_{\mathrm{liq}}{\lambda }_{\mathrm{eq},\mathrm{liq}}+s_{\mathrm{gas}}{\lambda }_{\mathrm{eq},\mathrm{gas}},$$

(50)

where ${\lambda }_{\mathrm{eq},\mathrm{liq}}$ and ${\lambda }_{\mathrm{eq},\mathrm{gas}}$ denote the hydration number when the membrane is liquid-equilibrated and vapor-equilibrated, respectively. Following Springer et al. [23],

$${\lambda }_{\mathrm{eq},\mathrm{gas}}=0.043+17.81\times \mathrm{RH}-39.85\times {\mathrm{RH}}^2+36.0\times {\mathrm{RH}}^3$$

(51)

and

$${\lambda }_{\mathrm{eq},\mathrm{liq}}=22.$$

(52)

The molar enthalpy of absorption/desorption is defined as

$$H_{\mathrm{ad}}=H_{\mathrm{ec}}+H_{\mathrm{mix}},$$

(53)

where Ref. [19] parameterized the mixing enthalpy as a function of temperature and membrane hydration as

$$H_{\mathrm{mix}}=\frac{1}{M_{{\mathrm{H}}_2\mathrm{O}}}\left(c_1\exp (-b_1\lambda )+c_2\lambda \exp (-b_2{\lambda }^2)\right),$$

(54)

where fit parameters $c_{1,2}$ and $b_{1,2}$ given in Ref. [19].

2.6. Numerical Implementation

We implemented our model in COMSOL Multiphysics Version 5.6 using the interfaces “Free and porous media flow”, “Heat transfer in porous media”, “Transport of Concentrated Species”, and “General form PDE”. The modeling domain is discretized using finite elements (Figure 2a) with linear shape functions. We have put special focus on resolving the interfaces of the HPIL with the AGFC and AGDL, going down to element sizes of less than 0.01 $\mathsf{\mu }$m (Figure 2e). The reason is that an evaporation film develops over which large gradients occur. At the same time we aim to reduce the computational costs by separating the dense mesh surrounding the HPIL to a coarser mesh in the rest of the model. The detailed meshing sequence is as follows. First, a mapped, quadrilateral mesh is applied to the interface between AGDL and AMPL. This mapped face is swept in x-direction to the CMPL/CGDL interface with a resolution as shown in Figure 2b. Next the AFCP end is meshed with quadrilateral elements. The resolution is chosen such that it resolves the transition from hydrophilic HPIL and hydrophobic AGDL as well from free flow in GFC and porous media flow in AGDL (Figure 2c–e). This surface mesh is then swept in x-direction to the bottom surface of the HPIL. The CGDL/CGFC interface is mapped as shown in Figure 2a. The bottom part of the AGDL and the full CGDL, CFCP and CGFC domains are meshed with free tetrahedra. Elements are anisotropic with a scale factor of 2 in x-direction. Finally, boundary layers are added to the top, bottom and side boundaries except of the symmetry plane. The transition from the dense mesh to the coarser mesh is not straightforward and causes many small elements in the bottom part of the AGDL. Mesh convergence is discussed in Appendix C. We use the Newton method to solve simultaneously for the fully coupled system of non-linear differential equations. Newton iterations are conducted until the solution or the residual reaches a relative error of ${10}^{-4}$. PARDISO (parallel sparse direct solver) is used as linear system solver. To reach convergence, the $\chi$ factor in Equation (27) is gradually ramped up from ${10}^{-4}$ to ${10}^2$. The solution time for calculating the base case model with 257,707 elements is around 2 days and 14 h with a memory usage of 130 to 140 GB on 2 sockets with 15 cores (Intel Xeon CPU E5-2680 v2 at 2.80 GHz).

2.7. Boundary and Initial Conditions

In the following, we describe the boundary conditions, and the initial conditions for the first iteration of the non-linear solver.

2.7.1. Transport of Liquid Water and Gas

We prescribe a fully developed flow boundary condition at the inlets and outlets of the LWC and GFCs (Figure 1b), that is, the flow has a developed velocity profile before it enters into the model domain. At the inlets, it is required that the tangential flow component is zero and that the average velocity is equal to the phase flow rate $\dot{V_{\beta }}$:

$$u_{\beta }-(u_{\beta }·n)n=0,\underset{A^{\mathrm{GFC},\mathrm{in}}}{\int \int }{u}_{\beta }·n \mathrm{d}x \mathrm{d}z=\dot{V_{\beta }},$$

(55)

where $A^{\mathrm{GFC},\mathrm{in}}=D^{\mathrm{GFC}}{L}^{\mathrm{GFC}}$ is the area of the inlet of the LWC or GFC and n is the unit interface normal vector. The phase flow rate $\dot{V_{\beta }}$ is defined as ${\dot{V}}_{\mathrm{liq}}=\frac{{\dot{m}}_{\mathrm{liq}}^{\mathrm{B}}}{{\rho }_{\mathrm{liq}}}$ and ${\dot{V}}_{\mathrm{gas}}={\dot{V}}_{\mathrm{gas}}^{\mathrm{AGFC}/\mathrm{CGFC},\mathrm{in}}$, respectively. At the outlets of the LWC and GFCs, it is required that

$$u_{\beta }-(u_{\beta }·n)n=0,\frac{1}{A^{\mathrm{GFC},\mathrm{in}}}\underset{A^{\mathrm{GFC},\mathrm{in}}}{\int \int }{P}_{\beta } \mathrm{d}x \mathrm{d}z={\overline{P}}_{\beta },$$

(56)

where the outlet average pressure is ${\overline{P}}_{\mathrm{liq}}=P_{\mathrm{gas}}^{\mathrm{B}}+P_{\mathrm{c}}^{\mathrm{B}}$ and ${\overline{P}}_{\mathrm{gas}}=P_{\mathrm{gas}}^{\mathrm{B}}$, respectively. For a detailed description of the fully developed flow boundary condition in COMSOL Multiphysics the reader is referred to the COMSOL documentation [39,40].

Symmetry conditions are applied along the symmetry plane of the model (see Figure 1b), which prescribes no penetration and vanishing shear stresses:

$$u_{\beta }·n=0,K_{\beta }n-(K_{\beta }n·n)n=0.$$

(57)

At all other outer boundaries of the LWC and GFC and of the MEA, no slip conditions are applied. Initial conditions are set as $P_{\mathrm{liq}}=P_{\mathrm{gas}}^{\mathrm{B}}+P_{\mathrm{c}}^{\mathrm{B}}$, $P_{\mathrm{gas}}=P_{\mathrm{gas}}^{\mathrm{B}}$ and zero velocities everywhere.

2.7.2. Transport of Gas Species

At the anode and cathode GFC inlets, the boundary condition for water vapor is defined as

$${\omega }_{{\mathrm{H}}_2\mathrm{O}}^{\mathrm{AGFC}/\mathrm{CGFC},\mathrm{in}}=y_{{\mathrm{H}}_2\mathrm{O}}^{\mathrm{AGFC}/\mathrm{CGFC},\mathrm{in}}\frac{M_{{\mathrm{H}}_2\mathrm{O}}}{M_{\mathrm{n},\mathrm{AGFC}/\mathrm{CGFC},\mathrm{in}}},y_{{\mathrm{H}}_2\mathrm{O}}^{\mathrm{AGFC}/\mathrm{CGFC},\mathrm{in}}={\mathrm{RH}}^{\mathrm{AGFC}/\mathrm{CGFC},\mathrm{in}}\frac{{P_{\mathrm{sat}}}_{|T^{\mathrm{B}}}}{P_{\mathrm{gas}}^{\mathrm{B}}}.$$

(58)

At the respective GFC channel outlets, an outflow boundary condition is applied, which requires no diffusive flux across the boundaries. At all other outer boundaries, a no flux condition is applied. Initial water vapor mass fraction are set to ${10}^{-4}$ considering the completely dry conditions at the GFC inlets.

2.7.3. Transport of Heat

The temperature is set to $T^{\mathrm{B}}$ at the anode and cathode FCP ends (see Figure 1b). At all other outer boundaries, a no flux condition is applied. The initial temperature for the first iteration of the nonlinear solver is set to $T^{\mathrm{B}}$ in the entire model domain as the temperature drop.

2.7.4. Transport of Dissolved Water

As an initial condition for the nonlinear solver, we set $\lambda =10$. A no flux condition is applied at all outer boundaries of the ACL, PEM and CCL domains.

3. Results and Discussion

First, the model results are presented for the base case (Section 3.1). Then, we investigate the role of the evaporation mass transfer coefficient (Section 3.2). In Section 3.3, model results are compared to the analytical model that is based on the experimental results of Cochet et al. [10]. Finally, model limitations are discussed in Section 3.4.

3.1. Base Case

3.1.1. Key Output Parameters

Figure 3 provides an overview of the base case simulation results in terms of key output parameters. The integrated evaporation rate adds up to 82 mg h${}^{-1}$ (Figure 3a). From this amount of generated water vapor, a majority of 75 mg h${}^{-1}$ leaves the cell via the AGFC and 7 mg h${}^{-1}$ are transported through the membrane and leave the cell via the CGFCs (Figure 3a). The average relative humidity is higher in the anode domains than in the respective cathode domains (Figure 3b). The gas leaves the cell via the AGFC at an average relative humidity of 0.2. The relative humidity increases only very slightly along the CGFC because the high gas channel speed keeps the CGFC dry. A through-plane gradient in dissolved water content develops, with a higher average value in the ACL than in the CCL (Figure 3c). The simulated heat flux is reported as an average over the anode and cathode FCP ends (Figure 3d), similarly to the measured values in the experiments of Cochet et al. [10]. It is around 0.9 W cm${}^{-2}$ over the AFCP and 0.3 W cm${}^{-2}$ over the CFCP. This sums up to around 1.2 W cm${}^{-2}$, which is approximately equal to the heat consumed by evaporation (Figure 3d). The heat generated/consumed by latent heat of ab-/desorption is essentially zero as the heat generated by absorption in the ACL cancels the heat consumed by desorption in the CCL. Mesh convergence of these results is demonstrated in Appendix C.

3.1.2. Water Transport

The water is transported in liquid, vapor and dissolved phase, and is subject to evaporation and ab-/desorption phase changes. Liquid water from the LWC enters the MEA via the HPIL (Figure 4d). This is due to the higher relative permeability of the HPIL in comparison to the hydrophobic AGDL as a consequence of the distribution of liquid water saturation. The reduced liquid water saturation is around 0.74 in the HPIL and sharply transitions to zero in the hydrophobic AGDL and the other domains of the MEA (Figure 4d). This follows from the distribution of the capillary pressure (Figure 4c), the experimentally obtained water retention curves of the hydrophilic and hydrophobic parts and the cutoff for saturations smaller than the immobile liquid water saturation (Figure A1). Within the HPIL, liquid water is transported towards the volume below the AGFC along the gradient in liquid water pressure (Figure 4a) induced by evaporation.

The evaporation rate is the highest at the AGFC/HPIL interface (Figure 5b,c). Furthermore, evaporation takes place along the HPIL/AGDL interfaces with highest rates close to the AGFC. Generally, the evaporation rate increases towards the AGFC inlet (Figure 5b,h). This is because the relative humidity increases in flow direction as the AGFC accommodates increasingly more water vapor (Figure 5a,e). Close to the interfaces of the HPIL with the AGFC and AGDL, an evaporation film forms with a thickness of 0.7 $\mathsf{\mu }$m (see zoom in Figure 5). Over this thickness, the evaporation rate decreases from its maximum at the respective interfaces to zero inside the HPIL (Figure 5g,h), while relative humidity drops from being fully saturated inside the HPIL to unsaturated at the interfaces (Figure 5d,e).

The transport paths of the generated water vapor are visualized with streamlines in Figure 6a and the spatial components of the total flux magnitude vector are shown in Figure 6b–m. Water vapor is transported not only in through-plane direction away from the AGFC/HPIL interface (Figure 6b,e) but also in y-direction (i.e., parallel to channel flow) away from the HPIL/AGDL interfaces (Figure 6c,i) and in z-direction (i.e., perpendicular to channel flow) towards the AGFC (Figure 6d,m). Correspondingly, the relative humidity decreases from the HPIL, where the gas is essentially fully saturated, towards the AGFC (Figure 5a). From the 75 mg h${}^{-1}$ leaving via the AGFC, 49 mg h${}^{-1}$ enter the AGFC via the interface with the HPIL (Figure 3a, orange streamlines in Figure 6a). The other 26 mg h${}^{-1}$ enter the AGFC via the two interfaces with the hydrophobic AGDL and originate from the HPIL side walls (cyan streamlines in Figure 6a) and bottom wall (red streamlines in Figure 6a). The majority of the water vapor flux crossing the two AGFC/AGDL interfaces is close to the HPIL and channel walls (Figure 6b,h). Inside the AGFC, water vapor is advected with the gas stream towards the AGFC outlet (Figure 6c,f). Diffusion leads to a contribution of the flux towards the inlet (Figure 6f). This contribution is dominant only very close to the walls of the AGFC (Figure 6c) where the gas speed is the lowest.

The part of the water vapor generated at the bottom HPIL wall in larger distance to the AGFC is transported to the ACL where it is absorbed into the ionomer. In the absence of liquid water, the absorption and the resulting dissolved water content is controlled by the vapor-equilibrated hydration number, which is a non-linear function of relative humidity (Equation (51)). Therefore, the distribution of dissolved water content follows mainly the distribution of relative humidity (Figure 5a,c). The highest values are located at the projected intersection of the HPIL with the ribs and LWC. The lowest values can be found in the region of the projected AGFC/AGDL interface. The dissolved water content is significantly smaller than the vapor-equilibrated hydration number (Figure 3c). This is due to the mass transport resistance introduced by the sorption source term (Equation (48)) and the low absorption mass transfer coefficient (Equation (49)). The diffusive flux of dissolved water in the PEM is the strongest in through-plane direction with the in-plane components being one order of magnitude smaller. In the CCL, the dissolved water content follows a similar distribution as in the ACL at a lower level (Figure 5c). Water vapor is desorbed from the ionomer in the CCL and transported towards the CGFCs where it is advected with the gas stream to the outlets (Figure 6a–d, bottom). The absolute difference between average values of the dissolved water content and vapor-equilibrated hydration number is smaller in the CCL than in the ACL (Figure 3c) due to higher overall mass transfer coefficients of desorption than absorption (Equation (49)). Therefore, the mass transport resistance due to absorption in the ACL is higher than due to desorption in the CCL. In Section 3.3 we discuss the role of the mass transfer coefficient for ab-/desorption for the mass transport resistance across the membrane. On the cathode side of the MEA, the relative humidity is the highest below the ribs due to the dry and fast gas flow in the CGFCs while the influence of the HPIL is not clearly visible.

The transport mechanisms are influenced by the generation of a large amount of water vapor, which has a significantly higher molar mass than the carrier gas hydrogen. Corresponding to the distribution of relative humidity (Figure 5a), the gas density and pressure are higher inside the MEA than in the AGFC (Figure 4b). Consequently, a convective flow develops from the HPIL towards the AGFC. The profiles of the water flux vector components in Figure 6e–m show that the convective contribution to the total water vapor flux is mostly smaller but of the same order of magnitude than the diffusive contribution. Due to the absorption of water vapor into the ionomer of the ACL, the gas pressure is locally decreased such that an additional convective flow is directed from the HPIL to the ACL. The desorption of water vapor leads to an increase in gas pressure in the CCL (Figure 4b). The variations in gas pressure, however, are significantly less than on the anode side due to the smaller amount of water vapor and smaller difference in molar mass between water vapor and nitrogen. Therefore, diffusion is the dominant transport mechanism for water vapor on the cathode side of the MEA towards the CGFCs.

3.1.3. Heat Transport and Feedbacks

The temperature decreases due to the consumption of latent heat of evaporation. It reaches a minimum temperature of around 62 °C, which is located at the center of the AGFC/HPIL interface (Figure 7). This temperature drop from the operating temperature of 80 °C leads to a reduction in the strongly temperature-dependent saturation pressure (Equation (30)) by a factor of 2.2. This leads in turn to a reduction in the local evaporation rate. Close to the channel side walls there is a sharp temperature drop due to the contact of the channel with the highly thermally conductive FCP (Figure 7d), resulting in a high local heat flux and maximum evaporation rates close to the channel walls (Figure 5i). Within the evaporation film, the temperature changes only slightly and monotonically (Figure 7b,c). A comparison to an isothermal simulation shows that the overall evaporation rate is reduced from 156 mg h${}^{-1}$ to 82 mg h${}^{-1}$. This demonstrates the importance of accounting for the transport of heat in an evaporative cooling cell.

The heat flux across the FCP ends is the highest in the regions within the ribs (Figure 7a) because the thermal conductivity in the FCP is much higher than thermal conductivity of the gas species in the GFCs. This channel-rib effect and local strong temperature gradients close to the AGFC/HPIL interface could be important for evaluating thermal stresses and related degradation mechanism.

3.2. Role of the Evaporation Mass Transfer Coefficient

The macroscopic formulation of the source term for evaporation (Equation (27)) relies on the evaporation mass transport coefficient ${\alpha }_{\mathrm{ec}}$, which we varied by means of the factor $\chi$ (Equation (29)). It has been shown in previous studies that at high enough mass transport coefficients, total evaporation rates saturate [16]. In our study, a plateau at high values of ${\alpha }_{\mathrm{ec}}$ is nearly reached for the water vapor flux across the AGFC/AGDL interfaces and across the membrane via the CGFC, the average relative humidity in the porous domains and the average dissolved water content (Figure 3a–c). At the bottom wall of the HPIL, relative humidity reaches 100% saturation at $\chi =10$ (Figure 8a). The through-plane profile of relative humidity underneath the HPIL does not change when increasing $\chi$ to 100. Altogether, this indicates that a mass transport limitation is reached at $\chi =10$ for water vapor in the MEA outside the HPIL.

In contrast, there is still a considerable increase with $\chi$ in the evaporation rate, water vapor flux across the AGFC/HPIL interface, the net water vapor flux via the AGFC, the average relative humidity in the AGFC and average heat flux across the FCP ends for $\chi$ larger than 10 (Figure 3a,b,d). The local evaporation rate at the AGFC/HPIL interface increases strongly with $\chi$ (Figure 8b). Furthermore, an increase in $\chi$ still leads to an increase in relative humidity in the AGFC above the HPIL (Figure 8a) and to an increase in the through-plane water vapor flux across the AGFC/HPIL interface (Figure 8c). Thus, higher values of $\chi$ than 100 would be required to reach a transport limitation for the water vapor generated at the AGFC/HPIL interface and, hence, a plateau in total evaporation rates and heat fluxes.

3.3. Comparison to the Analytical Model Based on Experiments

We have focused in this study on a small part of the experimental setup of Cochet et al. [10] to identify the main water and heat transport mechanisms. This limits the direct comparison to these experiments where evaporation takes place along more AGFC/HPIL interfaces. As we apply zero relative humidity at the inlets of the AGFCs, our base case model could be considered to represent the behavior of the first AGFC/HPIL interface in channel flow direction and along the HPIL. It therefore simulates an extreme sub-case of the experiments of Cochet et al. [10] because simulated local evaporation rates and heat fluxes are higher than the average experimental ones as relative humidity would increase in channel flow direction. Furthermore, the simulated membrane humidification, relative humidity on the cathode side and water vapor flux via the CGFC are lower than the experimental values since in the experiments, more water vapor is produced in total and the AGFC would become saturated at some distance from the channel inlets leading to a larger water vapor flux towards the CGFC.

The analytical model based on the experimental results of Cochet et al. [10] and its assumptions are described in Appendix A. To allow for an indirect comparison between numerical and experimental results, we use the same parameterization and geometry for the analytical model as for the base case numerical model. The results from this analytical model is in overall agreement with the numerical results (Figure 3). The mass transfer coefficient computed in the analytical model (Equation (A5)) has the same order of magnitude as in the base case model. The amount of water vapor transported to the cathode side is significantly higher than in the base case numerical model. To better understand the controls on the mass transport resistance of the ionomer, we investigated the role of the mass transfer coefficients for ab-/desorption. When increasing the both mass transfer coefficients by a factor of 11, water vapor mass flow rates leaving the CGFCs increase by 5.6 mg h${}^{-1}$ to 12.6 mg h${}^{-1}$ (Figure 9a) thereby reducing the difference between the numerical model and the analytical model. At the same time, total evaporation rates increase by 4.7 mg h${}^{-1}$ such that the amount of water vapor leaving the AGFC is only slightly reduced by less than 1 mg h${}^{-1}$. An increase in the mass transfer coefficient for ab-/desorption leads to higher average relative humidity in the porous domains on the cathode side (Figure 9b). Furthermore, the average dissolved water content increases while the absolute difference between the dissolved water content and its equilibrium state decreases, showing that the mass transfer resistance for water vapor across the membrane is reduced.

Our modeling results furthermore show that the analytical model could be improved considering that

• the temperature drop due to evaporation is likely not negligible also in the experiments and, hence, should be included in the calculation of the saturation pressure and the evaporation rate;
• evaporation at the side and bottom interfaces of the HPIL with the hydrophobic AGDL contributes to the total evaporation rate and to the water vapor flux towards the GFC outlets;
• the mass transport resistance due to ab-/desorption into/from the ionomer adds to the total mass transport resistance of water vapor towards the CGFC outlets.

3.4. Model Limitations

The applicability of the presented macro-homogeneous modeling approach to simulate the two-phase flow in the GDL has been debated in the past due to the lack of a significant length scale separation between the thickness and pore sizes of a GDL. In particular, it is an open questions whether it is valid to locally apply capillary pressure-saturation relationships [12], considering that these relationships are obtained for the entire hydrophobic/hydrophilic GDL at a given capillary pressure. A rather unrealistic result of using the capillary pressure-saturation relationship from [8] is the simulated maximum liquid water saturation of around 74% in the HPIL. This might be an artifact of the experimental determination in the hydrophilic materials [41] as 100% saturation is expected in the entire HPIL.

The presented results demonstrate the limitations of using the up-scaling of the Hertz-Knudsen-Schrage equation from the interfacial area to the continuum source term. This up-scaling has been typically applied for hydrophobic GDLs. The situation in an evaporative cooling cell is more complicated considering the transition between hydrophilic and hydrophobic portions of the modified GDL and a direct contact of the HPIL with the AGFC. In future, pore-level simulations such as those of Safi et al. [16] could help to improve this up-scaling in this specific situation.

A manually built complex mesh with very fine resolution close to the HPIL walls is necessary to resolve sharp gradients within the evaporation film whose thickness decreases less than 1 $\mathsf{\mu }$m at the largest investigated values of $\chi$ (Figure 8b,d). This high local resolution could be avoided by developing a thin-film approximation of the evaporation-induced mass and heat transport with the help of new ex-situ experiments. This thin-film approximation would allow to simulate larger domains with more HPILs and gas channels and to compare directly to the experiments of Cochet et al. [10].

Furthermore, the presented model does not account for the strong anisotropy of the gas diffusivity, permeability and thermal conductivity in GDLs. Future models of an evaporation cooling cell should include these anisotropic GDL properties to better represent the in-plane mass and heat transport.

4. Conclusions

We have developed a 3-D macro-homogeneous model of a non-operating evaporatively cooled PEMFC to better understand the mass and heat transport processes inside such a cell. This model includes all layers of the experimental thermal test cell of Cochet et al. [10] but focuses on one AGFC/HPIL interface. Model results are compared to the analytical model of Cochet et al. [10].

Our base case simulation results show that most of the water evaporates within a thin film adjacent to the AGFC/HPIL interface. The side and bottom walls of the HPIL adjacent to the hydrophobic AGDL contribute to the total evaporation rate and are important for the water vapor flux towards the AGFC and via the membrane towards the cathode side. This could be taken into account for an improvement of the analytical model of Cochet et al. [10], which considers only evaporation along the AGFC/HPIL interfaces.

The model sheds light on the water transport paths and mechanism as well as on the humidity distribution of a non-operating evaporative cooling cell. Most of the generated water vapor leaves the cell via the outlet of the AGFC while only a small amount is transported through the membrane via the outlets of the CGFCs. We demonstrate that this is partly due to the mass transfer coefficient for ab-/desorption inducing a mass transport resistance across the membrane. This could be included in the analytical model of Cochet et al. [10] to improve the calculation of the mass transport coefficient for the transport of water vapor towards the CGFC outlets. The cell is significantly more humid on the anode side than on the cathode side in terms of relative humidity in the porous domains and dissolved water content in the ionomer. Furthermore, the humidification on the anode side varies considerably between highest values below the HPIL and the rib and lowest values below the gas channel. Evaporation generates a large amount of water vapor, which has a larger molar mass than the carrier gas hydrogen. This causes significant variations in gas pressure and gas density and in turn induces non-negligible convective currents towards the AGFC. Therefore, the assumption of constant gas pressure and diffusion as a dominant transport mechanism in the GDL typically made in macroscopic PEMFC models is not valid. More work is needed to investigate the water transport paths and humidity distribution in a complete non-operating cell with several HPILs and GFCs.

The model demonstrates the need to account for the non-isothermal effects due to a large temperature drop induced by evaporation, including a strong feedback on the local evaporation rate. The total heat flux into the cell is 1.2 W cm${}^{-2}$. This is comparable to the heat produced by electrochemical reactions in an operating fuel cell. An extension of this model to include electrochemical reactions is desired to better understand the mass and heat transport in an evaporative cooling cell under operation. This is especially important at high current densities when high amount of water and heat is produced in the CCLs.

We varied the mass transport coefficient for evaporation. At the highest values investigated, a mass transport limitation is reached inside the MEA. In contrast, the transport of water vapor through the AGFC is still limited by the evaporation kinetics. Higher mass transport coefficients lead to thinner evaporation films, which requires finer meshes to resolve the stronger local gradients. This currently limits a further increase in mass transport coefficients due to computational limitations. The development of a thin-film approximation with the help of pore-level simulations and ex-situ experiments could improve the general treatment of evaporation and to allow for a direct comparison with experiments of evaporation cooling cells such as those of Cochet et al. [10]. Further model improvements could include anisotropic GDL properties.