An Open-Source Toolbox for Multiphase Flow Simulation in a PEM Fuel Cell

A proton exchange membrane (PEM) fuel cell is an electrolytic cell that can convert chemical energy of hydrogen reacting with oxygen into electrical energy with no greenhouse gases generated in the process. To satisfy increasingly demanding application needs, providing fuel cells with better performance and higher efficiency are of paramount importance. Computational fluid dynamics (CFD) analysis is an ideal method for fuel cell design optimization. In this paper, a comprehensive CFD-based numerical tool that can accurately simulate multiphase flow and the major transport phenomena occurring in a PEM fuel cell is presented. The tool is developed using the Open Source Field Operation and Manipulation (OpenFOAM) software (a free open-source CFD code). This makes it flexible and suitable for use by fuel cell manufacturers and researchers to get an in-depth understanding of the cell working processes to optimize the design. The distributions of velocity, pressure, chemical species, Nernst potential, current density, and temperature at case study conditions are as expected. The polarization curve follows the same trend as those of the I-V curves from numerical model and experimental data taken from the literature. Furthermore, a parametric study showed the key role played by the concentration constant in shaping the cell polarization curve. The developed toolbox is well-suited for research and development which is not always the case with commercial code. The work therefore contributes to achieving the objectives outlined in the International Energy Agency (IEA) Advanced Fuel Cell Annex 37 which promotes open-source code modelling of fuel cells. The source code can be accessed at http://dx.doi.org/10.17632/c743sh73j8.1.


Introduction
Owing to their higher power densities, lower operating temperatures, and zero emission, proton exchange membrane (PEM) fuel cells have become an integral part of the energy mix schemes of many countries across the globe. In fact, they are poised to replace some of the conventional power generation devices which depend on fossil fuels.
In this work, a 3-D, non-isothermal and multiphase flow OpenFOAM model of a complete cell (including all fuel cell components) is developed and validated. The work aims at filling the research knowledge gap in comprehensive open-source modelling approaches for PEM fuel cells. The model developed is an extension of an existing 3-D non-isothermal single-phase flow model introduced in a previous work (Kone, Zhang, Yan, Hu, & Ahmadi, 2017a) to include liquid water formation, transport, and their effects in PEM fuel cells. The moisture diffusion approach of multiphase flow modelling in a PEM fuel cell based on the unsaturated theory as described elsewhere (Kone et al., 2017b) has been adopted. It is believed that this modelling approach is more appropriate since it only requires one additional governing equation to the existing single-phase flow model. As for the mixture quantities, these are calculated according to the multiphase mixture model (Kone et al., 2017b). Like the previous model (Kone et al., 2017a), the present model has been developed using OpenFOAM version 4.0.

Multi-Regions and Multi-Physics
The 3-D geometry and mesh of the fuel cell are shown in Figure 1 and the detailed dimensions of the components are given in Table 1. The model consists of two bipolar plates and a membrane electrode assembly (MEA). The bipolar plates accommodate the gas flow channels. The MEA comprises of a membrane with a negatively charged electrode (anode) attached on one side and a positively charged electrode (cathode) attached on the other side. The electrodes are each made up of a gas diffusion layer and a catalyst layer. For converting chemical energy into electrical energy, hydrogen is fed to the cell on the anode side while oxygen is supplied on the cathode side. It is the electrochemical reactions between these two gases and the catalyst particles which provide the electrical energy.  Vol. 11, No. 3;2018  Since the physical quantities transported within the fluid and solid regions of a PEM fuel cell are different, the conservation equations that need to be solved in these regions can therefore vary. As shown in Figure 2, the fluid region encompasses the gas flow channels, the gas diffusion layers, and the catalyst layers, whereas the membrane and the bipolar plates make up the solid region. The reactants are transported in the gas flow channels to the gas diffusion layers and then to the surface of the catalyst layers. The product water moves from the surface of the catalyst layers into the pores of the gas diffusion layers where it is then transported into the gas flow channels. The porous structure of the gas diffusion and catalyst layers facilitates the removal of the heat generated by the electrochemical reactions towards

Assumptions
For the present model, the fuel cell components are treated as homogeneous and isotropic materials. No reactant gases can permeate the membrane. The electrochemical reactions are assumed to occur at the membrane-catalyst layer interface (Beale et al., 2016). Activation and concentration overpotentials are neglected in the anode side (Barbir, 2005). In the bipolar plates, joule heating is neglected because of high heat conductivity and the electrical potential distribution is assumed to be constant because of high electrical conductivity of the bipolar plate material (Ju, Meng, & Wang, 2005;Sivertsen & Djilali, 2005). The gas flow is treated as laminar and incompressible due to low velocities at steady-state operating conditions. The individual gases and the gas mixture are treated as ideal gases. No heat enters or leaves through the outer walls of the entire cell (Beale et al., 2016). The two-phase mixture is considered as a fine mist, and thus, the velocity of the liquid is equal to the velocity of the gas inside the gas flow channels where the two phases share the same pressure and temperature fields. Water vapour is assumed to be in equilibrium with liquid water at the interface separating the two. The liquid phase consists of liquid water cis.ccsenet.org Computer and Information Science Vol. 11, No. 3;2018 13 only.

