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}}_{M} = {\begin{cases} - i_{a} / (6 F) & ACL \\ - i_{p} / (6 F) & CCL \end{cases}$ , ${\dot{R}}_{O_{2}} = \begin{matrix} - i_{c} / (4 F) & CCL \end{matrix}$ , ${\dot{R}}_{MV} = \begin{matrix} {\tilde{R}}_{M} & ACL, MPL, ADL \end{matrix}$ , ${\dot{R}}_{WV} = \begin{matrix} {\tilde{R}}_{W} & GDLs, MPLs, CLs \end{matrix}$ , ${\dot{R}}_{C O_{2}} = {\begin{cases} i_{a} / (6 F) & ACL \\ i_{p} / (6 F) & CCL \end{cases}$

Generation rate of mass in gas and liquid phases	${\dot{m}}_{g} = {\begin{cases} M_{W} {\tilde{R}}_{W} + M_{M} {\tilde{R}}_{MV} & ADL, AMPL \\ M_{W} {\tilde{R}}_{W} + M_{M} {\tilde{R}}_{MV} + M_{C O_{2}} {\tilde{R}}_{C O_{2}} & ACL \\ M_{W} {\tilde{R}}_{W} - M_{O_{2}} {\tilde{R}}_{O_{2}} + M_{C O_{2}} {\tilde{R}}_{C O_{2}} & CCL \\ M_{W} {\tilde{R}}_{W} & CDL, CMPL \end{cases}$ , ${\dot{m}}_{l} = {\begin{cases} - {\dot{m}}_{g} & ADL, AMPL \\ - {\dot{m}}_{g} - i / F & ACL \\ - {\dot{m}}_{g} + i / F & CCL \\ - {\dot{m}}_{g} & CDL, CMPL \end{cases}$

Generation rate of heat	${\dot{R}}_{T} = {\begin{cases} {\| I \|}^{2} / σ_{s} + h_{v, W} {\dot{R}}_{W} + h_{v, M} {\dot{R}}_{M} & ADL, AMPL \\ i_{a} (η_{a} - (Δ H_{a} - Δ G_{a}) / 6 F) + {\| I \|}^{2} / σ_{s} + h_{v, W} {\dot{R}}_{W} + h_{v, M} {\dot{R}}_{M} & ACL \\ {\| I \|}^{2} / σ_{m} & PEM \\ (i_{c} - i_{p}) (η_{c} - (Δ H_{c} - Δ G_{c}) / 4 F) - i_{p} Δ H / 6 F + {\| I \|}^{2} / σ_{s} + h_{v, W} {\dot{R}}_{W} + h_{v, M} {\dot{R}}_{M} & CCL \\ {\| I \|}^{2} / σ_{s} + h_{v, W} {\dot{R}}_{W} + h_{v, M} {\dot{R}}_{M} & CDL, CMPL \end{cases}$

Generation rate of charges	$i_{m} = {\begin{cases} i_{a} & ACL \\ i_{c} - i_{p} & CCL \end{cases}$ , $i_{s} = {\begin{cases} - i_{a} & ACL \\ i_{p} - i_{c} & CCL \end{cases}$

Relative permeabilities	k_rl = s³, k_rg = (1 − s)³

Thickness of the GDL under inhomogeneous compression	$L (y) = {\begin{cases} L_{comp} = 2.5 \times 1 0^{- 4} m & under-rib \\ 19.30314 \log ((y - 0.0005) \times 1 0^{6} + 1) \times 1 0^{- 6} + L_{comp} & under-channel \end{cases}$

Porosity of the GDL under inhomogeneous compression	ε(y) = 1 − (1 − ε₀)(L₀/L(y))

In-plane permeability	$K_{in - p} = 5. \times 10^{- 13} (L {(y)}^{3} / {[1 - L (y)]}^{2})$

Electrical conductivity of the GDL	σ_in-p = − 1.159 × 10⁷L(y) + 6.896 × 10³, σ_th-p = − 8.385 × 10⁶L(y) + 3.285 × 10³

Interfacial electrical contact resistance	r_e,MPL/CL = 5.83 × 10⁻¹⁰exp⁡(2.06 × 10⁴L(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} (- \frac{K_{x} ρ_{g} k_{rg}}{μ_{g}} \frac{\partial p_{g}}{\partial x}) + \frac{\partial}{\partial y} (- \frac{K_{y} ρ_{g} k_{rg}}{μ_{g}} \frac{\partial p_{g}}{\partial y}) = {\dot{m}}_{g},

(1)

\frac{\partial}{\partial x} (- \frac{K_{x} ρ_{l} k_{rl}}{μ_{l}} \frac{\partial p_{l}}{\partial x}) + \frac{\partial}{\partial y} (- \frac{K_{y} ρ_{l} k_{rl}}{μ_{l}} \frac{\partial p_{l}}{\partial y}) = {\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} = σ \cos θ_{c} {(\frac{ε}{K})}^{0.5} J (s),

(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 (s) = {\begin{cases} 1.417 (1 - s) - 2.12 {(1 - s)}^{2} + 1.263 {(1 - s)}^{3} \\ 0^{°} \leq θ_{c} < 9 0^{°} \\ 1.417 s - 2.12 s^{2} + 1.263 s^{3} \\ 9 0^{°} \leq θ_{c} < 18 0^{°} . \end{cases}

(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} (u_{l} C_{M}) + \frac{\partial}{\partial y} (v_{l} C_{M}) \\ = \frac{\partial}{\partial x} (D_{M, x}^{eff} \frac{\partial C_{M}}{\partial x}) + \frac{\partial}{\partial y} (D_{M, y}^{eff} \frac{\partial C_{M}}{\partial y}) + {\dot{R}}_{M}, \end{array}

(5)

\begin{array}{l} \frac{\partial}{\partial x} (u_{g} C_{i, g}) + \frac{\partial}{\partial y} (v_{g} C_{i, g}) \\ = \frac{\partial}{\partial x} (D_{i, x}^{eff} \frac{\partial C_{i, g}}{\partial x}) + \frac{\partial}{\partial y} (D_{i, y}^{eff} \frac{\partial C_{i, g}}{\partial y}) + {\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}^{eff} = D_{j}^{0} ε^{1.5} {(1 . - s)}^{1.5} (CL, MPL),

(7)

D_{j}^{eff} = D_{j}^{0} ε {(\frac{ε - ε_{0}}{1 . - ε_{0}})}^{α} {(1 . - s)}^{1.5} (GDL),

(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{matrix} {\tilde{R}}_{W} = {\begin{cases} k_{e} \frac{ε s ρ_{l}}{M_{W}} (p_{W}^{sat} - y_{WV} p_{g}), & p_{W}^{sat} \geq y_{WV} p_{g}, \\ k_{c} \frac{ε (1 - s) y_{WV}}{R T} (p_{W}^{sat} - y_{WV} p_{g}), & p_{W}^{sat} < y_{WV} p_{g}, \end{cases} \\ {\tilde{R}}_{M} = A_{lg} h_{lg} s (1 . - s) \frac{p_{MV}^{sat} - p_{MV}}{R T}, \end{matrix}

(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} (ρ_{g} C_{p, g} u_{g} T + ρ_{l} C_{p, l} u_{l} T) + \frac{\partial}{\partial x} (ρ_{g} C_{p, g} v_{g} T + ρ_{l} C_{p, l} v_{l} T) \\ = \frac{\partial}{\partial x} (κ_{x}^{eff} \frac{\partial T}{\partial x}) + \frac{\partial}{\partial y} (κ_{y}^{eff} \frac{\partial T}{\partial y}) + {\dot{R}}_{T}, \end{array}

(10)

where C_p, $κ_{}^{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{matrix} \frac{\partial}{\partial x} (σ_{m, x} \frac{\partial φ_{m}}{\partial x}) + \frac{\partial}{\partial y} (σ_{m, y} \frac{\partial φ_{m}}{\partial y}) = i_{m}, \\ \frac{\partial}{\partial x} (σ_{s, x} \frac{\partial φ_{s}}{\partial x}) + \frac{\partial}{\partial y} (σ_{s, y} \frac{\partial φ_{s}}{\partial y}) = i_{s} . \end{matrix}

(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}^{ref} {(\frac{C_{M}}{C_{M}^{ref}})}^{γ} \exp (\frac{α_{a} F}{R T} η_{a}),

(12)

γ = {\begin{cases} 0, & C_{M} ξ_{M} > C_{M}^{ref}, \\ 1, & C_{M} ξ_{M} \leq C_{M}^{ref}, \end{cases}

(13)

i_{c} = (1 . - s) A_{c} i_{O_{2}}^{ref} (\frac{C_{O_{2}}}{C_{O_{2}}^{ref}}) \exp (\frac{α_{c} F}{R T} η_{c}) ξ_{O_{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{matrix} η_{a} = φ_{s, a} - φ_{m, a}, \\ η_{c} = V_{0} + φ_{m, c} - φ_{s, c} . \end{matrix}

(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_{ACL}^{} i_{a} d x d y}{h_{c} + h_{r}},

(16)

I_{c} = \frac{\iint_{CCL}^{} i_{c} d x d y}{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} - (\frac{K_{mem}}{μ_{l}} \nabla p_{l}) C_{M, N} + n_{d, M} \frac{I}{F} .

(18)

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

\begin{matrix} I_{Cell} = I_{a} = I_{c} - I_{p}, \\ V_{Cell} = φ_{s, c}^{0} . \end{matrix}

(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⁻⁴	m
MPL	l _GDL	0.3 × 10⁻⁴	m
CL	l _cl	0.2 × 10⁻⁴	m
PEM	l _pem	1.3 × 10⁻⁴	m
Height of a half-rib	h _r	0.5 × 10⁻³	m
Height of a half-channel	h _c	0.5 × 10⁻³	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	$(0.5139 λ - 0.326) e^{1268 ((1 / 303) - (1 / T))}$	Ω⁻¹·m⁻¹
Conductivity in CLs	σ_cl	300	Ω⁻¹·m⁻¹
Permeability of
GDL	K _th-p	3. × 10⁻¹²	m²
MPL	K _MPL	7. × 10⁻¹³
CLs	K _cl	3 × 10⁻¹⁴	m²
PEM	K _pem	7 × 10⁻¹⁸	m²
Contact angle of
GDLs	θ_c,GDL	120°
MPLs	θ_c,MPL	135°
CLs	θ_c,CL	95°
Viscosity of
gas phase	μ_g	2.03 × 10⁻⁵	kg m⁻¹ s⁻¹
liquid phase	μ_l	4.06 × 10⁻⁴	kg m⁻¹ s⁻¹
Electroosmotic coefficients of water	n _d	$2.5 \times \frac{λ}{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}	m² s⁻¹
Methanol in Nafion	D _M,N	$4.9 \times 1 0^{- 10} e^{[2436 (1 / 333 - 1 / T)]}$	m² s⁻¹
Oxygen in gas phase	$D_{O_{2}, g}$	1.775 × 10⁻⁵(T/273.15)^1.823	m² s⁻¹
Oxygen in Nafion	$D_{O_{2}, N}$	1.844 × 10⁻¹⁰	m² s⁻¹
Water vapor in gas	D _WV,g	2.56 × 10⁻⁵(T/307.15)^2.334	m² s⁻¹
Methanol vapor in gas	D _MV,g	$\begin{matrix} - 6.954 \times 1 0^{- 6} + 4.5986 \times 1 0^{- 8} T \\ + 9.4979 \times 1 0^{- 11} T^{2} \end{matrix}$	m² s⁻¹
Liquid water in the membrane	D _W,N	4.17 × 10⁻⁸(1 + 161e^−λ)e^−2436/T	m² s⁻¹
Interfacial transfer rate constant for methanol	h _lg	0.001	m² s⁻¹
Specific interfacial area between liquid and gas	A _lg	1 × 10⁵	m⁻¹
Henry law constant for methanol	k _H,M	0.096e^{0.04511(T − 273)}	atm
Henry law constant for oxygen	k _H	0.3125 × 101325/(RT)
Saturation pressure of vapor	log₁₀⁡p_W^sat⁡	$\begin{matrix} - 2.1794 + 0.02953 (T - 273) \\ - 9.1837 \times 10^{- 5} {(T - 273)}^{2} \\ + 1.4454 \times 10^{- 7} {(T - 273)}^{3} \end{matrix}$	atm
Evaporation rate constant of water	k _e	5 × 10⁻³	atm⁻¹ s⁻¹
Condensation rate constant of water	k _c	50	s⁻¹
Mole enthalpy change of the anode semireaction	ΔH_a	126.69	kJ mol⁻¹
Mole enthalpy change of the cathode semireaction	ΔH_c	−567.41	kJ mol⁻¹
Mole enthalpy change of the overall reaction	ΔH	−724.43	kJ mol⁻¹
Mole Gibbs free energy change of the anode semireaction	ΔG_a	9.35	kJ mol⁻¹
Mole Gibbs free energy change of the cathode semireaction	ΔG_c	−474.16	kJ mol⁻¹
Specific heat capacity of
liquid methanol	C _p,M	89.33	J mol⁻¹ K⁻¹
liquid water	C _p,W	75.48	J mol⁻¹ K⁻¹
Oxygen	C _p,O2	29.43	J mol⁻¹ K⁻¹
Carbon dioxide	C _p,CO2	37.21	J mol⁻¹ K⁻¹
Water vapor	C _p,WV	33.64	J mol⁻¹ K⁻¹
Electrochemical kinetics parameters
Exchange current density of ORR	$i_{O_{2}}^{ref}$	0.04222e_{(73200 × (1/353 – 1/T)/R)}	A m⁻³
Reference concentration of oxygen	$C_{O_{2}}^{ref}$	p_O2/RT	mol m⁻³
Transfer coefficient of cathode	α_c	1.0
Thermodynamic voltage	V ₀	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}^{eff} \frac{d C_{M}}{d x} = h_{M} (C_{M}^{\infty} - C_{M}),

(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]:

Sh = \frac{h_{M} L}{D_{M}} = 0.68 + \frac{0.67 R a_{M}^{0.25}}{{[1 + {(0.492 / S c_{M})}^{9 / 16}]}^{4 / 9}},

(21)

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

\begin{matrix} Ra = \frac{γ g (w_{M}^{\infty} - w_{M}) L^{3}}{υ D_{M}}, \\ Sc = \frac{υ}{D_{M}} . \end{matrix}

(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_{O_{2}}^{eff} \frac{d C_{O_{2}}}{d x} = h_{O_{2}} (C_{O_{2}}^{\infty} - C_{O_{2}}),

(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:

- λ_{GDL} \frac{d T}{d x} = h_{} (T - T^{\infty}),

(24)

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

\begin{matrix} Nu = \frac{h L}{λ} = 0.68 + \frac{0.67 R a_{}^{0.25}}{{[1 + {(0.492 / Pr)}^{9 / 16}]}^{4 / 9}} \\ Ra = GrPr \\ Gr = \frac{β g (T - T^{\infty}) L^{3}}{υ^{2}} . \end{matrix}

(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.

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.

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.

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.

Footnotes

Nomenclature

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (51306048), the National Natural Science Foundation of China of International Cooperation Project (51210011), and the Program for New Century Excellent Talents in University (NCET-13-0792).

References

Dyer

C. K.

, “Fuel cells for portable applications,” Journal of Power Sources, vol. 106, no. 1–2, pp. 31–34, 2002.

Zhao

T. S.

Chen

Yang

W. W.

, and Xu

, “Small direct methanol fuel cells with passive supply of reactants,” Journal of Power Sources, vol. 191, no. 2, pp. 185–202, 2009.

Bae

Kho

B. K.

Lim

T.-H.

I.-H.

Hong

S.-A.

, and Ha

H. Y.

, “Performance evaluation of passive DMFC single cells,” Journal of Power Sources, vol. 158, no. 2, pp. 1256–1261, 2006.

Zhao

T. S.

and Chen

, “Recent progress in understanding of coupled heat/mass transport and electrochemical reactions in fuel cells,” International Journal of Energy Research, vol. 35, no. 1, pp. 15–23, 2011.

Faghri

, and Bahrami

, “Recent advances in passive and semi-passive direct methanol fuel cells,” International Journal of Thermal Sciences, vol. 62, pp. 12–18, 2012.

Garvin

J. J.

and Meyers

J. P.

, “Modeling of coupled multiphase transport in direct methanol fuel cell diffusion layers,” Journal of the Electrochemical Society, vol. 158, no. 9, pp. B1119–B1127, 2011.

Zhu

Liao

, and Fu

, “Two-dimensional two-phase mass transport model for methanol and water crossover in air-breathing direct methanol fuel cells,” Journal of Power Sources, vol. 192, no. 2, pp. 502–514, 2009.

Guo

Chen

Y. P.

Xue

Y. Q.

, and Ma

C. F.

, “Three-dimensional transient modeling and analysis of two-phase mass transfer in air-breathing cathode of a fuel cell,” International Journal of Hydrogen Energy, vol. 38, pp. 11028–11037, 2013.

Chen

and Zhao

T. S.

, “Mathematical modeling of a passive-feed DMFC with heat transfer effect,” Journal of Power Sources, vol. 152, no. 1–2, pp. 122–130, 2005.

10.

Yeh

T.-K.

and Chen

C.-H.

, “Modeling and optimizing the performance of a passive direct methanol fuel cell,” Journal of Power Sources, vol. 175, no. 1, pp. 353–362, 2008.

11.

and Faghri

, “Water transport characteristics in a passive liquid-feed DMFC,” International Journal of Heat and Mass Transfer, vol. 53, no. 9–10, pp. 1951–1966, 2010.

12.

Oliveira

V. B.

Rangel

C. M.

, and Pinto

A. M. F. R.

, “One-dimensional and non-isothermal model for a passive DMFC,” Journal of Power Sources, vol. 196, no. 21, pp. 8973–8982, 2011.

13.

Yang

W. W.

Zhao

T. S.

, and Wu

Q. X.

, “Modeling of a passive DMFC operating with neat methanol,” International Journal of Hydrogen Energy, vol. 36, no. 11, pp. 6899–6913, 2011.

14.

and Liu

, “A three-dimensional mathematical model for liquid-fed direct methanol fuel cells,” Journal of Power Sources, vol. 160, no. 1, pp. 413–421, 2006.

15.

Kulikovsky

A. A.

, “Analytical model of the anode side of DMFC: the effect of non-Tafel kinetics on cell performance,” Electrochemistry Communications, vol. 5, no. 7, pp. 530–538, 2003.

16.

Zhao

T. S.

, and Xu

, “The role of under-rib convection in mass transport of methanol through the serpentine flow field and its neighboring porous layer in a DMFC,” Electrochimica Acta, vol. 51, no. 25, pp. 5420–5429, 2006.

17.

Yang

W. W.

and Zhao

T. S.

, “Two-phase, mass-transport model for direct methanol fuel cells with effect of non-equilibrium evaporation and condensation,” Journal of Power Sources, vol. 174, no. 1, pp. 136–147, 2007.

18.

Wang

Z. H.

and Wang

C. Y.

, “Mathematical modeling of liquid-feed direct methanol fuel cells,” Journal of the Electrochemical Society, vol. 150, no. 4, pp. A508–A519, 2003.

19.

and Liu

, “A three-dimensional two-phase flow model for a liquid-fed direct methanol fuel cell,” Journal of Power Sources, vol. 163, no. 2, pp. 907–915, 2007.

20.

Divisek

Fuhrmann

Gärtner

, and Jung

, “Performance modeling of a direct methanol fuel cell,” Journal of the Electrochemical Society, vol. 150, no. 6, pp. A811–A825, 2003.

21.

Yang

W. W.

and Zhao

T. S.

, “A two-dimensional, two-phase mass transport model for liquid-feed DMFCs,” Electrochimica Acta, vol. 52, no. 20, pp. 6125–6140, 2007.

22.

Miao

Y.-L.

X.-L.

, and Zou

J.-Q.

, “A two-dimensional two-phase mass transport model for direct methanol fuel cells adopting a modified agglomerate approach,” Journal of Power Sources, vol. 185, no. 2, pp. 1233–1246, 2008.

23.

Miao

Y.-L.

Zhao

T.-S.

, and Tao

W.-Q.

, “Numerical investigation of heat transport in a direct methanol fuel cell with anisotropic gas diffusion layers,” Frontiers in Heat and Mass Transfer, vol. 2, no. 1, Article ID 013001, 2011.

24.

Zhao

T. S.

, and Yang

W. W.

, “Modeling of water transport through the membrane electrode assembly for direct methanol fuel cells,” Journal of Power Sources, vol. 178, no. 1, pp. 291–308, 2008.

25.

Yang

W. W.

and Zhao

T. S.

, “Numerical investigations of effect of membrane electrode assembly structure on water crossover in a liquid-feed direct methanol fuel cell,” Journal of Power Sources, vol. 188, no. 2, pp. 433–446, 2009.

26.

Y.-L.

Miao

Zhao

T.-S.

, and Yang

W.-W.

, “Numerical study of the effect of the GDL structure on water crossover in a direct methanol fuel cell,” International Journal of Hydrogen Energy, vol. 37, no. 5, pp. 4422–4438, 2012.

27.

Pasaogullari

Mukherjee

P. P.

Wang

C.-Y.

, and Chen

K. S.

, “Anisotropic heat and water transport in a PEFC cathode gas diffusion layer,” Journal of the Electrochemical Society, vol. 154, no. 8, pp. B823–B834, 2007.

28.

Hottinen

and Himanen

, “PEMFC temperature distribution caused by inhomogeneous compression of GDL,” Electrochemistry Communications, vol. 9, no. 5, pp. 1047–1052, 2007.

29.

Hottinen

Himanen

Karvonen

, and Nitta

, “Inhomogeneous compression of PEMFC gas diffusion layer. Part II. Modeling the effect,” Journal of Power Sources, vol. 171, no. 1, pp. 113–121, 2007.

30.

Nitta

Himanen

, and Mikkola

, “Thermal conductivity and contact resistance of compressed gas diffusion layer of PEM fuel cell,” Fuel Cells, vol. 8, no. 2, pp. 111–119, 2008.

31.

Nitta

Karvonen

Himanen

, and Mikkola

, “Modelling the effect of inhomogeneous compression of GDL on local transport phenomena in a PEM fuel cell,” Fuel Cells, vol. 8, no. 6, pp. 410–421, 2008.

32.

, “Investigation of the effects of the anisotropy of gas-diffusion layers on heat and water transport in polymer electrolyte fuel cells,” Journal of Power Sources, vol. 191, no. 2, pp. 259–268, 2009.

33.

Yang

W. W.

Zhao

T. S.

, and He

Y. L.

, “Modelling of coupled electron and mass transport in anisotropic proton-exchange membrane fuel cell electrodes,” Journal of Power Sources, vol. 185, no. 2, pp. 765–775, 2008.

34.

Miao

Y.-L.

, and Zou

J.-Q.

, “Modeling the effect of anisotropy of gas diffusion layer on transport phenomena in a direct methanol fuel cell,” Journal of Power Sources, vol. 195, no. 11, pp. 3693–3708, 2010.

35.

Nitta

Himanen

, and Mikkola

, “Contact resistance between gas diffusion layer and catalyst layer of PEM fuel cell,” Electrochemistry Communications, vol. 10, no. 1, pp. 47–51, 2008.

36.

Nitta

Hottinen

Himanen

, and Mikkola

, “Inhomogeneous compression of PEMFC gas diffusion layer. Part I. Experimental,” Journal of Power Sources, vol. 171, no. 1, pp. 26–36, 2007.

37.

Nam

J. H.

and Kaviany

, “Effective diffusivity and water-saturation distribution in single- and two-layer PEMFC diffusion medium,” International Journal of Heat and Mass Transfer, vol. 46, no. 24, pp. 4595–4611, 2003.

38.

Tomadakis

M. M.

and Sotirchos

S. V.

, “Effective Knudsen diffusivities in structures of randomly overlapping fibers,” AIChE Journal, vol. 37, no. 1, pp. 74–86, 1991.

39.

Tomadakis

M. M.

and Sotirchos

S. V.

, “Ordinary and transition regime diffusion in random fiber structures,” AIChE Journal, vol. 39, no. 3, pp. 397–412, 1993.

40.

Murgia

Pisani

Shukla

A. K.

, and Scott

, “A numerical model of a liquid-feed solid polymer electrolyte DMFC and its experimental validation,” Journal of the Electrochemical Society, vol. 150, no. 9, pp. A1231–A1245, 2003.

41.

P.-W.

Zhang

Wang

Q.-M.

Schaefer

, and Chyu

M. K.

, “The performance of PEM fuel cells fed with oxygen through the free-convection mode,” Journal of Power Sources, vol. 114, no. 1, pp. 63–69, 2003.