Modeling of the Transport Phenomena in Passive Direct Methanol Fuel Cells Using a Two-Phase Anisotropic Model

Abstract

The transport phenomena in a passive direct methanol fuel cell (DMFC) were numerically simulated by the proposed two-dimensional two-phase nonisothermal mass transport model. The anisotropic transport characteristic and deformation of the gas diffusion layer (GDL) were considered in this model. The natural convection boundary conditions were adopted for the transport of methanol, oxygen, and heat at the GDL outer surface. The effect of methanol concentration in the reservoir on cell performance was examined. The distribution of multiphysical fields in the membrane electrode assembly (MEA), especially in the catalyst layers (CLs), was obtained and analyzed. The results indicated that transport resistance for the methanol mainly existed in the MEA while that for oxygen and heat was primarily due to natural convection at the GDL outer surface. Because of the relatively high methanol concentration, the local reaction rate in CLs was mainly determined by the overpotential. Methanol concentration between 3 M and 4 M was recommended for passive liquid feed DMFC in order to achieve a balance between the cell performance and the methanol crossover.

1. Introduction

Direct methanol fuel cells (DMFCs) are electrochemical devices which can directly covert chemical energy into electricity. The passive DMFC is regarded as one of the most promising candidates to power portable devices due to its system simplicity as it has no parasitic power losses caused by the liquid pump or oxidant blower [1–3]. In the operation of a passive DMFC, methanol solution is stored in a small tank, and the cathode side of the cell is directly exposed to the atmosphere. Thus, the supply of methanol solution and oxygen to the reaction sites in catalyst layers is primarily based on the natural convection. In the past two decades, numerous scientists devoted their efforts to the commercialization of the passive DMFC. However, the wide application of this technique is still hindered by many practical problems [4, 5], such as low electroactivity and high cost of the catalyst, methanol crossover, and water management. More efforts should be made to investigate the fundamentals of the passive DMFC system. And computational modeling is a powerful and economical tool to study the intrinsically coupled transport processes in the DMFC which are much more difficult to be studied by the experimental testing.

Many numerical models have been reported for DMFCs. And a few of them focused on the passive DMFCs [6–13]. The continued improvements of these models shed light on the mechanism of inherently coupled transport processes by removing assumptions ever used. Earlier models are mainly single-phase models [9, 14–16]. In recent years, many two-phase models have been developed for DMFCs [6–8, 17–22] which can more realistically reflect the two-phase mass transport processes in the MEA. Ye et al. [7] developed a 2D two-phase model to study the methanol and water crossover in an air-breathing DMFC. Z. H. Wang and C. Y. Wang's two-phase model [18] based on the mixture multiphase flow approach was regarded as a classical numerical model for the active DMFC and was followed by many researchers. In their model, the CLs were treated as interfaces for simplicity. Yang and Zhao [21] presented a two-phase model for liquid feed DMFCs with a microagglomerate model accounting for the transfer process of oxygen from the gas pore to the active sites of the catalyst. This model had the clear physical sense and was easy to be applied in numerical codes. As performance of the DMFC was sensitive to the cell temperature, heat transport was modeled in some papers [9, 12, 23]. In previous studies, the liquid saturation was simply treated as continuously distributed through different layer of the cell. However, when researchers studied the water management issues in DMFC, they revealed a sudden change of liquid saturation at the interface between two different layers [11, 24–26].

GDLs contribute most volume of the MEA and have great effects on cell performance. GDLs are typically made of highly porous carbon cloth or carbon paper. The random overlapping structure of these materials makes the GDL exhibit strong anisotropic transport properties. For example, electrons transfer much faster along the in-plane direction than the through-plane direction. Additionally, inhomogeneous deformation of GDLs under the compression of flow plates will also change species transport properties in it. Up to now, many numerical models for PEMFC [27–33] and active DMFCs [34] considered anisotropic effects of GDLs. Pasaogullari et al. [27] presented a nonisothermal two-phase model for the cathode side of a PEM fuel cell. Their research focuses on the coupled heat and water transport in the GDL. Their results showed that the inherent anisotropy has a great effect on the temperature distribution and water transport. Himanen et al. [28–31] carried out a series of work to study the effect of inhomogeneous compression of GDL. Deformation affected physical properties are tested from experimental rig. Then, coefficients obtained from experimental data were adopted in their numerical models. Calculation of fluid flow their models are limited in single-phase model.

In our former works, we also developed a 2D two-phase transport model for active DMFCs [22]. This model was built on the classical liquid and vapor two-phase flow model in porous media. All layers in the MEA, including the GDLs, CLs, proton exchange membrane (PEM), and the microporous layers (MPLs) were considered as porous media with finite width. We improved this model later by taking account of the anisotropic transfer properties of the GDLs and the detailed multistep reaction mechanism of the methanol oxidation reaction (MOR) in anode CL [34]. Then, this modified model was used to investigate the mass, heat, and charges transport phenomena in the active DMFC [23, 26].

In the present work, we want to extend our numerical model for active DMFC to the passive DMFCs. Several boundary conditions are modified. And the anisotropic transport properties of the GDL, including the inherent anisotropy and deformation of the GDL, are taken into account in the model. A modified agglomerate model is used to consider the transport process of oxygen in cathode CL. The physical fields related to key transport processes, such as distribution of overpotentials, species concentration, local reaction rate, and temperature, are numerically studied.

2. Mathematical Model

A 2D physical domain is sketched in Figure 1(a). It represents the typical seven-layer structure of a MEA in passive DMFC which consists of a pair of GDLs, MPLs, and CLs on both the anode and cathode and a PEM. The MEA is sandwiched between two current collecting plates. Since both channel width and rib width in the current collecting plate are symmetrical with respect to their middle points, only a half-rib width and a half-channel width were considered. Figure 1(b) shows the inhomogeneous deformation of the GDL under the clamping force. The mathematical expressions [28–31, 35, 36] of the deformed surface and the structure-dependent parameters, such as the porosity, permeability, electrical conductivity, and interfacial electrical contact resistance can be found in Table 1. More details on the model development can be found elsewhere [26, 34]. Figure 2 gives computational subdomains of each solving variable. Detailed information of governing equations for solving variables related to the subdomains is given as follows.

