- Email: [email protected]
Numerical analysis of a polymer electrolyte fuel cell
Available online at www.sciencedirect.com
International Journal of Hydrogen Energy 29 (2004) 945 – 954
www.elsevier.com/locate/ijhydene
Numerical analysis of a polymer electrolyte fuel cell H.-M. Junga , W.-Y. Leeb;∗ , J.-S. Parka , C.-S. Kimb a School
of Aerospace and Mechanical Engineering, Hankuk Aviation University, 200-1 Hwajon, Duckyang, Kyunggi, South Korea b Fuel Cell Research Center, Korea Institute of Energy Research, P.O. Box 103, Yusong, Daejon 305-343, South Korea Received 1 September 2002; received in revised form 1 January 2003; accepted 1 March 2003
Abstract Numerical simulation was carried out to predict the 5ow, temperature, and current distributions in a polymer electrolyte fuel cell (PEFC). The continuity, momentum, and energy equations with mass and heat source/sink terms produced by chemical reactions are solved using a general computational 5uid dynamic codes. A local current density at each point on the electrode surface is calculated as a function of gas pressure, cell temperature, humidity, and partial pressure as well as cell voltage using an empirical electro-chemical equation. The performance was simulated using Na:on 115 membrane with the active area of 100 cm2 with a serpentine 5ow channel, which is suited for uniform gas supply. The predictions indicate that 5ow distribution and current production are a;ected signi:cantly by each other. This approach can be used to understand and investigate the e;ects of various parameters, such as the pattern and dimensions of 5ow channel con:gurations, operating conditions, such as inlet humidity, reactant utilization ratio, and pressure on the PEFC performance for optimal design. ? 2003 International Association for Hydrogen Energy. Published by Elsevier Ltd. All rights reserved. Keywords: PEFC; Numerical analysis; Computational 5uid dynamics
1. Introduction Fuel cells are considered as alternative devices to generate electric power in transportation and stationary applications. Among di;erent types of fuel cells, PEFC is the most promising type for transportation and small-scale electric power generation applications due to its high-power density, simple and safe construction and fast start-up, even at low temperatures. Recently, considerable attention has been given to computational research related to fuel cell performance analysis and optimization. Many institutes and industries have started to use methods and techniques o;ered by computational 5uid dynamics to understand and optimize electro-chemical processes that take place in fuel cells better. Computational ∗
Corresponding author. Fax: +82-42-860-3739. E-mail address: [email protected] (W.-Y. Lee).
analyses are capable to predict the 5ow, pressure, temperature, and current distributions in a speci:ed fuel cell. However, most of the existing research involving PEFC are limited to a single channel and utilize symmetric assumptions on channel sides. It is necessary to extend to the entire cell or stack to understand practical PEFC performance. In this study, a 100 cm2 single cell is investigated numerically for performance analysis and optimal design [1]. Many parameters are associated with PEFC design. Specially, the voltage and current produced by the electro-chemical reaction in the stack depend on the con:guration of channels in 5ow :eld plate and the operating conditions, such as stack temperature, gas usage, gas partial pressures, and humidity. When a fuel cell operates, thermal energy is produced associated with the production of electricity and released by a cooling plate. The local heating rate at the stack is calculated as a function of the cell voltage and distributed current. The performance of PEFC
0360-3199/$ 30.00 ? 2003 International Association for Hydrogen Energy. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.ijhydene.2003.03.001
946
H.-M. Jung et al. / International Journal of Hydrogen Energy 29 (2004) 945 – 954 current collector bipolar plate gas outlet (manifold) gas diffusion layer
(a)
(b)
(c)
(d)
gas inlet (manifold)
MEA
Fig. 1. Schematic diagram of fuel cell (PEFC).
is a;ected by these parameters. An investigative technique is well suited for this task is numerical modeling. This reduces the time and cost associated with an optimal design [2]. Also, the theory used to develop the modeling technique includes the electrical and electro-chemical characteristics of the membrane electrode assembly (MEA) of the PEFC. Electrical characteristics used when developing the model include the equal potential nature of the electrodes. In order to simplify the calculation of the electro-chemical characteristics of the MEA an empirical equation that represents the relationship between the current and voltage derived from experimental data was used. These characteristics were combined to develop a modeling technique that calculates fuel cell performance. Implementation of the modeling and simulation technique is described and demonstrated by single cell with the active area of 100 cm2 . In this study, a number of di;erent 5ow con:guration were evaluated to :nd the proper 5ow channels and then using the selected 5ow channel, thermal and electro-chemical performances are investigated. In this paper we concentrate on analyzing the 5ow, pressure, temperature and current distributions for better understanding and optimal design of PEFC. This scheme can be used not only to investigate 5ow, thermal and electrical analysis but also to obtain information for the design and operation of PEFC systems (Fig. 1). 2. Numerical method 2.1. Numerical model First of all, several types of 5ow :eld geometry are analyzed to :nd the optimal channel con:gurations. The con:g-
Fig. 2. Several types of 5ow :eld plate geometry, (a) the parallel channels, (b) the discontinuous channels, (c) the serpentine channel, (d) the spiral channel.
urations of the 5ow channels are shown in Fig. 2. Because the e;ect of the channel con:guration is more important on the air side, the investigation on 5ow :eld focused on it. These plates have the same active area with 100 cm2 . Also, the width and height of 5ow channels on their surface is modeled to 1 mm. The e;ect of gas channel con:gurations was investigated with the 2-D model without consideration of chemical reaction. Fig. 3 shows the 3-D domain. A single cell has three 5ow :eld plates. Among several types of 5ow :eld plate geometry, the serpentine channel type to investigate thermal and electro-characteristic is used (Fig. 4). The length of 5ow path on three 5ow :eld plates is equivalent to 10 cm in the axial direction. These models include: (A) Rectangular channels of same width (1 mm) and height (1 mm) with hydrogen gas or water are used. The width of channel is 1:2 mm and height is 1 mm with air gas. The width of rib on three 5ow :eld plates is equivalent to 1 mm. (B) Electro-chemical model is divided into eight sub-regions (Fig. 3(a)): the anode gas channel, the anode gas di;usion layer, the anode catalyst layer, membrane, the cathode catalyst layer, the cathode gas di;usion layer, the cathode gas channel, and the cooling channel. Table 1 lists the details of the stack design and the operating conditions. (C) Porous media block with attached supply plenums. We assumed that 5ow resistance of porous media to be uniform. (D) The MEA model consists of a Na:on 115 proton exchange membrane and two carbon paper electrodes both
H.-M. Jung et al. / International Journal of Hydrogen Energy 29 (2004) 945 – 954
947
The end plate The channels of water plate
The ribs area of air plate
The channels of hydrogen plate
The channels of air plate
The ribs area of hydrogen plate
The cathode gas diffusion layer + catalyst layer
(a)
Nafion 115 membrane The anode gas diffusion layer + catalyst layer Hydrogen gas inlet
Air gas inlet
Cooling water inlet
Cooling water outlet
Air & residual gases outlet
(b) Hydrogen & residual gases outlet
Fig. 3. (a) Layout of the basic geometrical element, (b) grid arrangement used in a single cell (3-D).
impregnated with 0:6 mg Pt=cm2 . The porosity of carbon paper is assumed to 70%. In order to consider the in5uence of the manifold on 5ow :eld, manifolds are included in the domain of calculation. In the numerical analysis, variables are pressure, temperature and three velocities in the Cartesian coordinates. The very re:ned meshes for many degrees of freedom to get reasonable accuracy are used. Initially, 3,634,632 meshes are created and locally more re:ned on the catalyst/cathode dif-
fusion layer. The conservation of mass, momentum and energy equations in the 3-D 5ow solver software, STAR-CD, is modi:ed to include the electro-chemical aspects of a fuel cell [3–5]. 2.2. Governing equations The fuel cell operation is described by continuity and energy conservation principles. The mass conservation
948
H.-M. Jung et al. / International Journal of Hydrogen Energy 29 (2004) 945 – 954
The inlet of hydrogen
The outlet of residual gases
(a)
The inlet of air
The outlet of residual gases
(b)
The outlet of cooling water
The inlet of cooling water
(c) Fig. 4. Flow channel con:guration on 5ow-:eld plates (2-D), (a) the hydrogen 5ow channel con:guration on 5ow-:eld plate, (b) the air 5ow channel con:guration on 5ow-:eld plate, (c) the cooling 5ow channel con:guration on 5ow-:eld plate.
equation is @ (u˜ j ) = sm ; @xj
(1)
where u˜ j is the relative velocity in 5uid local coordinate frame and sm is the mass source term. Even though the source term in most computation domain is neglected, the source term should be considered in PEFC because hydrogen
H.-M. Jung et al. / International Journal of Hydrogen Energy 29 (2004) 945 – 954 Table 1 Physical parameters and properties for the simulation Quantity
Value
The anode di;usion layer width The cathode di;usion layer width Membrane width The anode catalyst layer width The cathode catalyst layer width Volume fraction membrane in catalyst layer Relative humidity of inlet air/fuel Pressure of reactant gases Stochiometry of hydrogen Stochiometry of air Fuel cell temperature Inlet temperature of the cooling water Cell voltage
0:000284 m 0:000284 m 0:000125 m 0:000065 m 0:00005 m 0.5 100%/100% 101; 325 Pa 1.42 2.5 343 K 333 K 0:6 V
where PT is the temperature di;erence between inlet temperature and outlet temperature of the cooling, Pe is electrical power of a single cell and the average voltage of single cell is Vc . For solids and constant density 5uids, such as water, this CFD code solves the transport equation for the speci:c internal energy, e, which is given as e ≡ cT N − c o To + m m Hm ; (6) where cN is the mean constant-volume speci:c heat. The usage of hydrogen and air are calculated by using Eqs. (6) and (7). These equations allow the gas usage of any fuel cell system of given power to be calculated: m˙ air = 3:57 × 10−7 ×
1 Pe × Ua Vc
m˙ hydrogen = 1:05 × 10−8 × and oxygen react and some of them are consumed. The mass source term can be considered as a sink and contains electro-chemical aspects of a fuel cell. More speci:cally, these terms correspond to the consumption of hydrogen and oxygen along the channels. In this paper, the sm term is calculated in the sub-program, which is connected to the main 5ow solver. The heat transfer is implemented through the following general form of the enthalpy conservation equation for a 5uid mixture. It should be noted that the static enthalpy is de:ned as the sum of the thermal and chemical components. In this paper, the static enthalpy is equal to the thermal components: @ @ui @P (u˜ j h − Fh; j ) = u˜ + ij + sh ; @xj @xj @xj h = cNp T − cpo To +
m m H m = ht +
(2) m m Hm ;
(3)
where h is static enthalpy, mm the mass fraction of mixture constituent m, Hm the heat of formation in the constituent m, and sh is energy source. The sh term can be considered as heat source in catalyst/cathode di;usion layer contact surface. The heat source at a single cell is calculated as shown in Eq. (4). This term is calculated in the sub-program: 1:21 qi = vi × ii × −1 ; (4) vi where qi is the heating rate, vi is the local voltage through each element, ii is the local current through each element. The heat value generated by a fuel cell can be calculated using Eq. (4) for deciding on the mass 5ow rate of cooling. If n% of this heat is removed by cooling water, the mass rate of the cooling water is nPe 1:25 m˙ cooling = −1 (kg=s); (5) cp PT Vc
949
Pe Vc
(kg=s); (kg=s);
(7) (8)
where Ua is the fraction of the oxygen usage of any fuel cell system and UH is the fraction of the usage of hydrogen. It is assumed that the usage of hydrogen is 70% and oxygen is 40% in the calculation. The usage of hydrogen means that the mass 5ow rate of hydrogen supplied by hydrogen inlet manifold is consumed 70% of gases, and 20% of residual gases is discharged to hydrogen outlet manifold. Also, the usage of oxygen means that 40% of supplied oxygen gases is consumed by chemical reaction. Oxygen is originated from the air inlet manifold, where 79% of air is nitrogen and 21% of air is oxygen. And 60% of residual gas, which is consist of oxygen and vapor water is discharged to the air outlet manifold. These terms are calculated in the sub-program. For the 5ow of the porous media, the void volume or pore can be regarded as the distributed resistance. The gas permeability equation can be obtained from the force balance between pressure-induced forces and resistance forces in the porous media: − K i ui =
@P ; @ i
(9)
where i is the orthotropic directions, Ki is the permeability and ui is the super:cial velocity in direction i . The permeability Ki is assumed to be a quasi-linear function of the super:cial velocity magnitude |˜v|: Ki = i |˜v| + !i ;
(10)
where i and !i represent permeability coeRcients in porous media. These constant values are determined by numerical experimentation in this calculation. By analogy with the equations of state (density), −1 mm = ; (11) m m where mm and m are mass fraction and density of constituent m.
950
H.-M. Jung et al. / International Journal of Hydrogen Energy 29 (2004) 945 – 954
The overall species conservation equation is mm = 1;
START
(12)
m
where subscript m is the component index. The MEA modeling technique used empirical equations derived experimentally at active area 100 cm2 . The empirical equation that is used for PEFC stack models relates the voltage of a MEA to its current. The other fuel cell relation can be developed to emphasize di;erent operational and design aspect. The following relation is useful to predict the performance numerically on the MEA. Vc = F(Ic ; Tc ; Hca ; Han ; Pca ; Pan ; PH2 ; PO2 ; Uh ; Ua );
1st STEP 1) A local geometric current density at each point along the channel is calculated. 2) The partial pressure of reactants is calculated. Sub-routine Program
CFD (STAR-CD)
initial current density calculated
1) The pressure drop between supplied gases inlet and discharged residual gases are calculated. 2) The temperature profile is calculated. The maximum and minimum stack temp. are used to determine the cooling flowrate.
2.3. Boundary conditions This numerical simulation is assumed to be steady state. A single fuel cell has a potential of about 0:6 V when generating current of approximately 400 mA=cm2 . Then, the operating cell temperature is 70◦ C. The humidity assumes constant value. As a base operational conditions, the mass 5ow rate of hydrogen is 4:2 × 10−7 kg=s and the mass 5ow rate of air is 3:57 × 10−5 kg=s. It is assumed that usage of hydrogen is 70% and air is 40%. Air is mainly composed of nitrogen and oxygen. The cathode-side gas mixture consists of 79% nitrogen and 21% oxygen. This numerical simulation is assumed to be steady state. Also, gases are assumed to be ideal and incompressible. Proper gas distribution and heat management are essential for better energy eRciencies in a PEFC. In order to use the active area e;ectively, the reactants are distributed homogeneously by the 5ow channel and manifold. The gases move from the 5ow channels of the 5ow :eld plate through the porous electrodes to the reaction sites near the membrane. The gases move from an area of greater concentration, the 5ow channel of the 5ow :eld plate, to an area of lesser concentration, the electrode–membrane junction [6,7].
The properties (pressure temperature, current etc.) of each element are assumed to be locally uniform.
2nd STEP
(13)
where Vc is the cell potential, Ic is the cell current, Pca is the cathode side pressure, Pan is the anode side pressure, PH2 is the partial pressure of hydrogen and PO2 is the partial pressure of oxygen. This empirical equation is calculated in the sub-program.
mass sink, heat source
3rd STEP 1) The calculation of current dentity distribution is continued until the initial current density and total current density are met. 2) The usages of oxygen and hydrogen are calculated once more. Sub-routine Program
The pressure and temperature are calculated using STAR-CD. The new current distribution is calculated.
mass sink, heat source
4th STEP CFD (STAR-CD)
The calculations of local mass flow rate and temperature through each element are repeated.
Fig. 5. Diagram of a single fuel cell modeling numerical procedure.
empirical electro-chemical equation. The numerical calculation procedures are shown in Fig. 5.
3. Results and discussions 2.4. Numerical solution procedure A :nite volume method is used to solve the numerical model. The continuity and energy equations with mass and heat source/sink terms produced by chemical reactions are solved using a general computational 5uid dynamic code. A local current density at each point on the electrode surface is calculated as a function of gas pressure, cell temperature, humidity, and partial pressure as well as voltage using the
Results are divided into two parts: the results based on 5ow distributions of four di;erent 5ow :eld con:gurations then the performance of a single cell are discussed. The di;erences for the two results are for the 5ow :eld analysis governing equations related velocity and pressure are used, and for a single cell performance whole governing equations are used including the electro-chemical equation to calculate the 5ow, pressure, temperature and current distributions.
The Divided Rate of Total Mass Flow (%)
H.-M. Jung et al. / International Journal of Hydrogen Energy 29 (2004) 945 – 954 Table 2 The pressure drop of four 5ow channel types
8 7
The channel type
Pressure drop (kPa)
6
The The The The
0.37 2.76 1.6 2.87
5 4
A
3
parallel channels type discontinuous channels type serpentine channel type spiral channel type
2 1 0 10
20
30
40
50
Channel Order The Divided Rate of Total Mass Flow (%)
951
40
B
C
D
20
0 1
2
3
4
5
1
2
3
4
5
1
2
3
4
5
Channel Order Fig. 6. The divided rate of total mass 5ow in each channel; A is parallel channel type, B is discontinuous channel type, C is serpentine channel type and D is spiral channel type.
3.1. Design of ;ow channel on ;ow
droplets. For this reason, BSuchi et al. (1996) recommend a pressure loss of 1–3 kPa in each channel on the cathode side [8]. As seen in Fig. 6, while on the parallel and discontinuous 5ow :eld, 5ow rates are not distributed evenly between the channels, on the spiral and serpentine 5ow :elds, uniform 5ow rates are obtained. The pressure drops of four channel types are described in Table 2. In each case, the pressure drops in air plate are calculated. The parallel 5ow :eld has the advantage of o;ering a lower pressure than the others. But due to many parallel paths, the 5ow distribution between the channels is not even and water removal may be ine;ective. The discontinuous 5ow channel has the same problem as the parallel 5ow :eld. The spiral channel con:guration has the same even 5ow distribution as the serpentine, but due to vertically upward channels, the effective water removal may not be obtained. In this study, the serpentine with :ve parallel 5ow channels is selected as an optimal 5ow :eld con:guration for the further investigation of thermal and current distribution at a single cell. 3.2. Single cell thermal ;ow and electrical characteristic analysis The performance of a PEFC has been analyzed on the basis of the validated 5ow models, including the integration of the energy balance and empirical electro-chemical equations. The model equations are solved at a given value of cell voltage, 0:6 V. A typical PEFC includes a circuit for cooling liquid water, which is taken into account in the simulation. Upon completion of the simulated the pressure, temperature, and current distributions of the channels and MEA were presented. Fig. 7 shows the pressure distribution on the hydrogen and air 5ow plate. As expected, the pressures are decreased with increasing distance along the channels. It is seen that the pressure drops which are proportional to 5ow rates are evenly distributed between the :ve channels. Unlike the 5ow calculation in Fig. 6, electro-chemical and thermal e;ects are included to calculate the pressure distribution in Fig. 7. Regardless of electro-chemical reactions, the 5ow(and/or pressure) uniformity on the 5ow :eld plates can be obtained when using the serpentine 5ow channel con:guration.
952
H.-M. Jung et al. / International Journal of Hydrogen Energy 29 (2004) 945 – 954
101400
Pressure [Pa]
101200
101000
100800
100600
2
10 4
8 6
6
Length [cm]
4
8
Length [cm]
2 10
Fig. 7. The pressure distribution of air 5ow :eld.
The computed distributions of 5ow, temperature, and current density are shown in Figs. 8–10. The velocity magnitude distribution of the air 5ow channel are represented in Fig. 8. Although the values of the velocity magnitude are di;erent depending on the depth of the channels, but even distribution can be seen between channels. The temperature pro:le and current pro:le on the MEA are shown in Figs. 9 and 10, respectively. The average tem-
perature of the MEA is 357 K and the di;erence between the maximum and minimum temperatures is 4 K. It is seen that a relatively uniform temperature can be obtained on the cell. In Fig. 10, the current density is decreased with increasing longitudinal distance. Obviously, the current density at the inlet area of air supply is higher than at the outlet area. This is clearly related to the values of oxygen concentration.
H.-M. Jung et al. / International Journal of Hydrogen Energy 29 (2004) 945 – 954
953
6.4
6.0
5.6
Velocity magnitude [m/s]
6.8
5.2 2 2 4 4
Length [m]
6
6 8
Length [m]
8
Fig. 8. The velocity magnitude of air channel.
4. Conclusions The numerical modeling and simulation is a powerful tool for optimal design by determining the in5uence of di;erent parameters in a PEFC. A scheme for a numerical analysis of a PEFC is presented. In this study, :rst of all, the 5ow and pressure distributions are analyzed by the 2-D model without consideration
of chemical reactions. Four di;erent 5ow :eld con:gurations are evaluated numerically to :nd proper 5ow channel geometry. The 5ow simulations reveal the 5ows of the parallel serpentine and spiral channels are more evenly distributed between channels than the straight parallel channel. Then, with a serpentine channel con:guration, the 3-D numerical analysis for a single cell has been performed to
954
H.-M. Jung et al. / International Journal of Hydrogen Energy 29 (2004) 945 – 954 370
Temperature [K]
360
350
340
330
320
310
300 0.00
0.02
0.04
0.06
0.08
Horizontal distance [m] Fig. 9. Temperature distribution along normal section of MEA.
Fig. 10. The current distribution for the PEFC.
compute the pressure, velocity, temperature, and current distribution using the mass, momentum, energy equations, and empirical electro-chemical equation.
The energy equation with heat transfer relations are used to calculate temperature distribution and a MEA electro-chemical relationship derived from experimental data is used to calculate current distribution on the cell as a function of operating conditions such as oxygen partial pressure, cell temperature. The simulation results can be used to evaluate the performance of a cell or stacks by the current distribution, and used to design and modify the gas and cooling 5ow plate for the better performance. The model proposed in this study can be used for cell performance analysis at various conditions, such as di;erent cell temperature, reactant pressures, and gas utilization ratio. The scheme would be helpful for investigating fuel cells and for optimal design fuel cell performance. References [1] Vladimir K, Rupak D. Three-dimensional modeling of a medium size PEM fuel cell stack: thermal e;ects and electrical performance. Fourth International ASME/JSME/ KSME Symposium on Computational Technologies for Fluid/Thermal/Chemical Systems with Industrial Applications; August 4–8, Vancover, BS, Canada, 2002. [2] Lee JH, Lark TR, Appleby AJ. Modeling electro-chemical performance in large scale proton exchange membrane fuel cell stacks. J Power Source 1998;70:258–68. [3] Tuomas M. Design and experimental characterization of polymer electrolyte membrane and fuel cells. Helsinki University of Technology; 2000. [4] Lee JH, Lark TR. Modeling fuel cell stack systems. J Power Source 1998;73:229–41. [5] Sandip D, Sirivatch S, Van Zee JW. Numerical prediction of mass-exchange between cathode and anode channels in a PEM fuel cell. J Heat Mass Transfer 2001;44:2029–42. [6] James L, Andrew D. Fuel cell systems explained. John Wiley and Sons Inc. 2000. [7] Um S, Wang C-Y, Chen KS. Computational 5uid dynamics modeling of proton exchange membrane fuel cells. J Electrochem Soc 2000;147(12):4485–93. [8] BSuchi FN, Tsukada A, Haas O, Scherer GG. Polymer electrolyte fuel cells: 5ow :eld for eRcient air operation. Paul Scherer Institute Annual Report. 1996.
International Journal of Hydrogen Energy 29 (2004) 945 – 954
www.elsevier.com/locate/ijhydene
Numerical analysis of a polymer electrolyte fuel cell H.-M. Junga , W.-Y. Leeb;∗ , J.-S. Parka , C.-S. Kimb a School
of Aerospace and Mechanical Engineering, Hankuk Aviation University, 200-1 Hwajon, Duckyang, Kyunggi, South Korea b Fuel Cell Research Center, Korea Institute of Energy Research, P.O. Box 103, Yusong, Daejon 305-343, South Korea Received 1 September 2002; received in revised form 1 January 2003; accepted 1 March 2003
Abstract Numerical simulation was carried out to predict the 5ow, temperature, and current distributions in a polymer electrolyte fuel cell (PEFC). The continuity, momentum, and energy equations with mass and heat source/sink terms produced by chemical reactions are solved using a general computational 5uid dynamic codes. A local current density at each point on the electrode surface is calculated as a function of gas pressure, cell temperature, humidity, and partial pressure as well as cell voltage using an empirical electro-chemical equation. The performance was simulated using Na:on 115 membrane with the active area of 100 cm2 with a serpentine 5ow channel, which is suited for uniform gas supply. The predictions indicate that 5ow distribution and current production are a;ected signi:cantly by each other. This approach can be used to understand and investigate the e;ects of various parameters, such as the pattern and dimensions of 5ow channel con:gurations, operating conditions, such as inlet humidity, reactant utilization ratio, and pressure on the PEFC performance for optimal design. ? 2003 International Association for Hydrogen Energy. Published by Elsevier Ltd. All rights reserved. Keywords: PEFC; Numerical analysis; Computational 5uid dynamics
1. Introduction Fuel cells are considered as alternative devices to generate electric power in transportation and stationary applications. Among di;erent types of fuel cells, PEFC is the most promising type for transportation and small-scale electric power generation applications due to its high-power density, simple and safe construction and fast start-up, even at low temperatures. Recently, considerable attention has been given to computational research related to fuel cell performance analysis and optimization. Many institutes and industries have started to use methods and techniques o;ered by computational 5uid dynamics to understand and optimize electro-chemical processes that take place in fuel cells better. Computational ∗
Corresponding author. Fax: +82-42-860-3739. E-mail address: [email protected] (W.-Y. Lee).
analyses are capable to predict the 5ow, pressure, temperature, and current distributions in a speci:ed fuel cell. However, most of the existing research involving PEFC are limited to a single channel and utilize symmetric assumptions on channel sides. It is necessary to extend to the entire cell or stack to understand practical PEFC performance. In this study, a 100 cm2 single cell is investigated numerically for performance analysis and optimal design [1]. Many parameters are associated with PEFC design. Specially, the voltage and current produced by the electro-chemical reaction in the stack depend on the con:guration of channels in 5ow :eld plate and the operating conditions, such as stack temperature, gas usage, gas partial pressures, and humidity. When a fuel cell operates, thermal energy is produced associated with the production of electricity and released by a cooling plate. The local heating rate at the stack is calculated as a function of the cell voltage and distributed current. The performance of PEFC
0360-3199/$ 30.00 ? 2003 International Association for Hydrogen Energy. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.ijhydene.2003.03.001
946
H.-M. Jung et al. / International Journal of Hydrogen Energy 29 (2004) 945 – 954 current collector bipolar plate gas outlet (manifold) gas diffusion layer
(a)
(b)
(c)
(d)
gas inlet (manifold)
MEA
Fig. 1. Schematic diagram of fuel cell (PEFC).
is a;ected by these parameters. An investigative technique is well suited for this task is numerical modeling. This reduces the time and cost associated with an optimal design [2]. Also, the theory used to develop the modeling technique includes the electrical and electro-chemical characteristics of the membrane electrode assembly (MEA) of the PEFC. Electrical characteristics used when developing the model include the equal potential nature of the electrodes. In order to simplify the calculation of the electro-chemical characteristics of the MEA an empirical equation that represents the relationship between the current and voltage derived from experimental data was used. These characteristics were combined to develop a modeling technique that calculates fuel cell performance. Implementation of the modeling and simulation technique is described and demonstrated by single cell with the active area of 100 cm2 . In this study, a number of di;erent 5ow con:guration were evaluated to :nd the proper 5ow channels and then using the selected 5ow channel, thermal and electro-chemical performances are investigated. In this paper we concentrate on analyzing the 5ow, pressure, temperature and current distributions for better understanding and optimal design of PEFC. This scheme can be used not only to investigate 5ow, thermal and electrical analysis but also to obtain information for the design and operation of PEFC systems (Fig. 1). 2. Numerical method 2.1. Numerical model First of all, several types of 5ow :eld geometry are analyzed to :nd the optimal channel con:gurations. The con:g-
Fig. 2. Several types of 5ow :eld plate geometry, (a) the parallel channels, (b) the discontinuous channels, (c) the serpentine channel, (d) the spiral channel.
urations of the 5ow channels are shown in Fig. 2. Because the e;ect of the channel con:guration is more important on the air side, the investigation on 5ow :eld focused on it. These plates have the same active area with 100 cm2 . Also, the width and height of 5ow channels on their surface is modeled to 1 mm. The e;ect of gas channel con:gurations was investigated with the 2-D model without consideration of chemical reaction. Fig. 3 shows the 3-D domain. A single cell has three 5ow :eld plates. Among several types of 5ow :eld plate geometry, the serpentine channel type to investigate thermal and electro-characteristic is used (Fig. 4). The length of 5ow path on three 5ow :eld plates is equivalent to 10 cm in the axial direction. These models include: (A) Rectangular channels of same width (1 mm) and height (1 mm) with hydrogen gas or water are used. The width of channel is 1:2 mm and height is 1 mm with air gas. The width of rib on three 5ow :eld plates is equivalent to 1 mm. (B) Electro-chemical model is divided into eight sub-regions (Fig. 3(a)): the anode gas channel, the anode gas di;usion layer, the anode catalyst layer, membrane, the cathode catalyst layer, the cathode gas di;usion layer, the cathode gas channel, and the cooling channel. Table 1 lists the details of the stack design and the operating conditions. (C) Porous media block with attached supply plenums. We assumed that 5ow resistance of porous media to be uniform. (D) The MEA model consists of a Na:on 115 proton exchange membrane and two carbon paper electrodes both
H.-M. Jung et al. / International Journal of Hydrogen Energy 29 (2004) 945 – 954
947
The end plate The channels of water plate
The ribs area of air plate
The channels of hydrogen plate
The channels of air plate
The ribs area of hydrogen plate
The cathode gas diffusion layer + catalyst layer
(a)
Nafion 115 membrane The anode gas diffusion layer + catalyst layer Hydrogen gas inlet
Air gas inlet
Cooling water inlet
Cooling water outlet
Air & residual gases outlet
(b) Hydrogen & residual gases outlet
Fig. 3. (a) Layout of the basic geometrical element, (b) grid arrangement used in a single cell (3-D).
impregnated with 0:6 mg Pt=cm2 . The porosity of carbon paper is assumed to 70%. In order to consider the in5uence of the manifold on 5ow :eld, manifolds are included in the domain of calculation. In the numerical analysis, variables are pressure, temperature and three velocities in the Cartesian coordinates. The very re:ned meshes for many degrees of freedom to get reasonable accuracy are used. Initially, 3,634,632 meshes are created and locally more re:ned on the catalyst/cathode dif-
fusion layer. The conservation of mass, momentum and energy equations in the 3-D 5ow solver software, STAR-CD, is modi:ed to include the electro-chemical aspects of a fuel cell [3–5]. 2.2. Governing equations The fuel cell operation is described by continuity and energy conservation principles. The mass conservation
948
H.-M. Jung et al. / International Journal of Hydrogen Energy 29 (2004) 945 – 954
The inlet of hydrogen
The outlet of residual gases
(a)
The inlet of air
The outlet of residual gases
(b)
The outlet of cooling water
The inlet of cooling water
(c) Fig. 4. Flow channel con:guration on 5ow-:eld plates (2-D), (a) the hydrogen 5ow channel con:guration on 5ow-:eld plate, (b) the air 5ow channel con:guration on 5ow-:eld plate, (c) the cooling 5ow channel con:guration on 5ow-:eld plate.
equation is @ (u˜ j ) = sm ; @xj
(1)
where u˜ j is the relative velocity in 5uid local coordinate frame and sm is the mass source term. Even though the source term in most computation domain is neglected, the source term should be considered in PEFC because hydrogen
H.-M. Jung et al. / International Journal of Hydrogen Energy 29 (2004) 945 – 954 Table 1 Physical parameters and properties for the simulation Quantity
Value
The anode di;usion layer width The cathode di;usion layer width Membrane width The anode catalyst layer width The cathode catalyst layer width Volume fraction membrane in catalyst layer Relative humidity of inlet air/fuel Pressure of reactant gases Stochiometry of hydrogen Stochiometry of air Fuel cell temperature Inlet temperature of the cooling water Cell voltage
0:000284 m 0:000284 m 0:000125 m 0:000065 m 0:00005 m 0.5 100%/100% 101; 325 Pa 1.42 2.5 343 K 333 K 0:6 V
where PT is the temperature di;erence between inlet temperature and outlet temperature of the cooling, Pe is electrical power of a single cell and the average voltage of single cell is Vc . For solids and constant density 5uids, such as water, this CFD code solves the transport equation for the speci:c internal energy, e, which is given as e ≡ cT N − c o To + m m Hm ; (6) where cN is the mean constant-volume speci:c heat. The usage of hydrogen and air are calculated by using Eqs. (6) and (7). These equations allow the gas usage of any fuel cell system of given power to be calculated: m˙ air = 3:57 × 10−7 ×
1 Pe × Ua Vc
m˙ hydrogen = 1:05 × 10−8 × and oxygen react and some of them are consumed. The mass source term can be considered as a sink and contains electro-chemical aspects of a fuel cell. More speci:cally, these terms correspond to the consumption of hydrogen and oxygen along the channels. In this paper, the sm term is calculated in the sub-program, which is connected to the main 5ow solver. The heat transfer is implemented through the following general form of the enthalpy conservation equation for a 5uid mixture. It should be noted that the static enthalpy is de:ned as the sum of the thermal and chemical components. In this paper, the static enthalpy is equal to the thermal components: @ @ui @P (u˜ j h − Fh; j ) = u˜ + ij + sh ; @xj @xj @xj h = cNp T − cpo To +
m m H m = ht +
(2) m m Hm ;
(3)
where h is static enthalpy, mm the mass fraction of mixture constituent m, Hm the heat of formation in the constituent m, and sh is energy source. The sh term can be considered as heat source in catalyst/cathode di;usion layer contact surface. The heat source at a single cell is calculated as shown in Eq. (4). This term is calculated in the sub-program: 1:21 qi = vi × ii × −1 ; (4) vi where qi is the heating rate, vi is the local voltage through each element, ii is the local current through each element. The heat value generated by a fuel cell can be calculated using Eq. (4) for deciding on the mass 5ow rate of cooling. If n% of this heat is removed by cooling water, the mass rate of the cooling water is nPe 1:25 m˙ cooling = −1 (kg=s); (5) cp PT Vc
949
Pe Vc
(kg=s); (kg=s);
(7) (8)
where Ua is the fraction of the oxygen usage of any fuel cell system and UH is the fraction of the usage of hydrogen. It is assumed that the usage of hydrogen is 70% and oxygen is 40% in the calculation. The usage of hydrogen means that the mass 5ow rate of hydrogen supplied by hydrogen inlet manifold is consumed 70% of gases, and 20% of residual gases is discharged to hydrogen outlet manifold. Also, the usage of oxygen means that 40% of supplied oxygen gases is consumed by chemical reaction. Oxygen is originated from the air inlet manifold, where 79% of air is nitrogen and 21% of air is oxygen. And 60% of residual gas, which is consist of oxygen and vapor water is discharged to the air outlet manifold. These terms are calculated in the sub-program. For the 5ow of the porous media, the void volume or pore can be regarded as the distributed resistance. The gas permeability equation can be obtained from the force balance between pressure-induced forces and resistance forces in the porous media: − K i ui =
@P ; @ i
(9)
where i is the orthotropic directions, Ki is the permeability and ui is the super:cial velocity in direction i . The permeability Ki is assumed to be a quasi-linear function of the super:cial velocity magnitude |˜v|: Ki = i |˜v| + !i ;
(10)
where i and !i represent permeability coeRcients in porous media. These constant values are determined by numerical experimentation in this calculation. By analogy with the equations of state (density), −1 mm = ; (11) m m where mm and m are mass fraction and density of constituent m.
950
H.-M. Jung et al. / International Journal of Hydrogen Energy 29 (2004) 945 – 954
The overall species conservation equation is mm = 1;
START
(12)
m
where subscript m is the component index. The MEA modeling technique used empirical equations derived experimentally at active area 100 cm2 . The empirical equation that is used for PEFC stack models relates the voltage of a MEA to its current. The other fuel cell relation can be developed to emphasize di;erent operational and design aspect. The following relation is useful to predict the performance numerically on the MEA. Vc = F(Ic ; Tc ; Hca ; Han ; Pca ; Pan ; PH2 ; PO2 ; Uh ; Ua );
1st STEP 1) A local geometric current density at each point along the channel is calculated. 2) The partial pressure of reactants is calculated. Sub-routine Program
CFD (STAR-CD)
initial current density calculated
1) The pressure drop between supplied gases inlet and discharged residual gases are calculated. 2) The temperature profile is calculated. The maximum and minimum stack temp. are used to determine the cooling flowrate.
2.3. Boundary conditions This numerical simulation is assumed to be steady state. A single fuel cell has a potential of about 0:6 V when generating current of approximately 400 mA=cm2 . Then, the operating cell temperature is 70◦ C. The humidity assumes constant value. As a base operational conditions, the mass 5ow rate of hydrogen is 4:2 × 10−7 kg=s and the mass 5ow rate of air is 3:57 × 10−5 kg=s. It is assumed that usage of hydrogen is 70% and air is 40%. Air is mainly composed of nitrogen and oxygen. The cathode-side gas mixture consists of 79% nitrogen and 21% oxygen. This numerical simulation is assumed to be steady state. Also, gases are assumed to be ideal and incompressible. Proper gas distribution and heat management are essential for better energy eRciencies in a PEFC. In order to use the active area e;ectively, the reactants are distributed homogeneously by the 5ow channel and manifold. The gases move from the 5ow channels of the 5ow :eld plate through the porous electrodes to the reaction sites near the membrane. The gases move from an area of greater concentration, the 5ow channel of the 5ow :eld plate, to an area of lesser concentration, the electrode–membrane junction [6,7].
The properties (pressure temperature, current etc.) of each element are assumed to be locally uniform.
2nd STEP
(13)
where Vc is the cell potential, Ic is the cell current, Pca is the cathode side pressure, Pan is the anode side pressure, PH2 is the partial pressure of hydrogen and PO2 is the partial pressure of oxygen. This empirical equation is calculated in the sub-program.
mass sink, heat source
3rd STEP 1) The calculation of current dentity distribution is continued until the initial current density and total current density are met. 2) The usages of oxygen and hydrogen are calculated once more. Sub-routine Program
The pressure and temperature are calculated using STAR-CD. The new current distribution is calculated.
mass sink, heat source
4th STEP CFD (STAR-CD)
The calculations of local mass flow rate and temperature through each element are repeated.
Fig. 5. Diagram of a single fuel cell modeling numerical procedure.
empirical electro-chemical equation. The numerical calculation procedures are shown in Fig. 5.
3. Results and discussions 2.4. Numerical solution procedure A :nite volume method is used to solve the numerical model. The continuity and energy equations with mass and heat source/sink terms produced by chemical reactions are solved using a general computational 5uid dynamic code. A local current density at each point on the electrode surface is calculated as a function of gas pressure, cell temperature, humidity, and partial pressure as well as voltage using the
Results are divided into two parts: the results based on 5ow distributions of four di;erent 5ow :eld con:gurations then the performance of a single cell are discussed. The di;erences for the two results are for the 5ow :eld analysis governing equations related velocity and pressure are used, and for a single cell performance whole governing equations are used including the electro-chemical equation to calculate the 5ow, pressure, temperature and current distributions.
The Divided Rate of Total Mass Flow (%)
H.-M. Jung et al. / International Journal of Hydrogen Energy 29 (2004) 945 – 954 Table 2 The pressure drop of four 5ow channel types
8 7
The channel type
Pressure drop (kPa)
6
The The The The
0.37 2.76 1.6 2.87
5 4
A
3
parallel channels type discontinuous channels type serpentine channel type spiral channel type
2 1 0 10
20
30
40
50
Channel Order The Divided Rate of Total Mass Flow (%)
951
40
B
C
D
20
0 1
2
3
4
5
1
2
3
4
5
1
2
3
4
5
Channel Order Fig. 6. The divided rate of total mass 5ow in each channel; A is parallel channel type, B is discontinuous channel type, C is serpentine channel type and D is spiral channel type.
3.1. Design of ;ow channel on ;ow
droplets. For this reason, BSuchi et al. (1996) recommend a pressure loss of 1–3 kPa in each channel on the cathode side [8]. As seen in Fig. 6, while on the parallel and discontinuous 5ow :eld, 5ow rates are not distributed evenly between the channels, on the spiral and serpentine 5ow :elds, uniform 5ow rates are obtained. The pressure drops of four channel types are described in Table 2. In each case, the pressure drops in air plate are calculated. The parallel 5ow :eld has the advantage of o;ering a lower pressure than the others. But due to many parallel paths, the 5ow distribution between the channels is not even and water removal may be ine;ective. The discontinuous 5ow channel has the same problem as the parallel 5ow :eld. The spiral channel con:guration has the same even 5ow distribution as the serpentine, but due to vertically upward channels, the effective water removal may not be obtained. In this study, the serpentine with :ve parallel 5ow channels is selected as an optimal 5ow :eld con:guration for the further investigation of thermal and current distribution at a single cell. 3.2. Single cell thermal ;ow and electrical characteristic analysis The performance of a PEFC has been analyzed on the basis of the validated 5ow models, including the integration of the energy balance and empirical electro-chemical equations. The model equations are solved at a given value of cell voltage, 0:6 V. A typical PEFC includes a circuit for cooling liquid water, which is taken into account in the simulation. Upon completion of the simulated the pressure, temperature, and current distributions of the channels and MEA were presented. Fig. 7 shows the pressure distribution on the hydrogen and air 5ow plate. As expected, the pressures are decreased with increasing distance along the channels. It is seen that the pressure drops which are proportional to 5ow rates are evenly distributed between the :ve channels. Unlike the 5ow calculation in Fig. 6, electro-chemical and thermal e;ects are included to calculate the pressure distribution in Fig. 7. Regardless of electro-chemical reactions, the 5ow(and/or pressure) uniformity on the 5ow :eld plates can be obtained when using the serpentine 5ow channel con:guration.
952
H.-M. Jung et al. / International Journal of Hydrogen Energy 29 (2004) 945 – 954
101400
Pressure [Pa]
101200
101000
100800
100600
2
10 4
8 6
6
Length [cm]
4
8
Length [cm]
2 10
Fig. 7. The pressure distribution of air 5ow :eld.
The computed distributions of 5ow, temperature, and current density are shown in Figs. 8–10. The velocity magnitude distribution of the air 5ow channel are represented in Fig. 8. Although the values of the velocity magnitude are di;erent depending on the depth of the channels, but even distribution can be seen between channels. The temperature pro:le and current pro:le on the MEA are shown in Figs. 9 and 10, respectively. The average tem-
perature of the MEA is 357 K and the di;erence between the maximum and minimum temperatures is 4 K. It is seen that a relatively uniform temperature can be obtained on the cell. In Fig. 10, the current density is decreased with increasing longitudinal distance. Obviously, the current density at the inlet area of air supply is higher than at the outlet area. This is clearly related to the values of oxygen concentration.
H.-M. Jung et al. / International Journal of Hydrogen Energy 29 (2004) 945 – 954
953
6.4
6.0
5.6
Velocity magnitude [m/s]
6.8
5.2 2 2 4 4
Length [m]
6
6 8
Length [m]
8
Fig. 8. The velocity magnitude of air channel.
4. Conclusions The numerical modeling and simulation is a powerful tool for optimal design by determining the in5uence of di;erent parameters in a PEFC. A scheme for a numerical analysis of a PEFC is presented. In this study, :rst of all, the 5ow and pressure distributions are analyzed by the 2-D model without consideration
of chemical reactions. Four di;erent 5ow :eld con:gurations are evaluated numerically to :nd proper 5ow channel geometry. The 5ow simulations reveal the 5ows of the parallel serpentine and spiral channels are more evenly distributed between channels than the straight parallel channel. Then, with a serpentine channel con:guration, the 3-D numerical analysis for a single cell has been performed to
954
H.-M. Jung et al. / International Journal of Hydrogen Energy 29 (2004) 945 – 954 370
Temperature [K]
360
350
340
330
320
310
300 0.00
0.02
0.04
0.06
0.08
Horizontal distance [m] Fig. 9. Temperature distribution along normal section of MEA.
Fig. 10. The current distribution for the PEFC.
compute the pressure, velocity, temperature, and current distribution using the mass, momentum, energy equations, and empirical electro-chemical equation.
The energy equation with heat transfer relations are used to calculate temperature distribution and a MEA electro-chemical relationship derived from experimental data is used to calculate current distribution on the cell as a function of operating conditions such as oxygen partial pressure, cell temperature. The simulation results can be used to evaluate the performance of a cell or stacks by the current distribution, and used to design and modify the gas and cooling 5ow plate for the better performance. The model proposed in this study can be used for cell performance analysis at various conditions, such as di;erent cell temperature, reactant pressures, and gas utilization ratio. The scheme would be helpful for investigating fuel cells and for optimal design fuel cell performance. References [1] Vladimir K, Rupak D. Three-dimensional modeling of a medium size PEM fuel cell stack: thermal e;ects and electrical performance. Fourth International ASME/JSME/ KSME Symposium on Computational Technologies for Fluid/Thermal/Chemical Systems with Industrial Applications; August 4–8, Vancover, BS, Canada, 2002. [2] Lee JH, Lark TR, Appleby AJ. Modeling electro-chemical performance in large scale proton exchange membrane fuel cell stacks. J Power Source 1998;70:258–68. [3] Tuomas M. Design and experimental characterization of polymer electrolyte membrane and fuel cells. Helsinki University of Technology; 2000. [4] Lee JH, Lark TR. Modeling fuel cell stack systems. J Power Source 1998;73:229–41. [5] Sandip D, Sirivatch S, Van Zee JW. Numerical prediction of mass-exchange between cathode and anode channels in a PEM fuel cell. J Heat Mass Transfer 2001;44:2029–42. [6] James L, Andrew D. Fuel cell systems explained. John Wiley and Sons Inc. 2000. [7] Um S, Wang C-Y, Chen KS. Computational 5uid dynamics modeling of proton exchange membrane fuel cells. J Electrochem Soc 2000;147(12):4485–93. [8] BSuchi FN, Tsukada A, Haas O, Scherer GG. Polymer electrolyte fuel cells: 5ow :eld for eRcient air operation. Paul Scherer Institute Annual Report. 1996.