Governing Equations
The conservation of mass is expressed as: The conservation of momentum using Navier-Stokes equations reads: where the momentum source term S M , is equal to Darcy resistance in the porous media given by: The conservation of chemical species is expressed as: where the sum of the mass fractions of the other species is used to calculate the mass fraction of the inert species in each electrode, = 1 − ∑ The conservation of energy reads: where and are the energy source terms due to the heat released by the electrochemical reactions and water phase change, respectively.
Liquid water transport is written as (Lei Xing et al., 2014): where D l and S l are the diffusivity of liquid water and the mass source term due to phase change, respectively The cell voltage is expressed as:

Constitutive Equations
The fluid densities are obtained using the ideal gas assumption as follows: The molar fraction is related to the mass fraction by: The specific heat capacity of the gas mixture is defined by: where is the molar heat capacity determined by (Todd & Young, 2002): The dynamic viscosity of the gas mixture is obtained using Reichenberg's expression (Poling et al., 2000): The diffusion coefficient of a gas species i with respect to the total gas mixture in an ordinary diffusion is calculated by (Fairbanks & Wilke, 1950): where the binary diffusion coefficient D ji for each component of the gas mixture is obtained by (Fuller, Schettle.Pd, & Giddings, 1966): For simultaneous ordinary and Knudsen diffusion in the porous structure of the gas diffusion layer and the catalyst layer, the effective global-diffusion coefficient is defined as (Welty, Wicks, Wilson, & Rorrer, c2008): is the effective diffusion coefficient of the individual species of the gas mixture in the porous medium and it is expressed as follows (Welty et al., c2008): Likewise, the effective Knudsen diffusion coefficient of the individual species is expressed as: where the Knudsen diffusion coefficient of a species i, D i-Knud , is given by (Geankoplis, 1993): The local open-circuit or ideal potential given by the Nernst equation, reads: where the standard electrode potential E 0 is defined as: The cell activation overpotential is expressed as (Barbir, 2005): where the cathode exchange current density 0 is given by (Y. Wang, Chen, & Cho, 2013): The cell ohmic overpotential is expressed as (Barbir, 2005): R Ω is the area specific resistance of the cell defined by: where the membrane ionic conductivity σ i which is a function of the membrane water content λ is determined by (Springer, Zawodzinski, & Gottesfeld, 1991): the membrane water content is defined as (Springer et al., 1991): where the activity of water vapour in the gas mixture a is calculated as (Springer et al., 1991): the saturation pressure p sat is given by (Springer et al., 1991): The species electrochemical mass flux is defined as: The capillary pressure of liquid water in the porous media is written as (Y. Wang et al., 2013): where J(s) denotes Leverett function defined by (Y. Wang et al., 2013): The effective diffusivity of the gas species is correlated with porosity and saturation (Nam & Kaviany, 2003): where f(ε) and g(s) are the normalized porosity and saturation functions, calculated by (Dawes, Hanspal, Family, & Turan, 2009): The density of the two-phase mixture is defined by: cis.ccsenet.org Computer and Information Science Vol. 11, No. 3;2018 The specific heat capacity of the two-phase mixture is defined by: The thermal conductivity of the two-phase mixture is defined by: The specific latent heat of evaporation or condensation of water is determined by (Nguyen & White, 1993): 9T + 3.44 × 10 −3 T 2 + 2.54 × 10 −6 T 3 − 8.98 × 10 −10 T 4 (46)

Boundary Conditions
Boundary conditions are applied at all the boundaries of the computational domain, as well as at the internal interfaces between the different components, as shown in Figure 3. Table 2 contains specified boundary and initial values used.
Dirichlet boundary conditions are applied at the fluid inlets meaning that the inlet values for velocity, temperature, and mass fractions are prescribed at the inlet of both the anode and the cathode gas flow channels. The inlet velocity is a function of the species stoichiometric flow ratio, the current density, the electrode active area, and the gas flow channel cross-section area. It is calculated according to (Liu, 2012): The outlet value is prescribed for pressure at the outlet of both the anode and the cathode gas flow channels. Neumann boundary conditions are applied for all other variables, with their gradients assumed to be zero in the flow direction.
Insulation conditions are set for the flow of mass and chemical species at the interfaces between the fluids and the membrane as the latter is impermeable to the reactant gases.
Impermeability, no-slip and no-flux boundary conditions are applied at all solid-fluid interfaces within the computational domain. A zero-gradient temperature is applied at all the external surfaces of the cell.