Table 1:

Expressions of source terms and coefficients of the governing equations.

Parameters	Expressions
Mole generation rate of species	${\dot{R}}_{\mathrm{M}}=\{\begin{array}{cc}-i_{\mathrm{a}}/(\mathrm{6}F)\hfill & \text{ACL}\hfill \\ -i_{\mathrm{p}}/(\mathrm{6}F)\hfill & \text{CCL}\hfill \end{array}$ , ${\dot{R}}_{{\text{O}}_{\text{2}}}=\begin{array}{cc}-i_{\mathrm{c}}/(\mathrm{4}F)& \text{CCL}\end{array}$ , ${\dot{R}}_{\text{MV}}=\begin{array}{cc}{\tilde{R}}_{\mathrm{M}}& \text{ACL},\text{MPL},\text{ADL}\end{array}$ , ${\dot{R}}_{\text{WV}}=\begin{array}{cc}{\tilde{R}}_{\mathrm{W}}& \text{GDLs},\text{MPLs},\text{CLs}\end{array}$ , ${\dot{R}}_{\text{C}{\text{O}}_{\text{2}}}=\{\begin{array}{cc}{i}_{\mathrm{a}}/(\mathrm{6}F)\hfill & \text{ACL}\hfill \\ i_{\mathrm{p}}/(\mathrm{6}F)\hfill & \text{CCL}\hfill \end{array}$
Generation rate of mass in gas and liquid phases	${\dot{m}}_{\mathrm{g}}=\{\begin{array}{cc}{M}_{\mathrm{W}}{\tilde{R}}_{\mathrm{W}}+M_{\mathrm{M}}{\tilde{R}}_{\text{MV}}\hfill & \text{ADL},\text{AMPL}\hfill \\ M_{\mathrm{W}}{\tilde{R}}_{\mathrm{W}}+M_{\mathrm{M}}{\tilde{R}}_{\text{MV}}+M_{\text{C}{\text{O}}_{\text{2}}}{\tilde{R}}_{\text{C}{\text{O}}_{\text{2}}}\hfill & \text{ACL}\hfill \\ M_{\mathrm{W}}{\tilde{R}}_{\mathrm{W}}-M_{{\text{O}}_{\text{2}}}{\tilde{R}}_{{\text{O}}_{\text{2}}}+M_{\text{C}{\text{O}}_{\text{2}}}{\tilde{R}}_{\text{C}{\text{O}}_{\text{2}}}\hfill & \text{CCL}\hfill \\ M_{\mathrm{W}}{\tilde{R}}_{\mathrm{W}}\hfill & \text{CDL},\text{CMPL}\hfill \end{array}$ , ${\dot{m}}_{\mathrm{l}}=\{\begin{array}{cc}-{\dot{m}}_{\mathrm{g}}\hfill & \text{ADL},\text{AMPL}\hfill \\ -{\dot{m}}_{\mathrm{g}}-i/F\hfill & \text{ACL}\hfill \\ -{\dot{m}}_{\mathrm{g}}+i/F\hfill & \text{CCL}\hfill \\ -{\dot{m}}_{\mathrm{g}}\hfill & \text{CDL},\text{CMPL}\hfill \end{array}$
Generation rate of heat	${\dot{R}}_T=\{\begin{array}{cc}{\left|\mathbf{I}\right|}^{\mathrm{2}}/{\sigma }_{\mathrm{s}}+h_{\mathrm{v},\mathrm{W}}{\dot{R}}_{\mathrm{W}}+h_{\mathrm{v},\mathrm{M}}{\dot{R}}_{\mathrm{M}}\hfill & \text{ADL},\text{AMPL}\hfill \\ i_{\mathrm{a}}({\eta }_{\mathrm{a}}-(\mathrm{\Delta }{H}_{\mathrm{a}}-\mathrm{\Delta }{G}_{\mathrm{a}})/\mathrm{6}F)+{\left|\mathbf{I}\right|}^{\mathrm{2}}/{\sigma }_{\mathrm{s}}+h_{\mathrm{v},\mathrm{W}}{\dot{R}}_{\mathrm{W}}+h_{\mathrm{v},\mathrm{M}}{\dot{R}}_{\mathrm{M}}\hfill & \text{ACL}\hfill \\ {\left|\mathbf{I}\right|}^{\mathrm{2}}/{\sigma }_{\mathrm{m}}\hfill & \text{PEM}\hfill \\ (i_{\mathrm{c}}-i_{\mathrm{p}})({\eta }_{\mathrm{c}}-(\mathrm{\Delta }{H}_{\mathrm{c}}-\mathrm{\Delta }{G}_{\mathrm{c}})/\mathrm{4}F)-i_{\mathrm{p}}\mathrm{\Delta }H/\mathrm{6}F+{\left|\mathbf{I}\right|}^{\mathrm{2}}/{\sigma }_{\mathrm{s}}+h_{\mathrm{v},\mathrm{W}}{\dot{R}}_{\mathrm{W}}+h_{\mathrm{v},\mathrm{M}}{\dot{R}}_{\mathrm{M}}\hfill & \text{CCL}\hfill \\ {\left|\mathbf{I}\right|}^{\mathrm{2}}/{\sigma }_{\mathrm{s}}+h_{\mathrm{v},\mathrm{W}}{\dot{R}}_{\mathrm{W}}+h_{\mathrm{v},\mathrm{M}}{\dot{R}}_{\mathrm{M}}\hfill & \text{CDL},\text{CMPL}\hfill \end{array}$
Generation rate of charges	$i_{\mathrm{m}}=\{\begin{array}{cc}{i}_{\mathrm{a}}\hfill & \text{ACL}\hfill \\ i_{\mathrm{c}}-i_{\mathrm{p}}\hfill & \text{CCL}\hfill \end{array}$ , $i_{\mathrm{s}}=\{\begin{array}{cc}-i_{\mathrm{a}}\hfill & \text{ACL}\hfill \\ i_{\mathrm{p}}-i_{\mathrm{c}}\hfill & \text{CCL}\hfill \end{array}$
Relative permeabilities	krl = s3, krg = (1 − s)3
Thickness of the GDL under inhomogeneous compression	$L\left(y\right)=\{\begin{array}{cc}{\mathrm{L}}_{\text{comp\hspace{0.17em}}}=\mathrm{2.5}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{4}}\text{\hspace{0.17em}m}\hfill & \text{under-rib}\hfill \\ \mathrm{19.30314}\log \left(\right(y-\mathrm{0.0005})\times \mathrm{1}{\mathrm{0}}^{\mathrm{6}}+\mathrm{1})\times \mathrm{1}{\mathrm{0}}^{-\mathrm{6}}+{\mathrm{L}}_{\text{comp\hspace{0.17em}}}\hfill & \text{under-channel}\hfill \end{array}$
Porosity of the GDL under inhomogeneous compression	ε(y) = 1 − (1 − ε0)(L0/L(y))
In-plane permeability	$K_{\text{in}-p}=5.\times {10}^{-13}\left(L{\left(y\right)}^3/{\left[1-L\left(y\right)\right]}^2\right)$
Electrical conductivity of the GDL	σin-p = − 1.159 × 107L(y) + 6.896 × 103, σth-p = − 8.385 × 106L(y) + 3.285 × 103
Interfacial electrical contact resistance	re,MPL/CL = 5.83 × 10−10exp⁡(2.06 × 104L(y))

Figure 1:

Schematic of (a) the computational domain and (b) deformation of the GDL under a clamping force.

Figure 2:

Specified computational subdomains for the solving variables.

2.1. Pressures of Liquid and Gas Phases

In the liquid and gas two-phase flow in the MEA of passive DMFCs, pressures are continuous in both in-plane and through-plane directions. The pressure governing equations of liquid and gas phases can be given as

$\frac{\partial }{\partial x}\left(-\frac{K_x{\rho }_g{k}_{\text{rg}}}{{\mu }_g}\frac{\partial p_g}{\partial x}\right)+\frac{\partial }{\partial y}\left(-\frac{K_y{\rho }_g{k}_{\text{rg}}}{{\mu }_g}\frac{\partial p_g}{\partial y}\right)={\dot{m}}_g,$

(1)

$\frac{\partial }{\partial x}\left(-\frac{K_x{\rho }_l{k}_{\text{rl}}}{{\mu }_l}\frac{\partial p_l}{\partial x}\right)+\frac{\partial }{\partial y}\left(-\frac{K_y{\rho }_l{k}_{\text{rl}}}{{\mu }_l}\frac{\partial p_l}{\partial y}\right)={\dot{m}}_l-\frac{n_d{M}_W}{F}{i}_m,$

(2)

where K _x and K _y denote absolute permeability of the GDL on through-plane and in-plane direction, respectively. And i _m in (2) is the local generation rate of proton in the anode and cathode CLs. In the classical two-phase flow theory in porous media, the liquid-phase and gas-phase pressures are related to the capillary pressure:

$p_c=p_g-p_l=\sigma \cos {\theta }_c{\left(\frac{\epsilon }{K}\right)}^{\mathrm{0.5}}J\left(s\right),$

(3)

where σ and θ _c denote the surface tension between liquid and gas phases and the contact angle of the porous media, respectively. The capillary pressure is very sensitive to liquid saturation. The Leverett function, J(s), is commonly employed to account for the effect of liquid saturation on the capillary pressure:

$J\left(s\right)=\{\begin{array}{cc}\mathrm{1.417}\left(\mathrm{1}-s\right)-\mathrm{2.12}{\left(\mathrm{1}-s\right)}^{\mathrm{2}}+\mathrm{1.263}{\left(\mathrm{1}-s\right)}^{\mathrm{3}}\hfill & \hfill \\ \hspace{1em}{\mathrm{0}}^{°}\le {\theta }_c<\mathrm{9}{\mathrm{0}}^{°}\hfill & \hfill \\ \mathrm{1.417}s-\mathrm{2.12}{s}^{\mathrm{2}}+\mathrm{1.263}{s}^{\mathrm{3}}\hfill & \hfill \\ \hspace{1em}\mathrm{9}{\mathrm{0}}^{°}\le {\theta }_c<\mathrm{18}{\mathrm{0}}^{°}.\hfill & \hfill \end{array}$

(4)

When solving the equations in numerical program, the liquid saturation can be calculated through (3)–(4) after pressures of liquid and gas phases have been obtained from (1)–(2).

2.2. Species Concentrations

The conservation of species transport in liquid and gas phases are, respectively, written as

$\begin{array}{l}\frac{\partial }{\partial x}\left(u_l{C}_M\right)+\frac{\partial }{\partial y}\left(v_l{C}_M\right)\\ =\frac{\partial }{\partial x}\left(D_{M,x}^{\text{eff}}\frac{\partial C_M}{\partial x}\right)+\frac{\partial }{\partial y}\left(D_{M,y}^{\text{eff}}\frac{\partial C_M}{\partial y}\right)+{\dot{R}}_M,\end{array}$

(5)

$\begin{array}{l}\frac{\partial }{\partial x}\left(u_g{C}_{i,g}\right)+\frac{\partial }{\partial y}\left(v_g{C}_{i,g}\right)\\ =\frac{\partial }{\partial x}\left(D_{i,x}^{\text{eff}}\frac{\partial C_{i,g}}{\partial x}\right)+\frac{\partial }{\partial y}\left(D_{i,y}^{\text{eff}}\frac{\partial C_{i,g}}{\partial y}\right)+{\dot{R}}_{i,g},\end{array}$

(6)

where i in (6) represents different gas species. Note that the velocities of liquid and gas phases can be calculated from the gradient of phase pressures. The diffusion coefficients are expressed as

$D_j^{\text{eff}}=D_j^{\mathrm{0}}{\epsilon }^{\mathrm{1.5}}{\left(\mathrm{1}.-s\right)}^{\mathrm{1.5}}\hspace{1em}\left(\text{CL},\text{MPL}\right),$

(7)

$D_j^{\text{eff}}=D_j^{\mathrm{0}}\epsilon {\left(\frac{\epsilon -{\epsilon }_{\mathrm{0}}}{\mathrm{1}.-{\epsilon }_{\mathrm{0}}}\right)}^{\alpha }{\left(\mathrm{1}.-s\right)}^{\mathrm{1.5}}\hspace{1em}\left(\text{GDL}\right),$

(8)

where the constant α is 0.521 and 0.785 for in-plane and through-plane directions, respectively. Equation (7) is limited in CLs and MPLs. For the calculation of diffusion coefficients in GDLs, the anisotropy of GDLs is considered, and the expression is modified as (8) [37–39].

The nonequilibrium interfacial mass transfer of water and methanol between the liquid and gas phases is also taken into account in the present model. Simply, the interfacial species transfer rates are expressed as follows [20, 40]:

$\begin{array}{c}{\tilde{R}}_W=\{\begin{array}{cc}{k}_e\frac{\epsilon s{\rho }_l}{M_W}(p_W^{\mathrm{sat}}-y_{\text{WV}}{p}_g),\hfill & p_W^{\mathrm{sat}}\ge y_{\text{WV}}{p}_g,\hfill \\ k_c\frac{\epsilon \left(\mathrm{1}-s\right)y_{\text{WV}}}{RT}(p_W^{\mathrm{sat}}-y_{\text{WV}}{p}_g),\hfill & p_W^{\mathrm{sat}}<y_{\text{WV}}{p}_g,\hfill \end{array}\\ {\tilde{R}}_M=A_{\text{lg}}{h}_{\text{lg}}s\left(\mathrm{1}.-s\right)\frac{p_{\text{MV}}^{\mathrm{sat}}-p_{\text{MV}}}{RT},\end{array}$

(9)

where p _W ^sat⁡ and p_MV^sat⁡, respectively, denote the saturation pressure of water vapor and methanol vapor.

2.3. Heat Transfer

For the heat transfer process through the porous region, the thermal equilibrium between the liquid and gas phases is assumed, the general conservation of energy in the porous regions can be written as

$\begin{array}{l}\frac{\partial }{\partial x}({\rho }_g{C}_{p,g}{u}_{g}T+{\rho }_l{C}_{p,l}{u}_{l}T)+\frac{\partial }{\partial x}\left({\rho }_g{C}_{p,g}{v}_{g}T+{\rho }_l{C}_{p,l}{v}_{l}T\right)\\ =\frac{\partial }{\partial x}\left({\kappa }_x^{\text{eff}}\frac{\partial T}{\partial x}\right)+\frac{\partial }{\partial y}\left({\kappa }_y^{\text{eff}}\frac{\partial T}{\partial y}\right)+{\dot{R}}_T,\end{array}$

(10)

where C _p , ${\kappa }_{}^{\text{eff}}$ , and ${\dot{R}}_T$ represent the heat capacity, thermal conductivity, and local heat generation rate in the DMFC. Note that (10) is also applied to model the heat transfer in the membrane region.

2.4. Charge Transport

Charges in a passive DMFC include protons and electrons. Governing equations for their transport are given as

$\begin{array}{c}\frac{\partial }{\partial x}\left({\sigma }_{m,x}\frac{\partial {\phi }_m}{\partial x}\right)+\frac{\partial }{\partial y}\left({\sigma }_{m,y}\frac{\partial {\phi }_m}{\partial y}\right)=i_m,\\ \frac{\partial }{\partial x}\left({\sigma }_{s,x}\frac{\partial {\phi }_s}{\partial x}\right)+\frac{\partial }{\partial y}\left({\sigma }_{s,y}\frac{\partial {\phi }_s}{\partial y}\right)=i_s.\end{array}$

(11)

2.5. Electrochemical Kinetics

In CLs of the passive DMFC, the Tafel-like expressions are used to describe the kinetics of the methanol oxidation reaction (MOR) and oxygen reduction reaction (ORR):

$i_a=s{A}_a{i}_M^{\text{ref}}{\left(\frac{C_M}{C_M^{\text{ref}}}\right)}^{\gamma }\exp \left(\frac{{\alpha }_{a}F}{RT}{\eta }_a\right),$

(12)

$\gamma =\{\begin{array}{cc}\mathrm{0},\hfill & C_M{\xi }_M>C_M^{\text{ref}},\hfill \\ \mathrm{1},\hfill & C_M{\xi }_M\le C_M^{\text{ref}},\hfill \end{array}$

(13)

$i_c=\left(\mathrm{1}.-s\right)A_c{i}_{{\text{O}}_{\mathrm{2}}}^{\text{ref}}\left(\frac{C_{{\text{O}}_{\mathrm{2}}}}{C_{{\text{O}}_{\mathrm{2}}}^{\text{ref}}}\right)\exp \left(\frac{{\alpha }_{c}F}{RT}{\eta }_c\right){\xi }_{{\text{O}}_{\mathrm{2}}}.$

(14)

The modification factor ξ_O2 in (14) is obtained from the agglomerate model [22]. Overpotentials in (12) and (14) are determined by profiles of solid phase potential and electrolyte phase potential in CLs, which are calculated from (11).

$\begin{array}{c}{\eta }_a={\phi }_{s,a}-{\phi }_{m,a},\\ {\eta }_c=V_{\mathrm{0}}+{\phi }_{m,c}-{\phi }_{s,c}.\end{array}$

(15)

2.6. Current Balance and Cell Voltage

In an operating fuel cell, current balance should be satisfied. The mean current densities at the anode and cathode electrodes are

$I_a=\frac{{\iint }_{\text{ACL}}^{}{i}_{a}dx\hspace{0.17em}dy}{h_c+h_r},$

(16)

$I_c=\frac{{\iint }_{\text{CCL}}^{}{i}_{c}dx\hspace{0.17em}dy}{h_c+h_r}.$

(17)

Note that the “parasitic” current density due to methanol crossover, given as I _p = 6FN _M , is involved in the integration expression in (17). The flux of methanol crossover is calculated from

$N_M=-D_{M,N}\nabla C_{M,N}-\left(\frac{K_{\text{mem}}}{{\mu }_l}\nabla p_l\right)C_{M,N}+n_{d,M}\frac{I}{F}.$

(18)

Hence, the cell output current density and cell voltage can be determined by

$\begin{array}{c}{I}_{\text{Cell}}=I_a=I_c-I_p,\\ V_{\text{Cell}}={\phi }_{s,c}^{\mathrm{0}}.\end{array}$

(19)

Up to this point, all the governing equations have been presented. Detailed expressions of the source terms and coefficients used in these equations are listed in Tables 1 and 2.

Table 2:

Physicochemical properties and parameters used in simulation.

Parameters	Symbols	Value	Unit
Thickness of
GDL	l GDL	3.8 × 10−4	m
MPL	l GDL	0.3 × 10−4	m
CL	l cl	0.2 × 10−4	m
PEM	l pem	1.3 × 10−4	m
Height of a half-rib	h r	0.5 × 10−3	m
Height of a half-channel	h c	0.5 × 10−3	m
Environmental temperature	T	298.15	K
Environmental pressure	p	101325	pa
Nafion volume fraction in ACL and CCL	εN,acl, εN,ccl	0.4, 0.4
Conductivity in membrane phase	σm	$\left(\mathrm{0.5139}\lambda -\mathrm{0.326}\right)e^{\mathrm{1268}\left(\left(\mathrm{1}/\mathrm{303}\right)-\left(\mathrm{1}/T\right)\right)}$	Ω−1·m−1
Conductivity in CLs	σ cl	300	Ω−1·m−1
Permeability of
GDL	K th-p	3. × 10−12	m2
MPL	K MPL	7. × 10−13
CLs	K cl	3 × 10−14	m2
PEM	K pem	7 × 10−18	m2
Contact angle of
GDLs	θc,GDL	120°
MPLs	θc,MPL	135°
CLs	θc,CL	95°
Viscosity of
gas phase	μg	2.03 × 10−5	kg m−1 s−1
liquid phase	μl	4.06 × 10−4	kg m−1 s−1
Electroosmotic coefficients of water	n d	$\mathrm{2.5}\times \frac{\lambda }{\mathrm{22}}$
and methanol	n d.M	n d x M
Diffusivity of
Methanol in liquid water	D M,l	10−5.4163 −999.778/T	m2 s−1
Methanol in Nafion	D M,N	$\mathrm{4.9}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{10}}{e}^{\left[\mathrm{2436}(\mathrm{1}/\mathrm{333}-\mathrm{1}/T)\right]}$	m2 s−1
Oxygen in gas phase	$D_{{\text{O}}_{\text{2}},\mathrm{g}}$	1.775 × 10−5(T/273.15)1.823	m2 s−1
Oxygen in Nafion	$D_{{\text{O}}_{\text{2}},\mathrm{N}}$	1.844 × 10−10	m2 s−1
Water vapor in gas	D WV,g	2.56 × 10−5(T/307.15)2.334	m2 s−1
Methanol vapor in gas	D MV,g	$\begin{array}{c}-\mathrm{6.954}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{6}}+\mathrm{4.5986}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{8}}T\\ +\hspace{0.17em}\mathrm{9.4979}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{11}}{T}^{\mathrm{2}}\end{array}$	m2 s−1
Liquid water in the membrane	D W,N	4.17 × 10−8(1 + 161e−λ)e−2436/T	m2 s−1
Interfacial transfer rate constant for methanol	h lg	0.001	m2 s−1
Specific interfacial area between liquid and gas	A lg	1 × 105	m−1
Henry law constant for methanol	k H,M	0.096e0.04511(T − 273)	atm
Henry law constant for oxygen	k H	0.3125 × 101325/(RT)
Saturation pressure of vapor	log10⁡pWsat⁡	$\begin{array}{c}-2.1794+0.02953\left(T-273\right)\\ -9.1837\times {10}^{-5}{\left(T-273\right)}^2\\ +1.4454\times {10}^{-7}{\left(T-273\right)}^3\end{array}$	atm
Evaporation rate constant of water	k e	5 × 10−3	atm−1 s−1
Condensation rate constant of water	k c	50	s−1
Mole enthalpy change of the anode semireaction	ΔHa	126.69	kJ mol−1
Mole enthalpy change of the cathode semireaction	ΔHc	−567.41	kJ mol−1
Mole enthalpy change of the overall reaction	ΔH	−724.43	kJ mol−1
Mole Gibbs free energy change of the anode semireaction	ΔGa	9.35	kJ mol−1
Mole Gibbs free energy change of the cathode semireaction	ΔGc	−474.16	kJ mol−1
Specific heat capacity of
liquid methanol	C p,M	89.33	J mol−1 K−1
liquid water	C p,W	75.48	J mol−1 K−1
Oxygen	C p,O2	29.43	J mol−1 K−1
Carbon dioxide	C p,CO2	37.21	J mol−1 K−1
Water vapor	C p,WV	33.64	J mol−1 K−1
Electrochemical kinetics parameters
Exchange current density of ORR	$i_{{\text{O}}_{\text{2}}}^{\text{ref}}$	0.04222e(73200 × (1/353 – 1/T)/R)	A m−3
Reference concentration of oxygen	$C_{{\text{O}}_{\text{2}}}^{\text{ref}}$	pO2/RT	mol m−3
Transfer coefficient of cathode	αc	1.0
Thermodynamic voltage	V 0	1.21	V

2.7. Boundary Conditions

The boundary conditions at twelve interfaces, shown in Figure 1, should be specified. The conditions at the internal interfaces between adjacent components of the MEA are specified based on the continuity principle and flux balance of mass, charges, and heat. Detailed information about these boundary conditions can be found in our former work [34].

Boundary 1 is the interface between the methanol reservoir and the anode GDL. It is assumed that natural convection controls the methanol transport through this interface, so the boundary condition can be given as

$-D_M^{\text{eff}}\frac{d{C}_M}{dx}=h_M\left(C_M^{\infty }-C_M\right),$

(20)

where C _M ^∞ is the methanol concentration in the reservoir and C _M is the methanol concentration at the GDL surface. The natural convection coefficient h _M can be obtained from the natural convection correlation on a vertical surface [41]:

$\text{Sh}=\frac{h_{M}L}{D_M}=\mathrm{0.68}+\frac{\mathrm{0.67}\text{R}{\text{a}}_M^{\mathrm{0.25}}}{{\left[\mathrm{1}+{\left(\mathrm{0.492}/\text{S}{\text{c}}_M\right)}^{\mathrm{9}/\mathrm{16}}\right]}^{\mathrm{4}/\mathrm{9}}},$

(21)

where Sh is the Sherwood number and Ra and Sc are the Rayleigh number and Schmidt number, respectively:

$\begin{array}{c}\text{Ra}=\frac{\gamma g\left(w_M^{\infty }-w_M\right)L^{\mathrm{3}}}{\upsilon D_M},\\ \text{Sc}=\frac{\upsilon }{D_M}.\end{array}$

(22)

Boundary 5 represents the interface between the environment and the cathode GDL. Oxygen in the air transfers through this interface to the GDL. At this boundary, the condition of natural convection mass transport is also adopted:

$D_{{\text{O}}_{\mathrm{2}}}^{\text{eff}}\frac{d{C}_{{\text{O}}_{\mathrm{2}}}}{dx}=h_{{\text{O}}_{\mathrm{2}}}\left(C_{{\text{O}}_{\mathrm{2}}}^{\infty }-C_{{\text{O}}_{\mathrm{2}}}\right),$

(23)

where the calculation of the natural convection coefficient h_O2 is similar to h _M but with the physical properties of oxygen.

With respect to the boundary conditions for heat transfer, it is assumed that all the heat releases from the cathode electrode to the environment because the thermal resistance at the anode side is much larger. So the boundaries 1 and 2 are specified as thermal isolating boundaries for heat transfer. And conditions at boundaries 4 and 5 are given based on the mechanism of natural convection heat transfer:

$-{\lambda }_{\text{GDL}}\frac{dT}{dx}=h_{}\left(T-T^{\infty }\right),$

(24)

where h is the heat transfer coefficient and can be determined from

$\begin{array}{c}\text{Nu}=\frac{hL}{\lambda }=\mathrm{0.68}+\frac{\mathrm{0.67}\text{R}{\text{a}}_{}^{\mathrm{0.25}}}{{\left[\mathrm{1}+{\left(\mathrm{0.492}/\text{Pr}\right)}^{\mathrm{9}/\mathrm{16}}\right]}^{\mathrm{4}/\mathrm{9}}}\\ \text{Ra}=\text{GrPr}\\ \text{Gr}=\frac{\beta g\left(T-T^{\infty }\right)L^{\mathrm{3}}}{{\upsilon }^{\mathrm{2}}}.\end{array}$

(25)

3. Results and Discussion

Using a self-written computer code based on the finite-volume method, all the governing equations are numerically solved under the baseline operating conditions listed in Table 2. It is noticed that the present model is extended from the DMFC model presented in our previous work [22, 34], in which the details of the model validation against the experimental data can be found. In the following, the model is used to investigate the transport mechanisms in the MEA of a passive DMFC.

3.1. Effect of Methanol Concentration on the Cell Performance

Figure 3 gives the polarization curve of the passive DMFC feed with different methanol concentration. When methanol concentration in the reservoir is below 3 M, the cell performance improves significantly with the increase in methanol concentration. This is because the concentration polarization in the anode electrode caused by the methanol transport plays the dominant role to the cell performance at low methanol concentration. The concentration polarization regions appear apparently on the 1 M and 2 M polarization curves when the cell voltage is below 0.4 V and 0.3 V, respectively. The limiting current density increases from 436 A/m² to 1154 A/m² and 2305 A/m² as the methanol concentration increases from 1 M to 2 M and 3 M. And the highest cell power density upgrades rapidly from 142 W/m² to 270 W/m² and 330 W/m² with methanol concentration varying from 1 M to 2 M and 3 M, as shown in Figure 4. For the cases of moderate methanol concentration (3 M-4 M) and high methanol concentration (5 M), the methanol diffusing to the ACL is sufficient enough for the methanol oxidation reaction (MOR). Thus, no concentration polarization region can be observed on the polarization curves. As methanol concentration increases from 3 M to 4 M, the improvement of cell performance is very limited: increase in current density at 0.1 V from 2305 A/m² to 2602 A/m² and increase in the highest power density from 330 W/m² to only 331 W/m². When the methanol concentration is higher than 4 M, the cell performance will go down due to the severe methanol crossover from the anode side to the cathode catalyst layer.

Figure 3:

Polarization curve of the passive DMFC for different methanol feeding concentration.

Figure 4:

Output power density of the passive DMFC for different methanol feeding concentration.

The methanol permeating through the PEM will directly react with oxygen in the cathode CL, and the reaction rate is represented by the so-called “crossover current density.” It can be seen in Figure 5 that the crossover current density at open circuit condition is approximately proportional to methanol concentration. With the increase in cell current density, the crossover current density decreases sharply when the methanol concentration is below 3 M. At the limiting current density, it even drops to zero because of the mass transport limitation of methanol. When the methanol concentration is higher than 3 M, the decrease in crossover current density becomes slow. And the crossover current density even gradually increases with cell current density when the methanol concentration of 5 M is filled in the methanol reservoir.

Figure 5:

Methanol crossover current density of the passive DMFC for different methanol feeding concentration.

For the operation of a passive DMFC, the methanol crossover is a double-edged sword. On the one hand, methanol crossover can lead to a larger mixed potential at the cathode electrode, which results in the reduction of cell voltage and cell performance. It is also a waste of fuel. On the other hand, the directly oxidation of methanol at the cathode CL releases a lot of heat, which makes the cell run at the relatively high temperature level. And the activation of the catalyst upgrades rapidly, tending to improve the cell performance. When the methanol concentration is low (<3 M), the crossover current density is relatively small. The combined effects of mixed potential and temperature are exhibited as an improvement of the cell performance. However, when methanol concentration is as high as 5 M, the crossover current density is too large for the normal operation of the cell. The cell current density and the highest power density both become lower compared with those under low methanol concentration condition, which can be clearly seen in Figures 3 and 4.

In the present study, the temperature at the middle point of the PEM is defined as the operating temperature of the cell. Figure 6 gives the operating temperature of the passive DMFC with methanol concentration varying from 1 M to 5 M. Cell temperature at open circuit condition mainly depends on the methanol crossover rate. So it can be seen that a higher cell temperature at the open circuit condition corresponds to a higher methanol concentration. Generally speaking, the cell temperature increases with the increase in the cell current density. In addition to the heat generated from methanol crossover, the anode MOR and cathode ORR also release a lot of heat. As the cell current density becomes higher, the heat generation from anode and cathode reactions increases significantly. As a result, the cell temperature gradually goes up. For the 3 M case, with the cell current density increasing from 0 A/m² to 2305 A/m², the cell temperature increases from 301.8 K to 316.7 K, about 14.9 K increment. The increase of cell temperature for 1 M, 2 M, 4 M, and 5 M is 3.1 K, 7.9 K, 18.0 K, and 18.8 K, respectively.

Figure 6:

Operating temperature in the PEM of the passive DMFC for different methanol feeding concentration.

3.2. Distribution of Key Physical Fields in the MEA

The operation of a passive DMFC is sensitive to several key factors, such as the specie concentration, temperature, and liquid saturation. Figure 7 gives the distribution of methanol and oxygen concentrations in the MEA with 4 M methanol concentration in the reservoir and cell voltage of 0.1 V. As natural convection primarily controls the mass transfer process through the out surface of the GDL, the concentration difference can be seen at both the anode and cathode side. Methanol concentration at the out surface of the anode GDL is approximately 3.6 M, 10% lower than methanol concentration in the reservoir. The environment pressure is 101.32 kPa, and the corresponding oxygen concentration is about 8.58 mol/m³. The oxygen concentration at the cathode GDL surface is about 5.5 mol/m³, 35.9% lower than the oxygen concentration in the air. The methanol concentration drops significantly from the GDL to the PEM and from the under-channel region to the under-rib region. Comparing the concentration drop of methanol at the GDL surface with that in the MEA, we can find that the main mass transport resistance for methanol exists in the MEA. And the distribution of methanol concentration along the in-plane direction is extremely nonuniform. This will lead to the nonuniform distribution of local reaction rate of methanol in the anode CL, which is discussed in the next section. And it can also be deduced that the methanol crossover mainly occurs in the under-channel region as the methanol concentration in this region is much higher than that in the under-rib region. In contrast, the distribution of oxygen concentration in the cathode porous layers reflects that the concentration drop of oxygen at the cathode side is much smaller than that at the GDL surface. The largest concentration difference of oxygen at the cathode electrode is about 0.65 mol/m³, far below the 3.08 mol/m³ at the cathode GDL surface. This indicates that the nature convection primarily controls the mass transport of oxygen.

Figure 7:

Distribution of (a) methanol concentration at the anode electrode, PEM, and CCL, (b) oxygen concentration at the cathode electrode (cell voltage: 0.1 V).

The distribution of temperature throughout the MEA at 0.1 V and with 4 M methanol concentration is shown in Figure 8. The thermal contact resistance at the adjacent layers is neglected. It can be seen that the highest temperature appears in the under-channel region of the cathode CL. However, the temperature difference in the MEA is very small. The largest temperature difference is about 1 K. So it is reasonable that we use the temperature at the middle point of the PEM to represent the cell temperature. The environmental temperature is 298.15 K. Thus, the temperature difference at the surface of the MEA is relatively large, about 21.5 K.

Figure 8:

Distribution of temperature in MEA of the passive DMFC (cell voltage: 0.1 V).

Figure 9 gives the distribution of liquid saturation in the anode and cathode electrode. In order to guarantee the robustness of the numerical model, the liquid saturation at the GDL surfaces of anode and cathode electrode is set as 0.95 and 0.05, respectively. Apparently, in Figure 9, the liquid saturation at the anode side is much higher than that at the cathode side as the cathode GDL directly faces the air. In the anode GDL, the liquid saturation in the under-channel region is higher than that in the under-rib region. This means the transport of methanol and water mainly occurs in the under-channel region of the MEA. Due to the hydrophobic feature of the MPL, the liquid saturation in the MPLs is quite lower than those in other layers. At the cathode side, the volume of pores in cathode electrode is mainly occupied by the gas phase. Even in the cathode CL, the liquid saturation is below 0.1. The liquid saturation in the under-rib region is higher than that in the under-channel region because of the longer transport path.

Figure 9:

Distribution of liquid saturation in (a) anode electrode and (b) cathode electrode of the passive DMFC (cell voltage: 0.1 V).

3.3. Analysis of the Coupled Multiphysical Filed in the CLs

The CLs are the key components in a fuel cell. Electrochemical reactions and the transport of species, charges, and heat all occur in these two layers. And these complicated processes are coupled with each other. So the analysis of the coupled multiphysical processes in the CLs can help to understand the mechanism of the fuel cell. As the difference of temperature in the CLs is much small, the reaction rates of electrochemical reactions are mainly determined by three key factors: species concentration, overpotential, and the liquid saturation. The reaction rate in the cathode CL is related to (1. − s), as seen in (14). So, in this paper, we name (1. − s) as “gas saturation” for convenience. In this section, we examine the relations between the local reaction rate and the three key factors.

Figure 10 shows the distribution of local reaction rate of MOR, methanol concentration, overpotential, and liquid saturation in the anode CL. According to the fundamental of the MOR, the reference methanol concentration C _M ^ref in (12)–(13) is set as 100 mol/m³. When the methanol concentration in catalyst layer is higher than 100 mol/m³, the reaction order γ is zero. This means the methanol concentration is high enough to supply adequate reactant for the active sites of the catalyst particles. So the local reaction rate is determined by the other two factors: overpotential and liquid saturation. In Figure 10, it can be seen that the methanol concentration in most area of the anode CL is generally higher than 100 mol/m³ due to the high methanol concentration of 4 M in the reservoir. Only methanol concentration in a very small area (color blue) of the under-rib region is lower than 100 mol/m³. In our former work [22], we have revealed that the local reaction rate in the anode CL of an active DMFC feed with 1 M or 2 M methanol concentration is mainly determined by the local methanol concentration. In the small area in Figure 10 where the methanol concentration is below 100 mol/m³, the same trend can be seen. The liquid saturation in anode CL generally decreases along the x direction and increase along the y direction while the totally contrary trends can be seen in the distribution of overpotential. And the highest local reaction rate appears in an area (color red) where the overpotential and liquid saturation are both on a moderate level.

Figure 10:

Distribution of local reaction rate, methanol concentration, overpotential, and liquid saturation in ACL of the passive DMFC (cell voltage: 0.1 V).

The distribution of local reaction rate, oxygen concentration, overpotential, and gas saturation in the cathode CL is given in Figure 11. Compared with the physical fields in the anode CL, the distributions of the four physical fields in Figure 11 are relatively simple as they all exhibit the monotonous trends. The variations of gas saturation and the oxygen concentration in cathode CL are quite small and have very limited effect on the local reaction rate. It is clearly seen that the distribution of local reaction rate of ORR is basically determined by the distribution of the overpotential. The very similar distributions of local reaction rate and overpotential can be seen. This phenomenon also greens well with conclusion in our former work [22, 34].

Figure 11:

Distribution of local reaction rate, oxygen concentration, overpotential, and gas saturation in CCL of the passive DMFC (cell voltage: 0.1 V).

As a comparison, the distribution of physical fields in anode and cathode CLs at cell voltage of 0.2 V, at which the highest power density is achieved, are given in Figures 12 and 13. As the methanol concentration in all local areas of the anode CL at cell voltage of 0.2 V are higher than 100 mol/m³, it is seen that both the anode and cathode local reaction rate depend on the distribution of overpotential. And it is very interesting that the local reaction rate of MOR in the under-rib region of the anode CL is higher than that in the under-channel region. And this has negative effect on the operation of the passive DMFC as it will aggravate the methanol crossover in the under-channel region.

Figure 12:

Distribution of local reaction rate, methanol concentration, overpotential, and liquid saturation in ACL of the passive DMFC (cell voltage: 0.2 V).

Figure 13:

Distribution of local reaction rate, oxygen concentration, overpotential, and gas saturation in CCL of the passive DMFC (cell voltage: 0.2 V).

The distributions of local reaction rate in ACL for different methanol concentration at cell voltage of 0.2 V are also compared in Figure 14. It has mentioned in Section 3.1 that mass transport limitation for methanol exists for the operation of passive DMFC with 1 M and 2 M methanol solution. At this condition, methanol concentration in the anode CL dominates the distribution of the local reaction rate of MOR, shown in Figure 14. With the increase in methanol concentration, the effect of methanol concentration on the local reaction rate becomes weak. For the 3 M and 4 M condition, the local reaction rates are mainly determined by the overpotential which responds to the transport of protons in the PEM. Besides that, the local reaction rates in 3 M and 4 M cases distribute more uniform compared to 1 M and 2 M cases. The uniform local reaction rates are beneficial for the high utilization of the catalyst. As a comprehensive consideration of the overall cell performance and the distribution of multiphysical fields in CL, the methanol concentration in reservoir between 3 M and 4 M is recommended to achieve a balance between the cell performance and methanol crossover.

Figure 14:

Distribution of local reaction rate in ACL of the passive DMFC for different methanol feeding concentration (cell voltage: 0.2 V).

4. Conclusions

A 2D two-phase nonisothermal mass transport model for a passive DMFC is proposed in the present work with the consideration of the anisotropic transport characteristic and deformation of the GDL. This model is used to simulate the coupled transport phenomena in the MEA. The natural convection boundary conditions are adopted for the transport of methanol, oxygen, and heat at the GDL outer surface. Based on the numerical results, some conclusions are summarized as follows.

1. The sudden change of species concentration and temperature concentration exists at the outer surfaces of the GDL due to the natural convection mechanism. The transport resistance for the methanol mainly exists in the MEA while that for oxygen and heat is primarily at the GDL outer surface.

2. When the methanol concentration in the reservoir is relatively high (>3 M), the local reaction rate in CLs is mainly determined by the overpotential and distributes more uniform which can improve the utilization efficiency of the catalyst.

3. The methanol concentrations between 3 M and 4 M are recommended for the passive liquid feed DMFC in order to achieve a balance between the cell performance and the methanol crossover.