Flow field optimization for proton exchange membrane fuel cells with varying channel heights and widths

Electrochimica Acta 54 (2009) 5522–5530

Contents lists available at ScienceDirect

Electrochimica Acta journal homepage: www.elsevier.com/locate/electacta

Flow field optimization for proton exchange membrane fuel cells with varying channel heights and widths Xiao-Dong Wang a , Yu-Xian Huang b , Chin-Hsiang Cheng c , Jiin-Yuh Jang b , Duu-Jong Lee d,∗ , Wei-Mon Yan e , Ay Su f a

Department of Thermal Engineering, School of Mechanical Engineering, University of Science and Technology Beijing, Beijing 100083, China Department of Mechanical Engineering, National Cheng Kung University, Tainan 70101, Taiwan Department of Aeronautics and Astronautics, National Cheng Kung University, Tainan 70101, Taiwan d Department of Chemical Engineering, National Taiwan University, Taipei 106, Taiwan e Department of Mechatronic Engineering, Huafan University, Taipei 22305, Taiwan f Department of Mechanical Engineering, Fuel Cell Center, Yuan Ze University, Taoyuan 300, Taiwan b c

a r t i c l e

i n f o

Article history: Received 1 March 2009 Received in revised form 14 April 2009 Accepted 18 April 2009 Available online 3 May 2009 Keywords: Optimization Flow field design Simplified conjugate-gradient method Serpentine flow field Sub-rib convection

a b s t r a c t The optimal cathode flow field design of a single serpentine proton exchange membrane fuel cell is obtained by adopting a combined optimization procedure including a simplified conjugate-gradient method (SCGM) and a completely three-dimensional, two-phase, non-isothermal fuel cell model. The cell output power density Pcell is the objective function to be maximized with channel heights, H1 –H5 , and channel widths, W2 –W5 as search variables. The optimal design has tapered channels 1, 3 and 4, and diverging channels 2 and 5, producing 22.51% increment compared with the basic design with all heights and widths setting as 1 mm. Reduced channel heights of channels 2–4 significantly enhance sub-rib convection to effectively transport oxygen to and liquid water out of diffusion layer. The final diverging channel prevents significant leakage of fuel to outlet via sub-rib convection from channel 4. Near-optimal design without huge loss in cell performance but is easily manufactured is discussed. © 2009 Elsevier Ltd. All rights reserved.

1. Introduction The flow field design in the bipolar plates affects the performance of the proton exchange membrane (PEM) fuel cell. Mathematical models and numerical simulations interpret the conjugated transport processes and electrochemical reactions occur in the PEM fuel cell [1–15]. Comprehensive reviews on the PEM fuel cell models are available [16–18]. Water flooding normally occurs at the cathode electrode of the PEM fuel cell. When a large amount of liquid water accumulates in the porous layer pores, the oxygen transport resistance increases and the oxygen mass flow rate decreases. The cathode flow field design is a key factor determining reactant and product transport rates and for removing liquid water from cell. Effect of geometric parameters for flow field design on cell performance was extensively studied [19–33]. Among the numerous flow field geometries, the serpentine flow field has long channels for oxygen (reactant) and liquid water (product) flows. Ideally the oxygen is fed at a flow rate that is just enough

∗ Corresponding author. Tel.: +886 2 23625632; fax: +886 2 23623040. E-mail address: [email protected] (D.-J. Lee). 0013-4686/$ – see front matter © 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.electacta.2009.04.051

to maintain sufficiently high chemical reaction rates in membrane without considerable leakage. A properly designed cathode flow field can effectively distribute fed oxygen to reaction sites and effective remove excess produced liquid water from the fuel cell, thereby yielding high power density at a given oxygen feed rate. Optimization methods are available [34–37]. The conjugategradient method, which stems from perturbation principles, transforms the optimal problem into three separate problems, namely, the direct problem, the sensitivity problem, and the adjoining problem, whose solutions the updated conjugate-gradient direction as well as the step size towards convergence [38,39]. Cheng and Chang [40] simplified the conjugated-gradient method to take constant step size in each iteration step and use a direct numerical differentiation method to determine the sensitivity and the search direction. This work proposed an optimization approach combining the simplified conjugate-gradient method refined from [40] and a three-dimensional, two-phase, non-isothermal fuel cell model refined form [41] to look for optimal flow field for single serpentine PEM fuel cell considering the minimum of reciprocal of cell power density as the objective function. Both channel heights and channel widths were the search variables.

