In a fuel cell, a wide variety of different physical phenomena occur which are interdependent. In the following, these physical phenomena will be presented with the respective governing equations which are used in the single channel model realized in ANSYS Fluent.
Electrochemistry Modeling
The electrochemistry model adopted in ANSYS Fluent is based on work of Kulikovsky et al., Mazumder and Cole, and Um et al. 9–11. The main task of the model is to compute the rate of reactions at the anode and cathode. This is done via the surface overpotential, which is the difference between the phase potential of the solid and the phase potential of the membrane. The resulting electron transport through the solid conductive materials and the protonic transport through the membrane is modeled as follows.
$$\nabla \cdot \left({\kappa }_{sol}\nabla {V}_{sol}\right)+{j}_{sol}=0$$
1
$$\nabla \cdot \left({\kappa }_{M}\nabla {V}_{M}\right)+{j}_{M}=0$$
2
Here \(\kappa\) is the electrical conductivity, \(V\) the electric potential and \(j\) the volumetric transfer current. The subscripts \(M\) and \(sol\) denote the values at the membrane and the solid respectively. for the membrane and solid material respectively. \(\nabla\) is the Nabla operator. For the solid phase, \({j}_{sol}=+ {j}_{cat}\) on the cathode side and \({j}_{sol}=- {j}_{an}\) for the anode side. For the membrane phase, \({j}_{M}= -{ j}_{cat}\) on the cathode side and \({j}_{M}=+ {j}_{an}\) for the anode side. These terms are so-called source terms and are modeled with the Butler-Volmer equation as follows.
$${j}_{an}=\left({{\zeta }}_{\text{a}\text{n}} {i}_{an}^{ref}\right) {\left(\frac{\left[A\right]}{{\left[A\right]}_{ref}}\right)}^{{\gamma }_{an}}{\cdot e}^{\frac{{a}_{an}F{\eta }_{an}}{\overline{R}T}}$$
3
$${j}_{cat}=\left({{\zeta }}_{cat} {i}_{cat}^{ref}\right) {\left(\frac{\left[C\right]}{{\left[C\right]}_{ref}}\right)}^{{\gamma }_{cat}}{\cdot e}^{\frac{{a}_{cat}F{\eta }_{cat}}{\overline{R}T}}$$
4
Here \({i}^{ref}\) is the reference exchange current density, \(\zeta\) the specific active surface area, \(\gamma\) the concentration dependence, \(a\) the charge transfer coefficient, \(\stackrel{-}{R}\) the universal gas constant, \(F\) the Faraday constant and \(T\) the temperature. The subscripts \(cat\) and \(an\) denote the values at the cathode and anode respectively. \(\left[A\right]\) and \(\left[C\right]\) represent the molar concentration of the species upon which the reaction rates depend for the anode and cathode respectively. The local surface overpotential \(\eta\) quantifies the activation loss. On the anode side, it can be modeled over the difference between the solid and membrane potentials. At the cathode side, the open-circuit voltage \({V}_{oc}\) is further subtracted from the difference to account for the gain in electrical potential.
$${\eta }_{an}={V}_{sol}-{V}_{M}$$
5
$${\eta }_{cat}={V}_{sol}-{V}_{M}-{V}_{oc}$$
6
Current and Mass Conservation
The chemical reaction on the anode and cathode side occurs at the so-called triple-phase-boundary of catalyst, membrane and species. In the triple-phase boundary at the catalysts on the anode and cathode sides, the creation of the different species is modeled as follows.
\({S}_{H2}=-\frac{{M}_{H2}}{2 F} {j}_{an}\) | (7) |
\({S}_{O2}=-\frac{{M}_{O2}}{4 F} {j}_{cat}\) | (8) |
\({S}_{H2O}=\frac{{M}_{H2O}}{2 F} {j}_{cat}\) | (9) |
Here \(S\) is the volumetric source term and \(M\) the molar mass of hydrogen (\({H}_{2}\)), oxygen (\({O}_{2}\)) and water (\({H}_{2}O\)) respectively. \({S}_{H2}\) and \({S}_{O2}\) are negative as they are being consumed at the triple-phase boundary while \({S}_{H2O}\) is positive as it is being created. As the total electrical current which is generated at the anode and cathode is equal, the current conversation can be stated as follows.
$${\int }_{an}^{ }{j}_{an}dV={\int }_{cat}^{ }{j}_{cat}dV$$
10
Heat Source
The heat generated in the fuel cell stems from the enthalpy change due to the electrochemical reactions \({h}_{react}\), and the voltage losses which are transformed to heat. The voltage losses can be differentiated into activation losses and ohmic losses. Furthermore, the enthalpy change caused by the condensation and vaporization of water \({h}_{L}\) must be considered. Thus, the heat source term \({S}_{h}\) can be expressed as follows.
$${S}_{h}={h}_{react}-{j}_{an,cat} {\eta }_{an,cat}+{I}^{2}{Z}_{ohm}+{h}_{L}$$
11
Here \(I\) is the electrical current and \({Z}_{ohm}\) the ohmic resistivity.
Liquid Water Formation and Transport
The liquid formation and transport model is based on work from 12,13. In this approach, the conservation equation for the volumetric water saturation \(s\) governs the water formation and transport.
$$\frac{\partial \left({\Phi }{\rho }_{l}s\right)}{\partial t}+\nabla \cdot \left({\rho }_{l}{\overrightarrow{u}}_{l}s\right)={r}_{H2O}$$
12
Here \({r}_{H2O}\) is the condensation rate for water, \(u\) the velocity, \({\Phi }\) the porosity and \(\rho\) the density. The subscript \(l\) denotes values for liquid water. The condensation rate is expressed as follows.
$${r}_{H2O}={c}_{r}\text{max}\left(\left[\left(1-s\right)\frac{{p}_{H2O}-{p}_{s}}{\stackrel{-}{R}T} {M}_{H2O}\right], \left[-s {\rho }_{l}\right] \right)$$
13
Here \({c}_{r}\) is the condensation rate constant, \({p}_{H2O}\) the water vapor pressure and \({p}_{s}\) the saturation pressure. For the porous zones, the convective term in Eq. (12) is replaced by a capillary diffusion term.
$$\frac{\partial \left({\Phi }{\rho }_{l}s\right)}{\partial t} +\nabla \cdot \left[\left({\rho }_{l}\frac{K {s}^{3}}{{\mu }_{l}}\frac{d{p}_{cap}}{ds} \nabla s\right)\right]= {r}_{H2O}$$
14
Here \(t\) is the time, \(K\) the absolute permeability, \(\mu\) is the viscosity and \({p}_{cap}\) the capillary pressure. The capillary pressure is computed through the Leverett function and thus is dependent on the water saturation and the wetting phase 14.
$${p}_{cap}= \left\{\begin{array}{c}\frac{{\Gamma }cos\theta }{{\left(\frac{K}{{\Phi }}\right)}^{0.5}} \left(1.417 \left(1-s\right)-2.21{\left(1-s\right)}^{2}+1.263{\left(1-s\right)}^{3}\right)\\ \frac{{\Gamma }cos\theta }{{\left(\frac{K}{{\Phi }}\right)}^{0.5}} \left(1.417 \left(1-s\right)-2.21{s}^{2}+1.263{s}^{3}\right)\end{array}\right.\theta <90^\circ \theta >90^\circ$$
15
Here \(\theta\) is the contact angle and \({\Gamma }\) the surface tension.
The model was validated through experimental data from Iranzo et al. 15. As material and geometric parameters of the ZBT fuel cell stack used in the experiments could not be obtained, the given fuel cell parameters from Iranzo et al. were used for the initial validation of the model. Consequently, the operational parameters were also chosen in compliance with the experimental data from Iranzo et al. An overview of the model parameters for the bipolar plates, gas diffusion layers (GDL), catalytic layers and membrane can be found in Table 1. A mesh sensitivity study was conducted and a mesh of 375,000 elements was chosen to minimize the error while still maintaining an adequate model solving time.
Table 1
Model parameters used in the ANSYS Fluent model taken from 15.
Parameter | Value | Units |
Bipolar plate density | \(1990\) | \(kg {m}^{-3}\) |
Bipolar plate specific heat capacity | \(710\) | \(J k{g}^{-1}{K}^{-1}\) |
Bipolar plate thermal conductivity | \(120\) | \(W{m}^{-1}{K}^{-1}\) |
Bipolar plate electric conductivity | \(92600\) | \({{\Omega }}^{-1}{m}^{-1}\) |
Bipolar plate thickness | \(9.5\) | \(mm\) |
Bipolar plate – GDL contact resistance (serpentine) | \(3.52\times {10}^{-7}\) | \({\Omega }{m}^{2}\) |
GDL density | \(321.5\) | \(kg {m}^{-3}\) |
GDL porosity | \(0.82\) | \(-\) |
GDL electric conductivity | \(280\) | \({{\Omega }}^{-1}{m}^{-1}\) |
GDL viscous resistance (anode) | \(1\times {10}^{12}\) | \({m}^{-2}\) |
GDL viscous resistance (cathode) | \(3.86\times {10}^{12}\) | \({m}^{-2}\) |
GDL wall contact angle | \(110\) | \(deg\) |
GDL thickness | \(420\) | \(\mu m\) |
Catalytic layer surface-to-volume ratio | \(1.25\times {10}^{7}\) | \({m}^{2} {m}^{-3}\) |
Catalytic layer thickness (anode) | \(6.0\) | \(\mu m\) |
Catalytic layer thickness (cathode) | \(12.0\) | \(\mu m\) |
Membrane density | \(1980\) | \(kg {m}^{-3}\) |
Membrane thermal conductivity | \(0.16\) | \(W {m}^{-1}{K}^{-1}\) |
Membrane equivalent weight | \(1100\) | \(kg kmo{l}^{-1}\) |
Membrane thickness | \(175\) | \(\mu m\) |
Concentration exponent (anode) | \(0.5\) | \(-\) |
Concentration exponent (cathode) | \(1.0\) | \(-\) |
Reference exchange current density (anode) | \(448\times {10}^{5}\) | \(\mu A c{m}_{Pt}^{-2}\) |
Reference exchange current density (cathode) | \(448\) | \(\mu A c{m}_{Pt}^{-2}\) |
The results of the ANSYS fuel cell model compared to the experimental data of Iranzo et al. can be seen in Fig. 2. The percentage wise error reaches a maximum at the end of the analyzed interval of 6.5%. This is satisfactory, especially since the fuel cell system may primarily be operated at lower current densities ranges which are better suited for aircraft applications 2. Additionally, the polarization curve of the ZBT stack is shown in Fig. 2 which is used in the experimental part of this study.