Chapter 2 Simulation Model
2.3 Governing Equations
2.3.4 Electrochemical Reaction
At the electrode catalyst layers, current is supplied from electrochemical reactions that reactant species participate at anode and cathode. Hydrogen is oxidized at the anode and oxygen is reduced at the cathode. These two reactions are driven by the potential difference between the solid phase and electrolyte phase, called the activation overpotential. The Bulter-Volmer
equation describes this phenomenon and is expressed as
where αa,j and αc,j are kinetic constants determined from experimentally generated Tafel plots, αk,j are the concentration exponents of the k-th species for the j-th step. [Λk] represents the average interfacial molar concentration of the k-th species, and the subscript ref refers to concentration values at the reference state at which the reference current density is prescribed. The electrode potential, η, is the potential difference between the solid phase (ΦS) and the ionic phase (ΦF) of the electrode.
The term j0,j in Eq. (2-18) is the reference current density with units of A/m3. As a matter of fact, the reference current density at the cathode is several orders of magnitude smaller than the value at the anode. Both the reference current densities at anode and cathode are expressed as [34]
Anode: j0=1×109 (2-19)
Cathode: j0=3×105 exp[0.014189(T-353)] (2-20) 2.3.5 Current Equation
The continuity equation of current in any material leads to:
0 i =
⋅
∇ (2-21)
where i is the current density vector. When the material is a porous electrode, the current density is split into two parts: one flowing though the polymer electrolyte (ionic phase) denoted as iF and the other flowing through the solid parts (electronic phase) of the porous matrix denoted as iS. Fig. 2.3 simply shows the transfer of current within porous media. Eq. (2-18) can be written as
0 i iF +∇⋅ S =
⋅
∇ (2-22)
During the electrochemical reactions, electrons are either transferred from ionic phase to solid phase or vice versa shown in figure 2.3. Thus, Eq. (2-19) can be rewritten as
It is generally assumed that each phase has its own electrical potential. In the ionic phase, the electrical current is denoted as electrolytic term while the current is denoted as electronic term in solid phase. Chemical reaction can only occur at the interface of electrolyte and solid when there’s interaction between these two current components. Under this assumption, application of ohm’s law to Eq. (2-23) yields
where ΦS and ΦF are the electrical potentials of the ionic and solid phases, respectively.
2.3.6 Liquid Water Transport
The foundation of this two-phase model is an additional governing equation for the transport and formation of liquid water, given by
l
where s is the saturation of liquid water, and is defined as the ratio of liquid volume to the volume of pores. ρl is the density of liquid water and ζ is the relative mobility of liquid water. εd is the dry porosity which represents the porosity of the medium in the dry state. In this model where liquid water effect
is considered, pores in porous media are occupied by liquid, and the value of effective porosity decreases as the amount of liquid water in the medium increases. Therefore, εeff has the form
)
The first term in Eq. (2-25) represents the storage of liquid water in a certain domain. And the second term represents transport of liquid water by pressure-driven advection.
Membrane water content, λ, is given as the ratio of the number of water molecules to the number of charge (SO−3H+) sites. It shows how intense the electrochemical reaction is, and dictates the electrical conductivity of membrane.
The value of membrane water content is determined by the activity of water vapor at the catalyst/membrane interface. The activity in the water phase, a, is xwP/Psat. xw is the mole fraction of water vapor and the saturation pressure of water is denoted by Psat in units of atm, which can be computed from the curve-fitted expressions provided by Springer et al. [5]
3
Membrane water content is function of water activity for the Nafion 117 membrane by weighing membranes equilibrated above aqueous solutions of various lithium chloride concentrations, and the experimental relationship of water content versus water activity, a, used in this model is
⎪⎩
Eq. (2-28) is employed to calculate the value of λ based on solution of water vapor outside the membrane domain.
When water profile in the membrane has been determined, the membrane electrical resistance and the potential drop across it can easily be obtained. The electrical conductivity of Nafion 117 was proposed by Springer et al. [5], and is given by
Once the value of water content drops less than one, the membrane conductivity is assumed to be constant.
In a PEMFC, electro-osmotic water drag moves liquid water from anode to cathode. As this water built up at cathode, back diffusion occurs, resulting in the transport of water from cathode back to anode. This phenomenon is strongly dominated by capillary diffusivity. Capillary diffusivity depends on factors like the contact angle and saturation, and can be written in terms of the capillary pressure head as
ds
where pc is the capillary pressure head, κl is the permeability of the liquid and μl is the dynamic viscosity of the liquid. Based on the correlation that Springer et al. [5] proposed, capillary diffusivity is function of water content and temperature, and is expressed as [5]
⎪⎩
T ]
In Eq. (2-32), s is considered as a constant with value of 0.0126.
2.3.7 Concentration Loss
During the electrochemical reaction, performance drops when reactant species concentrations are deficient on the reaction surfaces, especially when operating a PEMFC at low operating voltage due to high formation of water that causes clogging. The total concentration loss is expressed as
R reactant molar concentration in catalyst layers and α is the mass transfer coefficient expressing how variations in electrical potential across reaction interfaces changes the reaction rate. The value of α depends on the reaction and electrode material.
One of concentration loss is due to reactant species depletion in the catalyst layers called Nerst potential changes, and it has the form
R
The second way that concentration contributes to concentration loss is due to reaction kinetics, and it has the form
R
Different bend angles and widths of a channel lead to different minor losses;
therefore, fluid velocity through bend decreases. Figure 2.4 illustrates a schematic of the PEMFC mass transport model, indicating that the reactant
concentration on the reaction surface is dominated by fluid velocity in the flow channel. The correlation between gas velocity in the flow channel and gas density on reaction surfaces is expressed as [3]
C the Sherwood number. According to the Eq. (2-37), the geometry of flow patterns is one of the parameters that dictate the performance of PEMFC.
Hence, high performance can be achieved from flow pattern design.
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
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