Overview of the Toolbox
The work presented in this paper extends the 3-D, non-isothermal, and single-phase flow model for a PEM fuel cell previously introduced (Kone et al., 2017a), to include liquid water formation, transport, and their effects in the fuel cell. The toolbox developed follows the same structure and organisation as OpenFOAM, as shown in Figure  5. One limitation of the OpenFOAM platform is that it has no multiphase flow PEM fuel cell module. The present tool in this study actively include the multi-phase flow programme. The software has two major components: pemfcMultiphaseNonIsothermalSolver and run (see Figure 5). The pemfcMultiphaseNonIsothermalSolver directory contains the source code for the solvers and libraries, whereas the run directory stores the constructed simulation case files.
cis.ccsenet.org Computer and Information Science Vol. 11, No. 3;2018 Figure 5. File structure of the toolbox A pemfcMultiphaseNonIsothermalSolver.c file that includes the main loop of the program is stored in the applications sub-directory. The files for the necessary libraries to link to during runtime are stored and compiled in the lib sub-directory. These include classes for geometry and mesh manipulation, multicomponent gas diffusivities, etc. The src sub-directory stores various program files which contains specific instructions for the initialization, as well as algorithms for solving the different field variables. The 0 sub-directory is a time subdirectory that holds the initial field data files. The config sub-directory stores the data files for geometry construction and mesh creation. The constant sub-directory contains information files on mesh and boundary conditions, along with physical properties. The system sub-directory stores the control files for the solvers, schemes, and solutions.

Grid Refinement Study
A mesh independence study was carried out to investigate the dependency of the results of the model on the mesh. Thus, the grid that was used for the case study was refined twice by adding 20% of the cells and 40% of the cells in every direction. Computations at case study conditions were then performed on these refined meshes and the temperature field solutions were compared with each other. The computation results are shown in Table 6. The maximum difference of temperature solution is 0.08 % for the compared meshes, indicating that the case study mesh (Mesh 1) provides sufficient confidence.