X.-D. Wang et al. / Electrochimica Acta 54 (2009) 5522–5530

5523

Table 1 channel heights and widths for various search steps. Search step

H0 (mm)

H1 (mm)

H2 (mm)

H3 (mm)

H4 (mm)

H5 (mm)

W1 (mm)

W2 (mm)

W3 (mm)

W4 (mm)

W5 (mm)

1 3 6 10 22

1 1 1 1 1

1 0.84 0.54 0.28 0.25

1 0.82 0.52 0.46 0.35

1 0.82 0.52 0.26 0.10

1 0.82 0.52 0.12 0.10

1 0.95 0.79 0.80 0.69

1 1 1 1 1

1 0.92 0.68 0.87 1.19

1 0.92 0.81 1.21 1.52

1 0.90 0.60 0.63 0.85

1 1.26 1.91 1.28 0.44

2. Optimization method 2.1. The PEM fuel cell model The two-fluid method used in the present work was refined from that adopted in Wang et al. [41] to incorporate the heat effects using energy equations for entire cells. The model assumes that the system is steady; the inlet reactants are ideal gases; the flow is laminar; and the porous layers such as the diffusion layer, catalyst layer and PEM are isotropic. The model includes continuity, momentum and species equations for gaseous species, liquid water transport equations in the channels, gas diffusion layers, and catalyst layers, water transport equation in the membrane, electron and proton transport equations. The Bulter–Volmer equation was used to describe electrochemical reactions in the catalyst layers. The main governing equations are listed in Appendix B. The thermal conductivities of solid matrix in Eq. (B9) are 150 W m−1 K−1 for the gas diffusion layers, catalyst layers and membrane. The thermal conductivity of the gaseous mixture is calculated by mix kinetic theory [42] and the thermal conductivity of liquid water is assumed to be a function of temperature. The last three terms in Eq. (B8) represent electrical work, Joule heating, and latent heat due to phase change of water, respectively. The source terms (Si , Sj , and others) and other relevant physicochemical parameters in Eqs. (B1)–(B9) are listed in [41].

The magnitude of Ci can be different for each search variable in accordance with its sensitivity to the objective function. Note that the convergence speed of the iteration of the optimization process may be appreciably slowed down by a fixed step size; however, the need to determine the one-dimensional search of step size for each ˇi has been removed. Detailed information regarding the SCGM procedure is available in [40]. 2.3. Serpentine flow field with changing channel height and width A single-serpentine PEM fuel cell with an active area of 9 mm × 9 mm was selected as the testing example to display applicability of the present inverse problem method to optimization of multi-parameters of PEM fuel cell. The cell consists of cathode and anode flow channels, 0.4 mm cathode and anode gas diffusion layers, 0.005 mm cathode and anode catalyst layers, and a 0.035 mm proton exchange membrane. The cathode and anode have the same serpentine flow field with varying channel heights and widths subject to optimization (Fig. 1a). The origin and three axes were fixed at the left-bottom corner of the cell. Fig. 1b shows a schematic of the cathode flow field with serpentine flow field design. The ribs are numbered as ribs 1, 2, 3 and 4. The three x–z cross-sections, section A–C (y = 1.5, 4.5 and 7.5 mm) and four y–z cross-sections, section I–IV (at centers for ribs 1–4) are shown in Fig. 1b. The inlet height of the channel 1 is fixed (H0 = 1 mm) and the heights of all turns are maintained as those at the end of the channel before them

2.2. Simplified conjugate-gradient method The main equations of conventional conjugate-gradient method are listed in Appendix C. The objective function of a process to be optimized can be expressed as follows: J = F (X1 , X2 , · · ·XN )

(1)

where Xi (i = 1,2,...,N) denotes the parameters for optimization (referred to as search variable or design variable). In the present study, the objective function (J) is defined as the reciprocal of cell output power density, Pcell , the product of output voltage (Vcell ) and average current density (I), so an optimal set of geometric parameters (H1 – H5 and W2 –W5 described in Section 2.3) is searched for to reach the minimum of the objective function, or maximum of cell output power. (k) The value of ˇi can be obtained by implementing a one(k−1)

dimensional search with respect to ˇi along the negative gradient direction of the objective function under condition that (k−1) (k−1) (k) the other search variables are kept unchanged, or Xj + ˇj j

