Chapter 2 Simulation Model
2.4 Boundary Conditions
The governing equations for the current PEMFC model are elliptic, partial differential equations. Hence, boundary conditions are required for all boundaries in the computational domain.
The boundary conditions are presented as follows:
Anodic inlet of the channel:
u=uin, T=Toperating , Y= YH2+ YH2O=1, j=0 (2-38) Anodic outlet of the channel:
P=Patm, j=0 (2-39)
Cathodic inlet of the channel:
u=uin, T=Toperating , Y= YO2+YH2O=1, j=0 (2-40) Cathodic outlet of the channel:
P=Patm, j=0 (2-41)
Outer surface of bipolar plate at Anode:
T=Toperating, V=0 (2-42)
Outer surface of bipolar plate at Cathode:
T=Toperating, V= ηtot=Vre+Vcell (2-43) where ηtot is the total overpotential and Vre is the reversible voltage given by reference [33]:
Vre=0.2329+0.0025T (2-44)
Walls:
T=Toperating, j=0 (2-45)
(a)
(b) (c) Fig. 2.1 Flow field pattern in PEMFC (a) 45-135; (b) 60-120; (c) 90-90,
respectively.
Z X
Y
(a)
(b) (c) Fig. 2.2 Flow field pattern in PEMFC (a) 1007; (b) 1010; (c) 1020, respectively.
Z X
Y
Fig. 2.3 Schematic of electrochemical reaction and the transfer current within a porous catalyst-containing conductor [33].
Fig. 2. 4 Schematic of the PEMFC mass transport model [3].
uin
conv
JC diff
JG
A B
Y GDL
Flow channel
HC
HG
Z
Catalyst layer
Chapter 3
Numerical Methods
3.1 Introduction to CFD-ACE+ software
CFD-ACE+ is a set of computer applications for multi-physics computational analyses. The programs provide integrated geometry and grid generation software, a graphical user interface for preparing the model, a computational solver for performing the simulation, and an interactive visualization software for examining and analyzing the simulation results.
Based on the finite-volume approach [35], governing equations mentioned in chapter 2 are solved by discretization methods which will be described in the following sections.
3.2 Numerical Method for CFD-ACE+
The CFD-ACE+ employs the finite-volume method to discretize the partial differential equations and utilize SIMPLEC scheme to obtain the pressure and velocity fields by solving mass and momentum conservation equations.
After that, use the pressure and velocity fields to substitute into the rest governing equations to have results calculated in sequence.
3.2.1 Finite-volume method
The geometric center of the control volume denoted by P is referred to as the cell center. CFD-ACE+ employs a co-located cell-center variable arrangement that means that all dependent variables and material properties are stored at the cell center P. In other words, the average value of any quantity
within a control volume is given by its cell-center value.
The conservation equations illustrated in chapter 2 can be consisted of four terms which are the unsteady, the convection, the diffusion, and the source terms, respectively. The form of a generalized transport equation is expressed as
+ φ vector, D is the diffusion coefficient, and Sφ is the source term which contains all terms that can’t be included in the previous terms.
Eq. (3-1) is also known as the generic conservation equation for quantity φ.
Integrating this equation over a control-volume cell, yields Eq. (3-2)
∀
where ∀ represents the volume. According to the Divergence theorem of Gauss, Eq. (3-2) can be further written as
∀ where A represents each surface of the control volume P.
The discretized equation can be expressed as
∑
With assembly of transient, convection, diffusion and source terms from numerical integration, governing equations are resulted in the following linear equation.
where φ*P is the previous value from last iteration, the subscript P denotes the value of dependent variable that represent the control volume P, I is inertial relaxation factor, and anb are known as the link coefficients.
In general, both SP and SU represent the source term linearized and can be function of φ. They are evaluated using the latest available value of φ, which is the value at the end of previous iteration. The neighboring-cell term is dominated by the kind of spatial difference scheme.
The integral conservation Eq. (3-2) is applied to each control volume contained in the computational domain, so is Eq. (3-5). Eq. (3-5), the finite difference equation (FDE), is the discrete equivalent of the continuous flow transport equation we discussed before. And the FDE is considered linear by evaluating the anb with the values of φ from the end of previous iteration.
3.2.2 SIMPLEC scheme
Solutions of the three momentum equations yield the three Cartesian components of velocity. Even though pressure is an important flow variable.
No governing PDE for pressure is presented. Pressure-based methods utilize the continuity equation to formulate an equation for pressure, and that’s what SIMPLEC is for.
SIMPLEC stands for Semi-Implicit Method for Pressure-Linked Equations Consistent, and is an enhancement to the will-known SIMPLE algorithm. In SIMPLEC, originally proposed by Van Doormal and Raithby [36], an equation for pressure-correction is derived from the continuity equation.
The finite-difference form of the x-momentum equation can be written as
e e x,e P
where e denotes the value is dependent on the flow direction. The pressure field should be provided to solve for u. If the preceding equation is solved with a guessed pressure P*, it will yield x-component velocity u* which satisfies the following equation
Nevertheless, u* won’t generally satisfy continuity. The strategy is to find corrections to u* and P* so that an improved solution can be obtained. Consider
For pursuing incrementally more accurate solutions, approximation for u'P
is given as
where dp is expressed as
∑
Moreover, an expression for u'e, the velocity correction at the face, is obtained by averaging from the cell-centre values and of the form
' neighboring-cell value. Besides, the continuity equation is written in discretization as
where Ce is the mass flux and written as correction form
'
Therefore, Eq. (3-13) can be recast as:
∑
=where Sm represents the mass source in a control volume:
∑
From Eq. (3-14), the correction term of mass flux is written as
e
In consideration of incompressible flow such as current situation, ρ’ is zero which means that no correction term is referred to density. In the case that the ideal gas assumption is considered, an equation of state is used to estimate the first-order partial derivative of density with respect to pressure, and is written as followings:
The correction to face-normal velocity is obtained as:
'
With the face-normal velocity correction and the density correction thus expressed in terms of pressure correction, a pressure correction equation can be obtained from Eq. (3-15):
∑
+The sequence of operation for SIMPLEC scheme is summarized as follows:
1. Guess the pressure field p*.
2. Solve the discretized momentum equations to obtain u*, v* and w*. 3. Evaluate C* from ρ* and u*, using the procedure outlined before.
4. Evaluate mass source term.
5. Obtain p’ by solving Eq. (3-20).
6. Use p’ to correct the pressure and velocity fields using Eq. (3-8) and (3-9).
7. Solve the discretized equations for other flow variables such as enthalpy, species etc.
8. Repeat the procedure from step 2 until convergence is acquired.
3.3 Computational procedure of PEMFC simulation
The model presented in this investigation was starting at the construction of structured grids by pre-processing software CFD-GEOM from CFD-ACE+.
Fig. 3-1 shows the model with structured grids. After the grids were built, the model was submitted to the solver, CFD-ACE-GUI, for properties and conditions assignments, simulation running included. The model was sequentially assigned to post-processing software called CFD-View for results visualization and extraction.
CFD-ACE+ uses an iterative, segregated solution method wherein the
equation sets for each variable are solved sequentially and repeatedly until a converged solution is obtained. The overall procedure is shown in Fig. 3-2.
At the each iteration, the program calculates residual for each variable, which is the sum of the absolute value of the residual for that variable at each computational node. And the calculation is regarded to be converged as each variable’s normalized residual is less than 10-4 [28].
Numerical calculations are performed on a Pentium Core™2 2.67GHz PC with 4G RAM. After a series of grid independent tests, whose results are listed in Table 3.1, the 729,375 grid system is adopted with a 0.878% error. Detailed information of grid for different components is listed in Table 3.2. In addition, for model validation, based on the reported operating conditions, Fig. 3.3 shows a good agreement between the predicted data of present model and the experimental results from Santaralli and Torchio [37].
Table 3.1 Grid independence test Grid Operating Voltage Current density
(A/cm2) Relative error
641,850 0.4V 631.6254414 2.886%
729,375 0.4V 618.9385131 0.878%
1,503,534 0.4V 613.3949275
Table 3.2 Grid for different components Grid
Components Bipolar
Plates Channels GDLS Catalyst
Layers Membrane 641,850 116,700 155,600 194,500 58,350 116,700 729,375 136,150 175,050 213,950 77,800 126,425 1,503,534 241,948 345,640 414,768 207,384 293,794
Fig. 3.1 Grid structure of current model.
Fig. 3.2 Numerical flow chart of current study [33].
0 100 200 300 400 500 0.4
0.5 0.6 0.7 0.8 0.9 1.0 1.1
Voltage [V]
Current Density (mA/cm2) Current model
M. G. Santarelli & M. F. Torchio [36]
Fig. 3.3 Comparison of current model result with Santarelli and Torchio at 323K.
Chapter 4
Results and Discussion
In this chapter, the numerical results are presented. The effects of bend angle and width on the PEMFC performance at different operating temperatures are analyzed and discussed with the measurable dimensionless variable, Peclet number. Finally, the comparisons for several reactant species distributions are made. Since a series of parametric studies is carried out, a reference case, served as the comparison base, is chosen in advance. The reference case is selects as 90–90 flow pattern since it is the most common channel configuration studied in the majority of researches. The cell temperature is 353K because it can achieve the best performance. The other detailed operating conditions for this reference case are listed in Table 4.1.
Table 4.1 Operating conditions of the reference case.
Parameter Values Reversible voltage 1.1154 (V)
Operating Temperature 353 (K) Atmosphere Pressure 101.3 (kPa)
Outlet Pressure 101.3 (kPa)
Flow direction Counter current Anode inlet conditions
Gas H2
Relative humidity 100%
Mass fraction of H2/H2O 8.868/91.132
Temperature 353 (K)
Flow rate 796.08 (ml/min) Cathode inlet conditions
Cathode Gas O2
Relative humidity 100%
Mass fraction of O2/H2O 60.89/39.11
Temperature 353 (K)
Flow rate 658.992 (ml/min)
4.1 Fundamental Behaviors of Reference Case
The flow field, mass transfer and distributions of temperature and current density characteristics inside a PEMFC are analyzed in this section in order to obtain a basic understanding of how three-dimensional flow and transport phenomena influence the electrochemical process and investigate the interdependence between variables.
Figure 4.1 shows the polarization curve of the reference case. With cell temperature of 353K, the reversible voltage is set to be 1.1154V by the equation obtained from reference [33]. The corresponding output current densities are collected for every 0.1V for the operating voltage ranged from 0.4 to 1.0V. In the high operating voltage regime, the current density drops substantially due to activation losses. Activation overpotential caused by these losses represents the voltage sacrificed to overcome the activation barrier associated with the electrochemical reaction. In the middle of such curve where performance is dominated by ohmic losses, the current density increases linearly with the decrease of voltage because the ohmic overpotential is linearly proportional to the current production. Reactant gases depletion and products accumulation would adversely affect the performance at high current densities. Such voltage
drop is referred to concentration losses. However, performance in this case does not drop considerably around the end of this curve, indicating that liquid water does not contribute to the voltage drop that will be discussed later.
The mass diffusion, which is the main gas transport mechanism within porous media, is crucial to PEMFC performance since reactants must pass through GDL to reach catalyst layer for electrochemical reaction. However, in flow channels, insufficient mass convection may cause fuel depletion in downstream areas. Therefore, the performance improvement correlates closely with the mass diffusion and convection rate within the channels. In fluid dynamics, the Peclet number, a dimensionless group relating the rate of mass convection of a flow to its rate of mass diffusion, is defined as Pe = VL/D, where V is the flow velocity, L the characteristic length, and D the mass diffusivity. Therefore, flows with smaller Peclet number are having higher mass diffusion rate that is favorable for gas transport in porous media. Figure 4.2 shows the specified number for each bend entrance along anodic flow channel where the corresponding values of Peclet number are determined. The inverse of local Peclet number (Pe-1) in streamwise direction at each anodic bend entrance is shown in Fig. 4.3. It shows that the inverse of Peclet number (Pe-1) increases along the anodic channel from the inlet toward the outlet because of the decrease of gas velocity and the increase of gas temperature caused by the electrochemical reaction. Pe-1 at the first bend is 0.0113 and then increases to 0.0151 at the last bend, implying a significant enhancement of mass diffusion.
With the increase of mass diffusion rate, local uniformity of variables would be improved in the downstream channel. For further analyses, the distributions of current density, velocity and other variables at operating voltage of 0.4 will be discussed below.
Figure 4.4 shows the distribution of current density in the membrane. It shows that the current density decreases slightly from anodic inlet toward the outlet, implying that the current density distribution is principally dictated by the hydrogen concentration rather than the oxygen’s. This is because under fuel-rich supply, hydrogen becomes dominant due to considerably higher reference current density compared to oxygen. For reference case, the maximum current density is 1050 mA/cm2 located near the anodic inlet and minimum is 502 mA/cm2 occurred in the vicinity of anodic exit, and that presents a difference of 548 mA/cm2 in current density on global views.
However, the local current variation occurred at the marginal rib with a value of 185 mA/cm2 is much larger than that at channels with a value of 36 mA/cm2. It implies that there exists a higher non-uniformity current distribution at the rib.
Also, the non-uniformity in local distribution appears to be significant near the first bend and becomes less non-uniform as the current density decreases along the channel. Moreover, the marginal region of the membrane, where large variation of current density is observed, non-uniformity of temperature and water content would exist correspondingly and that will be discussed later.
Figure 4.5 shows the distributions of velocity in both anodic and cathodic channels. In the anodic channel, it can be noticed generally that the velocity magnitude decreases slightly along their flow directions from inlet toward outlet due to the water back-diffusion. A close look at the figure reveals that the channels give rise to the additional minor losses due to the presence of secondary flows at the bending areas. The losses would cause energy dissipation within the cell. However, it is desirable for the mass transfer point of view for enhancement.
The change in the current density can cause local variation in the cell
temperature, and the distribution uniformity of temperature is important for minimizing the thermal stress on the membrane so that its lifetime would be extended. The temperature distributions at anodic and cathodic channels and in the membrane are shown in Fig. 4.6. The uniformity appears to be good at anodic and cathodic channels with the total variations of 1.5K and 2.2K, respectively, indicating that the local thermal stress is not significant. From the comparison of temperature variation at both channels, it indicates that the variation of temperature on the cathode is slightly higher than on the anode since cathode is the place where most heat is produced. As shown in Fig. 4.6(b), the overall temperature distribution in the membrane basically follows the current density distribution shown in Fig. 4.4. This is because heat production and current generation are positively correlated. However, a small discrepancy still exists at the marginal region, where temperature is relatively low but current density is considerably high. This is because the isothermal conditions are imposed on the energy equation rather than zero-flux conditions. This figure also shows that the temperature reaches 356K at the second, third and fourth channels, then, decrease slowly to the exit region. The good uniform distributions of reactant gases are important to the improvement of PEMFC performance, since they can effectively alleviate the flooding effect, leading to the high consumptions of reactant gases.
Figure 4.7 shows the distribution of hydrogen concentration on the interface between anodic GDL and catalyst layer. It shows that the hydrogen concentration increases along the flow direction, indicating that its concentration increases with decreasing current density. The concentration is maintained at 0.4 in the area above the first one third of the anodic channel, and slowly increases to 0.51 near the anodic exit region. Similar phenomenon also was
predicted by [38].
The distribution of oxygen concentration on the interface between cathodic GDL and catalyst layer is shown in Fig. 4.8. Unlike the hydrogen one, the oxygen concentration decreases along its flow direction with a highest value of 0.92 near the inlet region and the lowest one of 0.2 at the outlet region. From these results, it indicates that the property variations in the area above the cathodic channels are greater than that of anodic ones. From Figs. 4.4, 4.7 and 4.8, it indicates that the distribution of current density is inversely proportional to both the hydrogen and oxygen local concentrations.
Water generated in cathode catalyst layer may back-diffuse from cathode to anode, if the cathode is well-hydrated. However, protons moving from the anode to cathode would drag water with them, and this force is called electro-osmotic drag force. This phenomenon implies that, especially at high current densities, the electrolyte at the anode would become dried out–even if the cathode is fully hydrated. Therefore, a uniform distribution of water in both vapor and liquid states is necessary to maintain both the anode and cathode in well-humidified conditions. Figure 4.9(a) presents the distribution of water vapor on the interface between cathodic GDL and catalyst layer. It shows that the water vapor concentration increases positively with current density since both are the products of reaction, and the highest value of water vapor concentration, 0.8, is observed near the anodic inlet where local current density is high. Also, the lowest value of water vapor concentration, 0.1, is located in the last one third of the anodic channel. The flux of water vapor in the membrane from anode to cathode is shown in Fig. 4.9(b). Along the anodic channel, flux of water vapor from anode to cathode increases, implying that back-diffusion becomes insignificant in anodic downstream areas due to low
water vapor concentration. Therefore, it is noticed that in the anode, an amount of water is consumed by penetrating to the cathode due to osmotic drag, resulting a decrease of water vapor concentration. Correspondingly, hydrogen concentration increases from anodic inlet toward the outlet for species conservation.
The distribution of liquid water on the interface of cathodic GDL and catalyst layer is illustrated in Fig. 4.10. As mentioned earlier, saturation of liquid water is defined as a ratio of liquid volume to the pores’ one. This figure shows that liquid water only occurs at the first rib on the left-hand side where the local current density is relatively high comparing to other areas.
According to the results from reference [39], liquid water formulation is inversely proportional to the cell temperature, and, once the cell temperature is increased beyond 348K, most of the water vapor is removed out of the channels with the gas flow without liquid water formation and droplet accumulation.
Therefore, only small amount of liquid water exists in the PEMFC as shown in Fig. 4.10. From the comparison of liquid water profile with the distributions of other variables, it shows no similarity between them, indicating that flooding effect is not significant within the cell.
Water content is given as the ratio of the number of water molecules to the one of charge (SO−3H+) sites. This ratio indicates how well the membrane is
Water content is given as the ratio of the number of water molecules to the one of charge (SO−3H+) sites. This ratio indicates how well the membrane is