Comparison with Literature Model Results and Experimental Data
Figure 6 compares the cell polarization or I-V curve obtained using the concentration constant values of 0.2, 0.25 and 0.3 at case study conditions, to those of the numerical model and experimental data reported by (Yuan et al., 2010). This is because some of the parameters and properties used in the present model were directly borrowed from the model in (Yuan et al., 2010). In addition, other parameters such as the transfer coefficient, the reference exchange current density, etc. were adjusted to achieve a good agreement between the compared models.
From Figure 6, it is shown that the compared polarization curves follow similar trends, confirming the validity of the present model. The difference between the I-V curves at high current densities can be attributed to, amongst other factors, the combined effects of the cell inlet flow rates, the cell design details and material properties, and the concentration losses. In fact, the effective value of the concentration constant seen in real fuel cell operation is often much larger than the theoretical value. Therefore, in many modelling work, this constant is obtained empirically (O'Hayre et al., 2006).
Based on the present comparison, it is clear that the choice of the value for the concentration constant is decisive in obtaining a good agreement between the results of the present work and the numerical and experimental results at high densities from a previous study (Yuan et al., 2010). A c value of 0.3 provides the best agreement between the results. Figure 6. Cell I-V curves compared to the numerical and experimental results from a previous study (Yuan et al., 2010) At low current densities, increasing the value of c has negligible effects on the cell potential, whereas at the current densities higher than 4000 A/m 2 , the effects get significant and results in decreased cell potential. Figure 7 compares the activation, ohmic, and concentration overpotentials. The values for the concentration constant used above were carefully selected to ensure realistic proportions among these three types of overpotentials in the fuel cell. The activation overpotential constitutes the largest potential loss at any current density because of slow reaction rates at the electrodes. The potential loss associated with the ohmic resistance become significant at moderate current densities due to resistance to the flow of electrons and protons. The potential loss caused by mass transport rises at high current densities as the electrodes are rapidly depleted of the reactants by the electrochemical reactions.  Vol. 11, No. 3;2018 23

Discussion of the Case Study Results
The post-processing and the visualization of the results were achieved using ParaView, which is an open-source code visualization software that is used for postprocessing in OpenFOAM through the paraFoam utility, supplied with OpenFOAM. Figure 8 (a) and (b) display the velocity profiles along the gas flow channels for the fuel and air, respectively. The highest velocity is observed at the centre lines of the channels at the cell outlet while the lowest velocity is seen at the walls. These contours agree with those of fully developed laminar flow.  Water in a PEM fuel cell originates from water vapour in humid reactant gases and the oxygen reduction reaction (in the present work, it is assumed that the product water is in vapour form). Water vapour condensates into liquid water when its partial pressure is higher than the saturation pressure which is greatly affected by the local temperature. Therefore, liquid water in the fuel cell results from water vapour condensation. From Figure 12 and 13, it is observed that the saturation of liquid water increases from the cell inlet to the cell outlet. The reason why saturation swells along the flow direction is twofold: increased water vapour partial pressure along the flow direction because of the build-up of water vapour, and the nature of liquid water flow. show the distribution of hydrogen and oxygen mass fractions, respectively. As the reactants are consumed during the electrochemical reactions, their mass fractions decrease from the cell inlet to the cell outlet, and from the gas flow channels to the catalyst layer. The concentration of reactant gases in the electrodes is directly proportional to pressure (see Figure 15 and 16). It is also noticed that, compared to the results of the singlephase flow model reported in (Kone et al., 2017a), the mass fractions of the reactant gases have further decreased in the multiphase flow case because the transport of the gas species is hindered by the presence of liquid water (see Figure 17 and 18). In addition, the clogging of the diffusion passages by water droplets results in reduced diffusion rates.  Figure 19 (a) illustrates the Nernst potential decreasing from the cell inlet to the cell outlet as the anode is depleted of hydrogen because of the combined effects of hydrogen consumption and water vapour production by the electrochemical reactions. Figure 19 (b) shows the distribution of the observed local current density. It decreases in the flow direction because of the reduction in the Nernst potential in the same direction. A significant reduction is noticed at the corner edges of the outlet since these areas are devoid of active reaction surface. Compared to results of the single-phase flow model reported in (Kone et al., 2017a), a further reduction in both the Nernst potential and the local current density is noticed (see Figure 20 and 21, respectively). This can be attributed to a reduced effective reaction surface area caused by the presence of water droplets. Figure 19 (c) displays the local temperature distribution. This is affected by ohmic heating, and the heat generated during the electrochemical reactions and water phase change. The additional heat released by water phase change results in a significant increase in the average local temperature, in comparison with the average local temperature obtained with the single-phase flow model (Kone et al., 2017a) (see Figure 22).  (Kone et al., 2017a). The difference between the two lines is not significant at low current densities. At moderate and high current densities, however, the gap between the two lines starts to widen. For an cis.ccsenet.org Computer and Information Science Vol. 11, No. 3;2018 identical current density, the cell potential is lower with the multiphase flow model than that with the single-phase flow model. This is because the multiphase flow model accounts for liquid water effects on the gaseous reactant transport, which causes higher concentration losses. On the other hand, the single-phase flow model cannot accurately predict the cell potential drop when the local current density draws closer to the limiting current density.
In other words, while the existence of liquid water has beneficial effects in terms of membrane hydration, it also blocks the gas diffusion passages, reduces the diffusion rate and the effective reaction surface area, and hence leads to a deterioration in fuel cell performance. Figure 23. Comparison of the cell polarization curve with the single-phase flow model curve in (Kone et al., 2017a)

Conclusion and Outlook
An open-source code toolbox for the numerical simulation of multiphase flow in a PEM fuel cell has been developed using OpenFOAM. The toolbox includes a main program, relevant library classes, and a constructed simulation case for a co-flow galvanostatic simulation. The mathematical model accounts for liquid water formation, transport, and their effects in the fuel cell.
The case study results are as expected. Detailed distributions of velocity, pressure, liquid water saturation, the species mass fractions, Nernst potential, local current density and temperature fields were presented. The results of a mesh independence study revealed that for the coarsest case study mesh, the solution proved to be grid independent. The cell polarization curve obtained with the case study conditions was plotted and compared to a single-phase flow model polarization curve and the polarization curves from a known literature model and experimental data. The comparison with the single-phase flow model confirms that the existence of liquid water in the system reduces the cell performance. The result shows a good agreement when compared with numerically and experimentally obtained data. Furthermore, the outcomes of a parametric study showed the key role played by the concentration constant in shaping the cell polarization curve.
The developed model can, by no means, be considered complete. The present study can still be improved. So far, the membrane water transport (e.g. dissolved water transport) has been simplified. Also, to facilitate numerical implementation, it was assumed that the electrochemical reactions occur at a single interface between the cathode catalyst layer and the membrane. The developed toolbox can therefore be extended to include other features, such as improved membrane and or catalyst models.