(j = 1, 2, · · ·N, j = / i). Therefore, when the objective function reaches a minimum in the process of the above one-dimensional search (k−1) procedures, the corresponding value of ˇi is the optimal search size of ˇi in the kth search step. Determination of optimal search step size with conventional conjugated-gradient procedures is a complicated process. In the simplified conjugate-gradient method [40], the values of the step size ˇi are claimed to fix at a constant value without loss of convergence. That is, ˇi = Ci

i = 1, 2, · · ·N

(2)

Fig. 1. The PEM fuel cell under investigation. (a) Schematic of the fuel cell with varying channel heights and channel widths. (b) Index for the studied fuel cell with ribs 1–4 and cross-sections for profiling the velocities distributions discussed in latter sections.

5524

X.-D. Wang et al. / Electrochimica Acta 54 (2009) 5522–5530

Fig. 2. Flow chart for the present optimization method.

(for instance, the height of turn between channel 1 and 2 is kept as that at the end of channel 1). The height of channel 1 varies linearly from H0 at y = 0 mm to H1 at y = 8 mm. And the height of the channel 2 varies linearly from H1 at y = 8 mm to H2 at y = 1 mm, and so on. The width of the channel 1 is fixed (W1 = 1 mm) and the widths of the channels 2–5, W2 –W5 , vary with W5 = 5 − W1 − W2 − W3 − W4 . Five heights, H1 –H5 , and three widths, W2 –W4 , are the objects to be optimized. The adopted numerical models are summarized in the following sections. 2.4. Optimization scheme The objective function for optimization including eight search variables, H1 –H5 and W2 –W4 , is defined as follows: J=

1 Pcell (H1 , H2 , H3 , H4 , H5 , W2 , W3 , W4 )

(3)

The goal is to look for the set of heights (H1 , H2 , H3 , H4 , H5 ) and widths (W2 , W3 , W4 ) to yield minimum of J. For the sake of convenience, we use X1 –X8 to denote H1 , H2 , H3 , H4 , H5 , W2 , W3 and W4 , respectively, in what follows.

Fig. 2 shows the flow charts of the adopted optimization scheme. The initial guess for each search variable is made first, and in the successive steps the conjugate-gradient coefficients and the search directions are evaluated to estimate the new search variables. This process repeated to reach the minimum of the objective function. Specifically, the procedures are listed as follows:

(1) Make the initial guess for the search variables X1 –X8 , assign the values to the step sizes ˇ1 –ˇ8 . (2) Create geometry and grids of PEM fuel cell model for the assigned X1 –X8 . Specify all boundary conditions, and then numerically solve Eqs. (B1)–(B9). (3) Calculate the objective function J(X1 ,X2 ,X3 ,X4 ,X5 ,X6 ,X7 ,X8 ) based on direct solver Eqs. (B1)–(B9) and Eq. (3). If the convergence criterion is satisfied, then terminate iteration; otherwise, proceed to step (4). (4) Calculate the sensitivity coefficients, ∂J/∂Xi , of the objective function for each search variable based on Eq. (C1). (k) (5) Calculate the conjugate-gradient coefficients i for the each search variable based on Eq. (C4). For the first step with k = 1, (1) i = 0 (i = 1, 2,. . ., 8).

X.-D. Wang et al. / Electrochimica Acta 54 (2009) 5522–5530

5525

Fig. 3. Power density (Pcell ) and search variables (H1 –H5 , W2 –W5 ) at various search steps. Fig. 5. Polarization curves and power density curves for basic and optimal designs.

were performed to ensure that the solution is independent of the adopted grid size. The tested operational conditions included: temperatures for inlet gas streams at anode and cathode channels were constant (323 K); the reactants on the anode side and on the cathode sides were 100% in humidity; the inlet flow rates of anode and cathode were 24 and 72 cm3 min−1 , respectively. An operation voltage Vcell = 0.4 V was selected in this optimization study. The surfaces of anode and cathode current collectors are subjected to convective condition with heat transfer coefficient h = 150 W m−2 K−1 . 3. Results 3.1. Geometrical optimization Fig. 4. Optimization process from different initial guess.

(k)

(6) Calculate the search directions, i , for the each search variable based on Eq. (C3). (k+1) (k) (k) = Hi − ˇi i and then (7) Update new search variables by Hi return to step (2). In steps (2) and (4), the coupled, non-linear model equations, Eqs. (B1)–(B9), were converted to the finite-difference form by the control volume method and were solved iteratively with an iteration criterion for convergence of 10−6 . Preliminary numerical tests

Fig. 3 shows the changes in cell output power density, Pcell = 1/J, and those in channel heights, H1 –H5 , and channel widths, W2 –W5 , at each search step. The initial guess assumed H0 = H1 = H2 = H3 = H4 = H5 = W1 = W2 = W3 = W4 = W5 = 1 mm, identical to the basic design with straight channels. The corresponding Pcell is 7260 W m−2 . The Pcell increases with the search steps. After 22 search steps the Pcell approaches a plateau of 8894 W m−2 , about 22.51% higher than that for basic design with straight channels. The same optimal design is achieved using another initial guess for channel heights (H1 = H2 = H3 = H4 = H5 = 1.5 mm, W1 = 1 mm, W2 = W3 = W4 = 0.5 mm, and W5 = 2.5 mm) (Fig. 4). The obtained

Fig. 6. Schematic of PEM fuel for various search steps. (a) Step 1 (basic design); (b) step 6; (c) step 10; and (d) step 22 (optimal design).

5526

X.-D. Wang et al. / Electrochimica Acta 54 (2009) 5522–5530

optimal design is thereby at least a local minimum solution for the objective function J. Fig. 5 shows the polarization curves and power density curves for the basic and optimal designs. For operating voltages greater than 0.7 V, both basic and optimal designs reveal similar cell performance. However, at operating voltages lower than 0.7 V, the optimal design performs much better than the basic design, with the difference increasing with decreasing operating voltage. Evolution of geometric shape of flow channel during optimization is demonstrated (Table 1 and Fig. 6). The optimal height set reads: H0 (1 mm)>H5 (0.69 mm)>H2 (0.35 mm)>H1 (0.25 mm)>H3 ∼ (0.10 mm), and W3 (1.52 mm)>W2 = H4 (1.19 mm)>W1 (1 mm)>W4 (0.85 mm)>W5 (0.44 mm). Restated, the

Fig. 7. Current density distributions along cell width direction in the membrane at various search steps. (a) Section A; (b) section B; and (c) section C.

channel displays a tapered characteristic for channels 1, 3 and 4, and a diverging characteristic in height for channels 2 and 5. 3.2. Power density On search step 6, H1 –H4 are all approximately 0.53 mm, while H5 is 0.79 mm. The widths of channels 2–4, W2 –W4 , decrease and that of channel 5 increases to 1.91 mm (W5 ) (Fig. 3). Restated, a wide and diverging final channel is displayed connecting to the outlet (Fig. 6b). The corresponding cell power density reached 8293 W m−2 , accounting for 14.23% increase based on basic design.

Fig. 8. Liquid water distributions along cell width direction on the interface between the cathode gas diffusion layer and catalyst layer. (a) Section A; (b) section B; and (c) section C.

X.-D. Wang et al. / Electrochimica Acta 54 (2009) 5522–5530

On the following steps (step 10 and step 22 as examples in Fig. 3), the channel 2–4 become wider, with the corresponding channel heights becoming lower, with a final, diverging channel with reduced width (Fig. 6c and d). The corresponding cell power densities on step 10 and step 22 are 8711 and 8894 W m−2 , respectively. The geometrical modification from that in Fig. 6b and c leads to ((8711–8293)/7260=) 5.76% further increase in power density based on basic design. Change of flow design from that in Fig. 6c and d only increase additional 2.52% power density based on basic design. 4. Discussion 4.1. Current and liquid water distributions Figs. 7a–c and 8a–c show respectively the distributions of current density in the membrane and liquid water at interface of diffusion layer and catalyst layer yielded by fuel cell with flow field design at various search steps. The current density for basic design (Fig. 6a) gradually decreases when moving downstream with local maximums occurring under the channels and local minimum under the ribs (Fig. 7). On sections A–C the liquid water concentrations increase along the cell width direction with maximums occurring under the ribs and minimums under the channels (Fig. 8). The local liquid content is reduced and the local current density is increased along with optimization (Figs. 7 and 8) with decreased magnitude of fluctuation. Restated, the flow field geometry is changed to lead to a uniform liquid water distribution at low content. Flow design optimization suggests to reduce the heights and width for channel 2–4 from basic design so to achieve low and uniform distributions of liquid water (Fig. 8). The reduction of chan-

Fig. 9. Pressure drops between the adjacent channels at various search steps. (a) Cross rib 1; and (b) cross rib 4.

5527

nel cross-sectional area for channels 2–4 and placing a big and diverging final channel produces electrical current distributions close to that for optimal design at x < 6 mm (step 6 in Fig. 7). This improvement leads to 14.23% increase of power density compared with the basic design. Further adjustment on geometry for channels 2–4 in step 10 and next only affects the local current distribution at x > 6 mm. This presents a regime of geometric parameters for channels 2–4 yielding similar local cell power density. This observation is beneficial to manufacturing aspect since a now too-difficult to manufacture flow channels but still with good enough cell power enhancement is desired. The optimal design suggests a diverging but narrow final channel (Fig. 6d). The width reduction for channel 5 significantly improves its liquid water removal capacity to outlet (Fig. 8), hence enhancing the local power density (Fig. 7). As Wang et al. [42] addressed, a diverging, final channel is preferred to minimize reactant leakage to outlet owing to the occurrence of strong sub-rib convection. 4.2. Sub-rib flow field For a PEM fuel cell, when there is pressure difference between two adjacent channels, the pressure difference will force the reactants in the channel to cross over the rib and reach to adjacent channel, which is referred to as the sub-rib convection. For the single serpentine flow field, the reactants enter in a single channel and after several turns flow out of the cell; thus, there is larger pressure difference between any two adjacent channels and the sub-rib convection occurs under every rib for the single serpentine flow field.

Fig. 10. Velocity of gaseous species on the y–z cross-section in the cathode diffusion layer. (a) Under rib 1; and (b) under rib 4.

5528

X.-D. Wang et al. / Electrochimica Acta 54 (2009) 5522–5530

Fig. 11. Schematic of PEM fuel cells with various flow fields. (a) Optimal design but setting H5 = 0.1 mm; and (b) optimal design but setting all Wi = 1 mm.

Fig. 9(a)–(b) shows the pressure drops across rub 1 and 4 at various search steps. Fig. 10(a)–(b) shows the corresponding distributions of sub-rib convection. For the basic design, the cross-rib pressure differences range 1–10 Pa, and the sub-rib convection velocities range 0.01–0.11 m s−1 . The portion of sub-rib diffusion layer close to channel turns has weak sub-rib convection, so the liquid water is easily accumulated to block the oxygen transport. Reactants are largely consumed and liquid water is accumulated along the flow path, the local current density is low in the second half-cell with the basic design. During optimization with reduced channel heights and widths, the pressure drops and convection cross the ribs are significantly increased. The velocities under ribs 1 and 4 are 0.14–0.54 m s−1 and 0.27–1.06 m s−1 respectively for the optimal design, accounting for 5–6 times increase by mean sub-rib convection velocities. The strong sub-rib convection flow effectively removes liquid water and promotes oxygen transport into the catalyst layer. The channel 5 linked to outlet should be designed to yield a not-too-strong sub-rib convection to preventing serious short-circuit flow to outlet.

Fig. 12. Current density distributions in the membrane for basic design (Fig. 5a), optimal design (Fig. 5d), design with H5 = 0.1 mm (Fig. 10a), and design with all Wi = 1 mm (Fig. 10b).

4.3. Near-optimal flow field design The optimal design suggested by the present combined method presents a complicated geometry with various channel height and channel width specification. The corresponding average electric current densities at 0.4 V of applied voltage for basic and optimal design are 18,150 and 22,236 A m−2 , respectively. The “near-optimal” design without big loss in cell performance but is easily manufactured is looked for. The flow field with optimal design but setting H5 = 1 mm (Fig. 11a), and that with all Wi = 1 mm were tested (Fig. 11b). The corresponding average electric current densities are 20,597 and 21,826 A m−2 , respectively. Restated, the use of a straight, final channel of 0.1 mm height has led to ((22,236–20,597)/22,236=) 7.37% power loss compared with that of an optimal design. Setting all channel widths to be 1 mm with channel heights of (H0 = 1 mm, H1 = 0.25 mm, H2 = 0.35 mm, H3 ∼ = H4 = 0.10 mm, H5 = 0.69 mm) yields only ((22,236–21,826)/22,236=) 1.68% loss of current density. Restated, a design with channel heights obtained by the combined optimization method but with constant channel widths (=1 mm) enhances power density by 20.25%, close to that provided by the optimal design (22.51%). Apparently, the presence of a final, diverging channel has greater impact on cell performance than the fine adjustment of channel width at the simulation conditions set herein studied. The design with H5 = 0.1 mm and that with all Wi = 1 mm produce current distributions close to that with optimal design at y < 7 mm (Fig. 12). The current density for the former is low at the far end of the channel. The corresponding sub-rib convection under rib 4 is strong (2.27–3.37 ms−1 ) near the outlet (y > 8 mm) to direct large quantity of oxygen to outlet. Hence, the sub-rib convection at far end of the rib 4 (y < 6 mm) becomes weak to effective remove produced liquid water, hence yielding a sharp decrease in the local current density (Fig. 13). Conversely, the design with all the same channel widths of 1 mm produces a current density distribution close to that with optimal design over the entire cell. This observation reveals that the channel

Fig. 13. Current density distributions under rib 4 and channel 5 in the membrane for various flow fields. (A) Basic design; (B) optimal design; (C) design with H5 = 0.1 mm; and (D) design with all Wi = 1 mm.

height affects sub-rib convection flow more significantly than does the channel width. The optimization of the channel width at 0.4 V applied voltage has a contribution to performance improvement of < 2.3%. 5. Conclusions This work adopted a combined optimization procedure, including a simplified conjugate-gradient method and a completely three-dimensional, two-phase, non-isothermal fuel cell model, to look for optimal flow field design for a single serpentine fuel cell of size 9 mm × 9 mm with five channels. The cell output power density Pcell is maximized subjected to an optimal set of channel heights, H1 –H5 , and channel widths, W2 –W5 . The basic case with all channel heights and widths set at 1 mm yields a Pcell = 7260 W m−2 . The optimal design displays a tapered characteristic for channels 1, 3 and 4, and a diverging characteristic in height for channels 2 and 5, pro-

X.-D. Wang et al. / Electrochimica Acta 54 (2009) 5522–5530

ducing a Pcell = 8894 W m−2 , about 22.5% increment. The reduced channel heights of channels 2–4 significantly increase the sub-rib convection and widths for effectively removing liquid water and oxygen transport in gas diffusion layer. The final diverging channel minimizes the leakage of fuel to outlet via sub-rib convection from channel 4 to channel 5. Near-optimal design without huge loss in cell performance but is easily manufactured is tested. The use of a straight, final channel of 0.1 mm height has led to 7.37% power loss, while the design with all channel widths to be 1 mm with optimal channel heights obtained above yields only 1.68% loss of current density. The presence of a final, diverging channel has greater impact on cell performance than the fine adjustment of channel width at the simulation conditions set herein studied. Acknowledgment This study was supported by the National Natural Science Foundation of China (No. 50636020 and No. 50876009). Appendix A. Nomenclature

5529

membrane dry density (kg m−3 ) proton conductivity (S m−1 ) electron conductivity (S m−1 )

dry m s

subscripts eff effective g gaseous phase H2 O water k kth species of the mixture l liquid phase Appendix B. The PEM fuel cell model includes continuity, momentum and species equations for gaseous species, liquid water transport equations in the channels, gas diffusion layers, and catalyst layers, water transport equation in the membrane, electron and proton transport equations. They include: (i) Continuity equation for the gaseous species:





∇ · εg u g = −SL cp C D D F hfg Hi I j J krl kp M Mm nd  pc S Sc Sj SL Su T  u Vcell X Xi Y Z

specified heat (J kg−1 K−1 ) mass fraction mass diffusivity (m2 s−1 ) water diffusivity in the membrane (m2 s−1 ) Faraday constant (96,487 C/mol) evaporation latent heat of water (J kg−1 ) channel height (m) average current density in the fuel cell (A m−2 ) transfer current density (A m−3 ) object function relative permeability of the liquid water permeability (m2 ) molecular weight (kg mol−1 ) membrane equivalent weight (kg mol−1 ) electro-osmotic drag coefficient pressure (atm) capillary pressure (atm) ratio of the liquid water volume to the pore volume source term in the species equation source term in the phase potential equation source term accounting due to phase change of water source term in the momentum equation cell temperature (K) velocity vector (m s−1 ) operating voltage (V) x-direction coordinate (m) optimized search variable y-direction coordinate (m) z-direction coordinates (m)

Greek letters ˚m ionic phase potential (V) ˚m electronic phase potential (V)  water content in membrane search step size ˇi i conjugate-gradient coefficient ε porosity  search direction eff effective thermal conductivity (W m−1 K−1 ) s thermal conductivity of solid matrix (W m−1 K−1 ) thermal conductivity of fluids in the pores (W m−1 K−1 ) f  viscosity (kg m−1 s−1 )  density (kg m−3 )

(B1)

(ii) Momentum equation for the gaseous species:



ε 2

(1−s)



∇ · g u g u g =−ε∇ pg +

  ε ∇ · g ∇ u g +Su (1 − s)

(B2)

(iii) Species equation for the gaseous species:









∇ · εg u g Ck = ∇ · g Dk,eff ∇ Ck + Sc − SL

(B3)

(iv) Liquid water transport equation in the flow channels, gas diffusion layers and catalyst layers:



∇·

l kp krl ∂pc ∇s l ∂s



+∇ ·

nd MH2 O  im F





−∇ ·



l kp krl ∇ pg l



= SL

(B4)

(v) Liquid water transport equation in the membrane:



∇·

˛d MH2 O  im F





−

MH2 O dry Mm



D



∇

=0

(B5)

(vi) Proton and electron transport equations:

∇ · ( m ∇ ˚m ) = Sj

(B6)

∇ · ( s ∇ ˚s ) = −Sj

(B7)

(vii) The energy equation for the PEM fuel cell:







∇ · ε (1 − s) g u g Cp,g T + ∇ · εsl u l Cp,l T = ∇ · (eff ∇ T ) + j +

i2 + hfg SL

 (B8)

where j is the transfer current density calculated by Bulter-Volmer equation, eff is effective thermal conductivity, accounting for contribution of solid matrix and fluids in the pores of porous media, and can be expressed as: eff = −2s +

1 ε/(2s + f ) + (1 − ε)/3s

(B9)

Appendix C. The conjugate-gradient method evaluates the gradients of the objective function and sets up a conjugate directionfor the

5530

X.-D. Wang et al. / Electrochimica Acta 54 (2009) 5522–5530

updated search variables with the help of a sensitivity analysis. The negative gradient direction of objective function is specified as the first search direction in conjugate-gradient method, that is:



−∇ J =

−

∂F ∂F ∂F ,− , · · ·, − ∂X1 ∂X2 ∂XN



(C1)

where ∂F/∂Xi is referred to as the sensitivity coefficient. The sensitivity coefficient is calculated by introducing a small perturbation (Xi ) to the search variable, Xi . The new search variable, Xi , is updated by: (k+1)

Xi

(k)

= Xi (k)

where Xi

(k) (k+1)

+ ˇi i

(k+1)

and Xi

i = 1, 2, · · ·N

(C2)

denote values of Xi , in the kth and the (k+1)th (k)

search steps, respectively; ˇi

denotes the search step size of Xi in

(k+1) i

the kth search step; denotes the search direction of Xi , which can be expressed as a linear combination of the previous search (k) direction, i and the negative gradient direction of new objective



function, − ∂F/∂Xi (k+1)

i

=−

(k+1)

, that is,

∂F (k+1) (k+1) (k) + i i ∂Xi

(k+1)

where i

(k)

 (k+1)

=

(C3)

is referred to as the conjugate-gradient coefficient,

which must guarantee that i expressed as: i

i = 1, 2, · · ·N

(∂F/∂Xi )

(k+1)

(∂F/∂Xi )



is conjugated to − ∂F/∂Xi

(k+1)

and

2 i = 1, 2, · · ·N

(k)

(C4)

In basic conjugate-gradient method, the optimal search step (k) size, ˇi (i = 1, 2,. . ., N), needs to be determined. After the kth search, the objective function becomes: J (k) = F =F

 

(k)

(k)

(k)

(k−1)

+ ˇ1

X1 , X2 , · · ·XN X1



(k−1) (k) (k−1) 1 , X2

(k−1) (k) (k−1) 2 , · · ·, XN

+ˇ2

(k−1) (k) N

+ ˇN

 (C5)

References [1] T.E. Springer, T.A. Zawodinski, S. Gottesfeld, J. Electrochem. Soc. 136 (1991) 2334. [2] D.M. Bernardi, M.W. Verbrugge, J. Electrochem. Soc. 139 (1992) 2477. [3] T.F. Fuller, J. Newman, J. Electrochem. Soc. 140 (1993) 1218. [4] V. Gurau, F. Barbir, H. Liu, J. Electrochem. Soc. 47 (2000) 2468. [5] N. Djilali, D. Lu, Int. J. Therm. Sci. 41 (2002) 29. [6] D. Singh, D.M. Lu, N. Djilali, Int. J. Eng. Sci. 37 (1999) 431. [7] J.S. Yi, T.V. Nguyen, J. Electrochem. Soc. 146 (1999) 38. [8] X.D. Wang, Y.Y. Duan, W.M. Yan, J. Power Sources 172 (2007) 265. [9] S. Dutta, S. Shimpalee, J.W. Van, Zee, Int. J. Heat Mass Transfer 44 (2001) 2029. [10] S. Um, C.Y. Wang, J. Power Sources 125 (2004) 40. [11] T. Okada, X. Gang, M. Meeg, Electrochim. Acta 43 (1998) 2141. [12] J.J. Baschuk, X. Li, J. Power Sources 86 (2000) 181. [13] S. Mazumder, J.V. Cole, J. Electrochem. Soc. 150 (2003) A1510. [14] J.J. Hwang, J. Power Sources 164 (2007) 174. [15] H. Meng, J. Power Sources 168 (2007) 218. [16] C.Y. Wang, Chem. Rev. 104 (2004) 4727. [17] W.Q. Tao, C.H. Min, X.L. Liu, Y.L. He, B.H. Yin, W. Jiang, J. Power Sources 160 (2006) 359. [18] H. Li, Y. Tang, Z. Wang, Z. Shi, S. Wu, D. Song, J. Zhang, K. Fatih, J. Zhang, H. Wang, Z. Liu, R. Abouatallah, A. Mazza, J. Power Sources 178 (2008) 103. [19] A.C. West, T.F. Fuller, J. Appl. Electrochem. 6 (1996) 557. [20] T.V. Nguyen, J. Electrochem. Soc. 143 (1996) L103. [21] S. Shimpalee, J.W. Van, Zee, Int. J. Hydrogen Energy 32 (2007) 842. [22] X.D. Wang, Y.Y. Duan, W.M. Yan, X.F. Peng, J. Power Sources 175 (2008) 397. [23] A. Kumar, G.R. Ramana, J. Power Sources 113 (2003) 11. [24] S. Shimpalee, S. Greenway, J.W. Van Zee, J. Power Sources 160 (2006) 398. [25] S.W. Cha, R. O’Hayre, Y. Saito, F.B. Prinz, J. Power Sources 134 (2004) 57. [26] A. Turhan, K. Heller, J.S. Brenizer, M.M. Mench, J. Power Sources 160 (2006) 1195. [27] S.S. Hsieh, S.H. Yang, C.L. Feng, J. Power Sources 162 (2006) 262. [28] C. Xu, T.S. Zhao, Electrochem. Commun. 9 (2007) 497. [29] X.D. Wang, Y.Y. Duan, W.M. Yan, J. Power Sources 173 (2007) 210. [30] W.C. Weng, W.M. Yan, H.Y. Li, X.D. Wang, J. Electrochem. Soc. 155 (2008) B877. [31] C.H. Huanga, L.Y. Chena, S. Kimb, J. Power Sources 187 (2009) 136. [32] C.H. Min, J. Power Sources 186 (2009) 370. [33] H.K. Ma, S.H. Huang, Y.Z. Kuo, J. Power Sources 185 (2008) 1154. [34] I. Mohamed, N. Jenkins, J. Power Source 131 (2004) 142. [35] M. Prud’homme, T.H. Nguyen, Int. Commun. Heat Mass Transfer 25 (1998) 999. [36] Y. Hung, N.G. Zamani, Appl. Math. Comput. 80 (1996) 155. [37] P.C. Tuan, M.C. Ju, Numer. Heat Transfer, Part B 37 (2000) 247. [38] C.H. Cheng, M.H. Chang, J. Power Source 139 (2005) 115. [39] M.H. Chang, C.H. Cheng, J. Power Source 142 (2005) 200. [40] C.H. Cheng, M.H. Chang, Numer. Heat Transfer, Part B 43 (2003) 489. [41] X.D. Wang, X.X. Zhang, W.M. Yan, D.J. Lee, A. Su, Int. J. Hydrogen Energy (2009), doi:10.1016/j.ijhydene.2008.12.049. [42] H.K. Versteeg, W. Malasekera, An Introduction to Computational Fluid Dynamics, John Wiley & Sons, Inc., New York, 1995.