Three-dimensional multiphase flow computational fluid dynamics models for proton exchange membrane fuel cell: A theoretical development

Abstract

A review of published three-dimensional, computational fluid dynamics models for proton exchange membrane fuel cells that accounts for multiphase flow is presented. The models can be categorized as models for transport phenomena, geometry or operating condition effects, and thermal effects. The influences of heat and water management on the fuel cell performance have been repeatedly addressed, and these still remain two central issues in proton exchange membrane fuel cell technology. The strengths and weaknesses of the models, the modelling assumptions, and the model validation are discussed. The salient numerical features of the models are examined, and an overview of the most commonly used computational fluid dynamic codes for the numerical modelling of proton exchange membrane fuel cells is given. Comprehensive three-dimensional multiphase flow computational fluid dynamic models accounting for the major transport phenomena inside a complete cell have been developed. However, it has been noted that more research is required to develop models that include among other things, the detailed composition and structure of the catalyst layers, the effects of water droplets movement in the gas flow channels, the consideration of phase change in both the anode and the cathode sides of the fuel cell, and dissolved water transport.

Keywords

Computational fluid dynamics modelling multiphase flow numerical proton exchange membrane fuel cell review

Introduction

In recent years, fuel cells have become an important clean energy technology, and thereby a serious contender to replacing some of the traditional power systems which rely on fossil fuels. Fossil fuels produce a significant amount of pollutants and they are rapidly depleting resources. Thus, it is becoming increasingly difficult to ignore the critical role that fuel cells can play in the different energy mix scenarios of many countries.

Besides their environmental advantages, in general, fuel cells offer many other advantages over conventional energy conversion devices. They directly convert chemical energy of the fuel into useful work, without requiring any thermodynamic cycle. Thus, their practical efficiency in direct electrical energy conversion can reach as high as 60%.¹

Proton exchange membrane (PEM) fuel cells which are the focus of this review article have high power density and operate at relatively low temperatures, which make them well suited for automotive power system, as well as power generation systems for buildings and portable electronics.

However, despite this, PEM fuel cell economics, especially the high capital and working expenses associated with their fabrication and testing, constitute a major obstacle to their rapid development. To tackle this issue, computational fluid dynamic (CFD) modelling and simulation can be used to rapidly gain important insights into PEM fuel cell working processes. These include fluid flow, mass and heat transfers, and chemical reactions inside the fuel cell, which can provide critical information needed for the optimization of PEM fuel cells.

Most reviews on PEM fuel cell CFD modelling provide a general discussion of a wide range of CFD models for PEM fuel cell. This review focuses on published 3D multiphase flow CFD models. The idea is to create a focused review, aiming at critically analysing the quality of the models meeting these criteria by discussing their strengths and weaknesses and identifying outstanding issues instead of conducting a broad historical overview of the published literature.

The review is organized as follows. The next section defines four basic concepts that are essential to understanding the subsequent sections. ‘Rationale for PEM fuel cell CFD modelling’ section gives the rationale for the use of CFD in PEM fuel cell modelling. ‘Literature 3D multiphase flow CFD models’ section reviews the literature’s PEM fuel cell 3D multiphase flow CFD models according to their areas of investigation. ‘Multiphase flow models and phase change’ section examines multiphase flow modelling and the implementation of phase change in the reviewed papers. ‘Numerical procedures’ section outlines a number of salient numerical features of the reviewed models. ‘CFD codes for PEM fuel cell modelling’ section provides an overview of the most commonly used CFD codes for PEM fuel cell modelling. ‘Summary and outlook’ section summarizes the review highlighting future trends in the PEM fuel cell modelling.

Essential concepts

Transport phenomena

The study of transport phenomena is a multidisciplinary subject that draws on science and engineering concepts as diverse as fluid mechanics, mass transfer, heat transfer, and electromagnetism. Each one of these concepts deals with a specific transport phenomenon. The transport of momentum is dealt with in fluid mechanics, the transport of mass of various chemical species is considered in mass transfer, the transport of heat is explained in heat transfer, and the transport of charges is depicted in electromagnetism.²

The physical quantities that are transported within a PEM fuel cell are mass, momentum, chemical species, thermal energy, electrical current, and ionic current. To describe the transport of these physical quantities in the different components, the conservation laws in terms of rate of accumulation and generation/consumption are used.

Thus, the transport phenomena occurring within the cell can be represented by the solution of the conservation equations. The generic governing equations for PEM fuel cell models and applicable component are given in Table 1.

Table 1.

Generic governing equations for PEM fuel cell models.³

Equation	Applicable conservation laws	Component
(1)	Mass $\frac{\partial (ɛ ρ)}{\partial t} + \nabla \cdot (ɛ ρ \vec{V}) = S_{m}$	Gas flow channel, gas diffusion layer, catalyst layer
(2)	Momentum $\frac{\partial (ɛ ρ \vec{V})}{\partial t} + \nabla \cdot (ɛ ρ \vec{V} \vec{V}) = - ɛ \nabla P + ɛ μ \nabla^{2} \vec{V} + S_{M}$	Gas flow channel, gas diffusion layer, catalyst layer
(3)	Species $\frac{\partial (ɛ ρ w_{k})}{\partial t} + \nabla \cdot (ɛ ρ w_{k} \vec{V}) = - \nabla \cdot [D_{k} \nabla (ρ w_{k})] + S_{k}$	Gas flow channel, gas diffusion layer, catalyst layer
(4)	Energy $\frac{\partial (ɛ ρ c_{p} T)}{\partial t} + \nabla \cdot (ɛ ρ c_{p} T \vec{V}) = \nabla \cdot (k \nabla T) + S_{E}$	All components
(5)	Electrical charges $- \nabla \cdot (σ_{e} \nabla Φ_{e}) = S_{e}$	All components except gas flow channel
(6)	Ionic charges $- \nabla \cdot (σ_{i} \nabla Φ_{i}) = S_{i}$	Catalyst layer, membrane

PEM: proton exchange membrane.

Multiphase flow

Multiphase flow refers to fluid flows consisting of two or more phases, with a phase being a solid, a liquid, or a gas that coexist with another solid, liquid, or gas.

The most common types of multiphase flows are the two-phase flows and the three-phase flows. The two-phase flows include gas–liquid flow, gas–solid flow, liquid–liquid flow (e.g. oil–water mixtures in pipelines), and liquid–solid flow. The three-phase flows encountered in engineering applications include gas–liquid–solid flow, gas–liquid–liquid flows (e.g. gas–oil–water flows in oil recovery systems), and solid–liquid–liquid flows (e.g. immiscible liquid–liquid reaction in which a solid phase is formed).

The different phases forming a multiphase flow are usually separated by an identifiable boundary which constitutes the interface between two phases. The interfaces are where the transfer of material, momentum, and energy may occur between the phases. The modelling of multiphase flows is immensely challenging and often requires a sound understanding of the various multi-physics phenomena that are involved.

In PEM fuel cells, multiphase transport originates from the production of liquid water by the oxygen reduction reaction and the phase change processes (evaporation and condensation).

CFD

CFD is a sub-discipline of fluid dynamics that employs a set of highly sophisticated numerical methods and algorithms to solve and analyse problems of fluid dynamics and beyond.⁴ The problem under consideration is usually formulated mathematically using partial differential equations.

Since the CFD method is based on computational analysis, before any of the CFD numerical methods and algorithms can be used to compute practical numerical solutions, these must be translated into a computer program often referred to as CFD codes. CFD codes are usually written in a high-level computer programming language such as C or C++.⁵

In general, a CFD-based analysis is performed in five major steps: (1) problem definition, (2) mathematical model formulation, (3) pre-processing and mesh generation, (4) solving, and (5) post-processing. Some of these steps must be repeated multiple times in order to obtain the desired results.

The problem is usually defined in its simplest possible form that still accurately describes the actual real world system under consideration. The mathematical model is used to formulate the problem mathematically, describing the details of the flow. The pre-processing and mesh generation sets the various fields of interest before the start of the simulation and discretize the flow domain into computational mesh consisting of volumes or cells. In the solving step, the differential model gets replaced by a system of linear algebraic equations which are solved using algorithms. Finally, the post-processing offers the means to visualize the simulation data in order to inspect the details of the flow.

PEM fuel cell

The basic operating principle of fuel cells which consist of reversing water electrolysis to generate electricity from hydrogen and oxygen was discovered by William Grove in 1839.⁶ Although the technology has matured significantly, this same principle still defines the operation of current PEM fuel cells.

Thus, simply put, a PEM fuel cell is an electrolytic cell that runs on hydrogen and oxygen. Its working process can be a multiphase reaction process, generating heat while converting chemical energy into electrical energy without producing any greenhouse gases.

Working principle

There are three key components that form a PEM fuel cell. A solid polymer membrane sandwiched between two electrodes, a negatively charged electrode (anode), and a positively charged electrode (cathode). Both the anode and the cathode contain a catalyst layer (CL), a gas diffusion layer (GDL) with micro-porous layers, and a bipolar plate (BP) with gas flow channels (GFC). The schematic representation of the operation of a single cell is shown in Figure 1.

Figure 1.

An illustration of PEM fuel cell operation. PEM: proton exchange membrane.

The cell is fed with hydrogen at the anode side and oxygen at the cathode side. These two gases subsequently react together to provide chemical energy for conversion to electrical energy. However, due to the inefficiency in the energy conversion, a portion of chemical energy is converted to waste heat instead of electricity.⁷

The overall electrochemical reaction as shown in equation (9) is the sum of the reactions occurring at the electrodes (anode and cathode), given in equations (7) and (8), respectively.

At the anode

H_{2} \to 2 H^{+} + 2 e^{-}

(7)

At the cathode

\frac{1}{2} O_{2} + 2 H^{+} + 2 e^{-} \to H_{2} O

(8)

Overall

H_{2} + \frac{1}{2} O_{2} \to H_{2} O

(9)

Electrical current is created by the movement of electrons, resulting from the splitting of hydrogen into protons and electrons at the surface of the catalyst. The electrons travel from the negatively charged anode side to the positively charged cathode side, following an external circuit that connects the two electrodes through a load, whereas the protons travel through the membrane. At the catalyst of the cathode side, they meet with oxygen that is fed on that side to form water and produce heat.³

The product water may be in vapour or liquid form. Liquid water formation depends on the water vapour saturation pressure and temperature. The presence of excess liquid water can be problematic in PEM fuel cells, as it can create a phenomenon known as flooding. Flooding has a negative effect on mass transfer. Hence, to protect the integrity of the fuel cell, the product water and waste heat generated during the operation must be continuously removed.

Liquid water may be transported out of the cathode CL through the cathode GDL and then the cathode GFCs. Fuel cells may be fitted with cooling units at the BPs to deal with the waste heat.

The operating principle of PEM fuel cells may well be simple and basic, but the phenomena occurring during their operation are highly complex. They involve electrochemical reactions, species and charge transport, mass transfer, heat transfer, and multiphase flows. Sound knowledge of these multi-physics phenomena is therefore critical in PEM fuel cells optimization.

Components

As described in ‘Working principle’ section, a single cell consists of a polymer electrolyte membrane sandwiched between two electrodes. As illustrated in Figure 2, each electrode has a BP with GFCs, a GDL, and a CL.

Figure 2.

Schematic of a single-cell PEM fuel cell with straight channels. PEM: proton exchange membrane.

BP

The BPs (also known as the collector or separator plate) are electrically conductive plates which separate the gases in adjacent cells in the stack while connecting the individual cells electrically. They provide housing for the GFCs and a structural support for the stack.³ The heat generated within the membrane electrode assembly (MEA) is transferred to the cooling system through the BPs. Electrons are transferred from/to the GDLs to/from the BPs.

GFCs

The GFCs in the BPs provide the pathways for the transport of the reactant gases to the GDL and product water out of the electrodes. To ensure an even distribution of the fuel and oxidant within the cell, different flow field designs with different advantages exist. The most commonly seen are straight parallel, serpentine, and interdigitated flow fields.^1,8

GDL

The GDLs are porous media located between the BPs and the CLs. Thus, not only they allow uniform distribution of the reactants across the surface of the CLs, but they also provide an electrical connection (electrons transport) between the CLs and the BPs, as well as creating a structural support for the CLs.¹

CL

The CLs are the thin layers attached to both sides of the membrane. It is on the surface of the CLs that the electrochemical reactions that convert chemical energy in electrical energy take place within a PEM fuel cell.

Membrane

The membrane is a solid polymer electrolyte that separates the anode from the cathode in a PEM fuel cell. In other words, it is an electronic insulator that is also impermeable to the reactant gases but allows protons transport. Thus, the desired characteristics of a PEM fuel cell membrane are good electronic insulation, good separation of hydrogen in the anode side and oxygen in the cathode side, and high proton conductivity.

Rationale for PEM fuel cell CFD modelling

As mentioned in ‘CFD’ section, CFD is foremost a numerical method that is used to solve partial differential equations of fluid dynamics. These partial differential equations are non-linear and analytical solutions can only be obtained in a few simple special cases. In order to obtain solutions to most problem of practical interest, a CFD implementation is therefore necessary.⁵

This is applicable to PEM fuel cells as electrochemical devices, because the motion of the reactant gases and product water is governed by non-linear conservation equations derived from first principles of fluid mechanics making analytical solutions only possible in 1D and highly simplified cases.

In addition to the analytical challenge of solving non-linear equations in PEM fuel cell models, the cost of manufacturing fuel cells and the cost associated with their performance testing often render many experimental techniques uneconomical. Thus, much research work on fuel cells tries to improve the cell performance by increasing its efficiency while decreasing manufacturing and testing costs, through CFD techniques.

Indeed, over the past decade or so, CFD techniques have been employed extensively to model PEM fuel cells. The approach is based on numerical approximation of the solutions of the conservation equations of mass, momentum, energy, current, and species transport on a computational domain that is discretized using the finite-volume, finite-difference, or finite-element schemes.

Through CFD analysis, a deeper understanding of the problem under consideration can be obtained. In fact, the importance of CFD modelling and simulation of fuel cell working processes has been repeatedly addressed throughout the literature. For instance, CFD can help understand water management within a PEM fuel cell. Recalling ‘Working principle’ section, water transport is one of the central issues in PEM fuel cell technology.

CFD has many other important applications in PEM fuel cell system analysis and is consequently of immense technological and economic significance. Nevertheless, due to the complexity of the multi-physics phenomena occurring in a complete fuel cell, some simplifying assumptions are still necessary in order to create a numerically tractable model of a complete fuel cell.

Literature 3D multiphase flow CFD models

The past few decades have seen the number of published modelling and simulation work on PEM fuel cells increased dramatically. This has essentially been made possible by recent advances in computational techniques and the availability of high-performance computing resources.

The majority of the literature’s 1D computational PEM fuel cell models consider the direction orthogonal to the membrane, with varying degrees of complexity and details. The computational 2D models not only treat the transport phenomena that occur inside the fuel cell in the direction orthogonal to the membrane, but they also consider another additional direction, usually either the dimension along the GFC or that perpendicular to it.

Although the pioneering 1D models and subsequent 2D models have significantly contributed to the advancement of PEM fuel cell modelling, these models can only be used for particular phenomena in PEM fuel cells due to the 3D nature of the transport processes occurring in a PEM fuel cell. Consequently, they are not suitable for parametric studies or optimization of fuel cell designs.⁹ Moreover, the transports of mass, heat, and charges in porous structures are major issues for these models.

To accurately predict cell performance, a realistic 3D description of the cell geometry is necessary. Thus, 3D computational models have been developed to better understand the phenomena occurring within PEM fuel cells and then predict their performance though a full 3D description of a complete cell or a fuel cell stack requires a considerable computational effort. Three-dimensional models consider effects in both across and along the gas flow channel, in addition to the direction normal to the membrane.

This review intends to examine the majority of the literature’s PEM fuel cell 3D multiphase flow CFD models though it is likely that it may not include all of the published models. The focus is placed on comprehensive models that are relevant in the context of the present review work. The time frame is set to cover any of such models published through the second quarter of 2016, as of the time the present review has been completed. Thus, only the models that satisfy all criteria of 3D, multiphase flow and CFD-based modelling are included.

A direct comparison of the different models is problematic because they vary in their approach and complexity to modelling multiphase flow and phase change, as well as their validation methods. Also, one model may be well suited for specific operating conditions in a specific region, whereas another one would not be. Therefore, the majority of this review delves into discussing the advantages and potential drawbacks of each model.

The bulk of the models treat phenomena such as water and reactant transport, electrochemical reactions, species and charge transport, mass transfer, heat transfer occurring during the cell operation, with the effects of geometry or operating conditions.

Thus, depending on the characteristics of the phenomena being investigated, and the level of details included, the models can be divided into three distinct groups: models for transport phenomena, models for geometry or operating condition effects, and models for thermal effects.

The pie chart in Figure 3 shows that transport phenomenon is the most studied subject in the reviewed articles followed by geometry or operating condition effects and then thermal effects. This trend can somehow be expected since the 3D architecture of PEM fuel cells makes transport phenomena within them a 3D problem inherently. Moreover, liquid water transport and its effects in the different fuel cell components continue to capture the attention of many PEM fuel cell researchers.

Figure 3.

The reviewed papers according to their areas of investigation.

As was pointed out in previous paragraphs, 3D models are more complete and they can be used for design and optimization purposes. The focal point of the literature models for geometry or operating condition effect has been the design of the flow field; namely, the influence of various flow-field designs on cell performance. The majority of the models for thermal effects are concerned with the temperature distribution inside the cell and the effects of this on the fuel cell performance.

Table 2 summarizes the reviewed research papers. It should be mentioned that none of the models listed in Table 2 include all aspects of PEM fuel cell modelling. Overall, the models vary in complexity, because of the modelling of multiphase flow and the implementation of phase change, as well as other modelling assumptions that are used.

Table 2.

Summary of the reviewed research papers.

No.	Authors and year	Areas of investigation	Basic assumptions	Model validation
1.	Shimpalee and Dutta¹⁰	Thermal effects	Steady-state Non-isothermal	Comparison with previous model results
2.	Berning and Djilali¹¹	Transport phenomena	Steady-state Non-isothermal	None
3.	Mazumder and Cole¹²	Transport phenomena	Transient Non-isothermal	Comparison with experimental data
4.	Coppo et al.¹³	Thermal effects	Transient Non-isothermal	Comparison with experimental data
5.	Matamoros and Bruggemann¹⁴	Thermal effects	Steady-state Non-isothermal	Comparison with experimental data
6.	Al-Baghdadi and Al-Janabi¹⁵	Geometry/operating condition effects	Steady-state Non-isothermal	Comparison with experimental data
7.	Schwarz and Djilali¹⁶	Transport phenomena	Steady-state Non-isothermal	Comparison with experimental data
8.	Shimpalee et al.¹⁷	Thermal effects	Transient Non-isothermal	Comparison with previous model results
9.	Jang et al.¹⁸	Geometry/operating condition effects	Steady-state Isothermal	Comparison with experimental data
10.	Ren et al.¹⁹	Transport phenomena	Steady-state Isothermal	Comparison with experimental data
11.	Wang²⁰	Transport phenomena	Transient Isothermal	Comparison with experimental data
12.	Wang et al.²¹	Transport phenomena	Transient Non-isothermal	Comparison with experimental data
13.	Berning et al.²²	Transport phenomena	Steady-state Non-isothermal	Comparison with previous model results
14.	Schwarz and Beale²³	Transport phenomena	Steady-state Non-isothermal	None
15.	Schwarz and Djilali²⁴	Geometry/operating condition effects	Steady-state Non-isothermal	Comparison with previous model results
16.	Berning et al.²⁵	Geometry/operating condition effects	Steady-state Non-isothermal	Comparison with experimental data
17.	Yuan et al.²⁶	Geometry/operating condition effects	Steady-state Non-isothermal	Comparison with experimental data
18.	Wu et al.²⁷	Transport phenomena	Transient Non-isothermal	Comparison with experimental data
19.	Cordiner et al.²⁸	Transport phenomena	Steady-state Isothermal	Comparison with experimental data
20.	Fink and Fouquet²⁹	Transport phenomena	Steady-state Non-isothermal	Comparison with experimental data
21.	Kang et al.³⁰	Geometry/operating condition effects	Transient Isothermal	Comparison with experimental data
22.	Chiu et al.³¹	Geometry/ operating condition effects	Steady-state Isothermal	Comparison with experimental data
23.	Sun et al.³²	Transport phenomena	Steady-state Isothermal	Comparison with experimental data
24.	Cao et al.³³	Thermal effects	Steady-state Non-isothermal	Comparison with experimental data
25.	Hossain et al.³⁴	Transport phenomena	Steady-state Isothermal	Comparison with experimental data
26.	Jiang and Wang³⁵	Transport phenomena	Transient Non-isothermal	None
27.	Mancusi et al.³⁶	Geometry/operating condition effects	Steady-state Non-isothermal	Comparison with experimental data
28.	Choopanya and Yang³⁷	Geometry/operating condition effects	Transient Non-isothermal	Comparison with experimental data
29.	Ferreira et al.³⁸	Transport phenomena	Transient Isothermal	Comparison with experimental data
30.	Jo and Kim³⁹	Transport phenomena	Transient Isothermal	None
31.	Khazaee and Sabadbafan⁴⁰	Geometry/operating condition effects	Steady-state Non-isothermal	Comparison with experimental data

Modelling assumptions

Regardless of the type of multiphase flow model being used, it is often necessary to make some modelling assumptions in order to simplify the problem under consideration and thereby facilitate its numerical implementation without reducing the accuracy of the model. Indeed, these assumptions must be based on well-established theories and practices.

Some commonly used assumptions in PEM fuel cell modelling are steady state, ideal gases, laminar and incompressible gas flow, local thermal equilibrium, and isotropic and homogeneous components materials. For the state of water produced by electrochemical reaction in the cathode CL, both vapour and liquid assumptions are used.

Besides these assumptions, there are many other component and operating condition-specific simplifying assumptions that are also used. Some of these may, however, limit the capability of the model or even lead to erroneous results.

An example is the isothermal assumption, which means that the transport of thermal energy given by Table 1, equation (4) is not solved. This can produce results that are not physically representative when phase change is accounted for. In fact, it has been shown that the impact of temperature distribution on the amount of water that undergoes phase change is very significant.⁴¹

Another example is considering the CL as an interface without thickness. Although this assumption greatly simplifies the CL modelling and facilitates its numerical implementation, it may however introduce inaccuracy due to the neglected contribution of ionic resistance, particularly at low humidity.²⁰ Furthermore, experimental results suggest an agglomerate structure for the CL,^42,43 and it has been shown that agglomerate-type models agree better with the physical picture of reactant transport processes in the CL and they best fit the experimental data.^44–46

Also, many models neglect the influences of the motion of liquid water in the GFCs. As pointed out in ‘Components’ section, the GFCs in the BPs deliver reactant gases and help remove product water. This results in gas–liquid two-phase flow in the GFCs. Consequently, the reactant gases transport can be hindered when considerable amount of liquid water accumulate in the channels leading to poor cell performance. Thus, the effects of liquid water in the GFCs should be accounted for.

Models for transport phenomena

Berning and Djilali¹¹ presented a model that accounts for the major transport processes, as well as phase change. The results show that phase change exists in both the anode and cathode sides of the fuel cell. The anode and cathode transports are however decoupled in this model, the CLs are treated as thin interfaces, and liquid water transport by electro-osmotic drag is neglected.

Mazumder and Cole¹² modified a model originally proposed by Wang and Cheng⁴⁷ by adding an additional governing equation for the formation and transport of liquid water. Although the numerical predictions of the cell performance quantitatively over-predict experimentally measured polarization in the high current density regime, the results nonetheless indicate that at critical current density, the saturation levels can exceed 50% and are highest in the cathode.

Schwarz and Djilali¹⁶ incorporated the 1D multiple thin-film agglomerate catalyst model of Baschuk and Li⁴⁸ into a comprehensive 3D model to investigate the effects of transport limitations on cell performance in the CLs. It was found that transport limitations and ohmic losses in the CLs have substantial negative effects on the fuel cell performance. The model however does not consider the effects of liquid water in the GFCs.

Ren et al.¹⁹ developed a gas–liquid two-phase flow and transport model to address the need to account for all major flow and transport phenomena with an emphasis on various factors influencing cell performance. The results show that substituting the air with oxygen and increasing the inlet gas velocity while decreasing the thickness of the membrane and the width of the rib will improve the cell performance. This model assumes isothermal condition and water phase change is only considered in the cathode CL. It has been shown that phase change exists in both the anode and cathode sides of the fuel cell.¹¹

Wang²⁰ investigated multiphase flow, species transport, and electrochemical processes and their interactions using a two-phase flow model but with the assumption of isothermal condition. The results indicate that two-phase flow can occur in both the anode and cathode diffusion media, with the two phases coexisting at low humidity. This model does not however account for water phase change.

Wang et al.²¹ implemented a two-phase flow model to understand three critical issues in the channels: liquid water build-up towards the outlets, saturation spike in the vicinity of flow cross-sectional heterogeneity, and two-phase pressure drop. The results reveal that liquid water builds up quickly at the entrance region of the gas channel with the predicted saturation reaching as high as 20%. Water trapping around the geometrical heterogeneity was also found. This model is limited to a single channel that is not integrated into a full cell model.

Berning et al.²² presented a multiphase model that describes the flux of liquid water through the porous media and into the GFCs. In this model, the overall level of the predicted liquid saturation depends predominantly on the fraction of hydrophilic pores which can be accounted for in the irreducible saturation. This can be a limiting factor for the accuracy of the predictions since the results show that the irreducible saturation has an important impact on the limiting current density. Also, the interaction between the gas and liquid phases in the GFCs is not considered.

Schwarz and Beale²³ implemented a comprehensive model to perform multiphase transport calculations and to investigate the effects of ohmic losses and transport limitations on current density distributions for the cell design. It was found that the current density in the cathode CL is concentrated near the gas channels inlets due to the resistances of air and flooding to oxygen transport. The model however neglects the effects of liquid water in the gas channels.

Wu et al.²⁷ developed a model that fully couples the major transport processes and performed both liquid and vapour production modelling. The results show that the dynamic response of the fuel cell in vapour production modelling is significantly overestimated compared to liquid production modelling. A potential drawback of this model is that it was validated against experimental data that were obtained using a different type of membrane material.

Cordiner et al.²⁸ extended the validity of the existing isothermal two-phase model of Cordiner et al.^49,50 by implementing a decoupled GDL model into a CFD solver to describe the effect of the representation of liquid water flooding on cell performance. It was concluded that flooding must be treated as a 3D phenomenon, as it has different impacts in different regions of the fuel cell.

Fink and Fouquet²⁹ presented a model accounting for the major transport phenomena and their coupling. The simulation results show that liquid water mainly accumulates below the channel lands where condensation takes place resulting in liquid water production. A shortcoming of this model is that it does not study the influence of liquid water droplet movements on performance in the GFCs.

Sun et al.³² introduced Kirchhoff transformation technique and a combined finite element-upwind finite volume approach for efficiently achieving a fast convergence and reasonable solutions for a two-phase transport model. The results demonstrate the efficiency of the numerical technique used. The model used is, however, isothermal and does not take into account water back diffusion effect through the membrane.

Hossain et al.³⁴ developed a two-phase model to investigate the transport of species in the GDL considering liquid water saturation. It was found that higher cell performance can be achieved by optimizing the permeability of the GDLs. The main limits of this model are the assumption of isothermal condition and 1D water transport in the membrane by back diffusion only.

Jiang and Wang³⁵ developed a two-phase flow model by fully coupling species transport, heat transport, and electrochemical processes to investigate the formation and the transport of liquid water. It was found that liquid water build-up in the cathode channel dominates liquid water spreading inside the cell and dictates the performance of the fuel cell.

Ferreira et al.³⁸ used the volume of fluid (VOF) method to numerically investigate two-phase flow in the anode gas channel and analyse water movement. The results reveal that water moves as films for hydrophilic channel walls, whereas it moves as a droplet when the channel walls are hydrophobic. This model is isothermal and phase change or electrochemical reactions are not considered.

Jo and Kim³⁹ also used the VOF method to numerically investigate the dynamics of liquid water emerging from a micro pore on the GDL into the GFC. It was concluded that with decreased GDL surface contact angle, droplets from the outer and inner pores move along the side walls and the movement of droplets from the centre pore shows complex patterns. Furthermore, as the hydrophobicity of the side and top walls increases, the GDL surface water coverage ratio increases while the water volume fraction decreases. This model is isothermal and phase change or electrochemical reactions are not considered.

Models for geometry or operating condition effects

Al-Baghdadi and Al-Janabi¹⁵ conducted an optimization study using a multiphase model that incorporates the significant physical processes and the key parameters affecting cell performance. The results show that the model is capable of accurately quantifying the impacts of operating, design, and material parameters on the fuel cell performance. The main limits of this model are the treatment of the CL as a thin interface and the consideration of phase change in the GDLs only.

Jang et al.¹⁸ developed a model with conventional flow-field designs to study the influences of flow-field designs on fuel utilization, water removal, and cell performance. It was found that the cell performance of the serpentine flow field is the best followed by the Z-type flow field and the parallel flow field. The model assumes isothermal condition and does not account for the motion of liquid water.

Schwarz and Djilali²⁴ used the multiphase flow model of Schwarz and Djilali¹⁶ to explore the implementation of spatially distributed catalyst loadings in all three directions for both the anode and the cathode. The results show that at high current densities, the over-potential variations in the anode and the cathode CLs are fairly large. This study does not examine the use of variable composition for the CL.

Berning et al.²⁵ extended a previously published model by Berning et al.²² in order to account for phase change kinetics and to conduct a comparison between the interdigitated flow-field design and a conventional straight channel design. It was found that the interdigitated flow field yields higher and more uniform oxygen concentration as well as a lower overall liquid saturation at the CL. The interaction between the gas and liquid phases is, however, not considered in this model.

Yuan et al.²⁶ developed a multiphase model to predict the effects of operating parameters such as operating pressure, cell temperature, relative humidity of reactant gases, and air stoichiometric ratio on cell performance. The results indicate that increased operating pressure and temperature can enhance the cell performance. Moreover, it was found that the best performance occurs at moderate air relative humidity while hydrogen is fully humidified.

Kang et al.³⁰ used a multiphase porous cathode side model to study liquid water flooding in an interdigitated flow-field design. The results reveal the existence of a so-called liquid water avalanche phenomenon in this type of cathode design. Also, the existence of three distinct phases for liquid water flooding process is found; these are porous layer phase, channel phase, and drainage phase. This model is limited to the cathode side only and heat transfer and electrochemical reactions are not considered. Furthermore, the model does not account for water phase change.

Chiu et al.³¹ presented a two-phase transport model with parallel, interdigitated, and serpentine flow-field designs to investigate the effects of flow-field design on cell performance and water removal. The results show that the cell performance with parallel and serpentine flow fields can be improved by reducing the height or the width of the channel; of the three flow fields, the interdigitated flow field is the most effective on liquid water removal. Potential drawbacks of this model are the assumption of isothermal condition and the use of a highly simplified multiphase model.

Mancusi et al.³⁶ used a combined multi-fluid and VOF approach to study two-phase flow in a tapered channel design. The results show that tapering the channel downstream enhances water removal due to increased airflow velocity.

Choopanya and Yang³⁷ developed a multiphase model to investigate the effect of the parallel and the serpentine flow-field geometries on the cell performance under both steady-state and transient operations. It was found that in the steady-state run, the average current density increases with the use of a serpentine flow field due to superior water removal ability. For the transient operation, the use of a serpentine flow field combined with dry reactant gases help the removal of product water and accelerate the transport of reacting species to the reaction site. This model does not however account for water phase change, and the dynamics and movement of liquid water are not studied.

Khazaee and Sabadbafan⁴⁰ developed a two-phase model to investigate water management and the performance of PEM fuel cell with 1-serpentine and 4-serpentine channels flow-field configurations at different operating conditions. The results show that under identical conditions, the 4-serpentine configuration outperforms 1-serpentine configuration. However, the effects of liquid water in the channels are not studied in detail.

Models for thermal effects

Shimpalee and Dutta¹⁰ modified a model already presented by Dutta et al.⁵¹ by including the energy equation to predict the temperature distribution inside a straight channel and the effect of heat produced by the electrochemical reactions on cell performance. The predictions show that the cell performance depends not merely on geometric or operating condition parameters but also on the temperature rise inside the fuel cell. A shortcoming of this model is that it neglects the heat generated by ohmic resistance in the porous media and the membrane.

Coppo et al.¹³ used the 3D implementation of a previously developed 2D model by Siegel et al.^44,52 to analyse the effects of temperature on the cell operation. It was found that higher temperature positively affects the reaction kinetics. However, only the heat source due to phase change in the cathode CLs is added in the energy equation.

Matamoros and Bruggemann¹⁴ developed a model that specifically computes water and heat management. The simulation results show that the operating temperature has negative effects on the hydrating properties of the membrane. The only heat sources included in the model are Joule heating in the membrane and latent heat in the cathode.

Shimpalee et al.¹⁷ extended the existing model of Shimpalee et al.^53,54 to include energy and water phase change equations for studying phase change and heat transfer. It was found that a single-phase and isothermal model shows higher current overshoot than a water phase change and non-isothermal model. The model does not account for the heat sources due to ohmic resistance in the porous media and the membrane.

Cao et al.³³ presented a two-phase model to investigate amongst other things the interaction between water and thermal transport processes, and different boundary temperature of flow plate. The results show that the boundary temperature greatly affects the cell temperature distribution and indirectly influences water saturation distribution. The model however neglects liquid saturation in the channels.

Model validation

Validation is an essential aspect of any CFD-based analysis. During the validation process, the simulation results are compared with known experimental data or the results from previous numerical models.

Almost all of the models listed in Table 2 have been validated using some kind of experimental data. These include a detailed comparison of the predicted polarization curve or V–I curve with experimentally measured one. A comparison of the simulation results with the results from previous models was also performed in a few cases. Overall, the results were found to be in good agreement with either experimental or previous model results.

Multiphase flow models and phase change

As described in ‘Multiphase flow’ section, in PEM fuel cells, the production of liquid water by the oxygen reduction reaction and the phase change process are the origin of two-phase transport. The gas phase has multi-components (hydrogen and water vapour in the anode side; oxygen, nitrogen, and water vapour in the cathode side), while the liquid phase consists of only pure water.

In the reviewed articles, multiphase flow is modelled using the multi-fluid and multiphase mixture approaches essentially although the VOF model as well as some highly simplified models such as the moisture diffusion model, the porosity correction model, and a model in which the liquid phase is treated as a component of the gas phase is also used. Among other things, the multi-fluid model requires the specification of phase change kinetics, whereas local equilibrium between the gas phase and the liquid phase is inherently assumed in the multiphase mixture model. For more details on these topics, the reader is referred to the reference texts by Sunden and Faghri⁵⁵ and Wang et al.⁷

In this section, a brief overview of two-phase flow modelling is presented with a summarized outline of the selected multiphase flow models given in Table 3. The main water phase change mechanisms used are condensation of water vapour to liquid water and evaporation of liquid water to water vapour.

Table 3.

Summary of multiphase flow models and implementation of phase change.

Ref.	Multiphase flow model	Implementation of phase change (evaporation/condensation)
Shimpalee and Dutta¹⁰	Liquid phase is a component of gas phase	$S_{{wv}_{p}} = \frac{M_{H_{2} O} \sum_{n of v} {mass}_{n of v} / M_{n of v}}{(I - \frac{P_{wv}^{sat}}{P})} [\frac{P_{wp}^{sat} - P_{wv}}{P}] * C$
Berning and Djilali¹¹	Multi-fluid	${\overset{\cdot}{N}}_{W} = ψ N_{D, CV} k_{xm} π D^{2} \frac{x_{w 0} - x_{w \infty}}{1 - x_{w 0}}, RH \leq 100 %$ ${\overset{\cdot}{m}}_{cond} = - ψ C \frac{x_{w 0} - x_{w \infty}}{1 - x_{w 0}}, RH > 100 %$
Mazumder and Cole¹²	Multiphase mixture	${\overset{\cdot}{m}}_{l} = {M_{l} k_{c} \frac{ɛ x_{W}}{RT} (x_{W} P - P_{sat}) if x_{W} P > P_{sat} M_{l} k_{e} \frac{ɛ s ρ_{l}}{M_{l}} (x_{W} P - P_{sat}) if x_{W} P < P_{sat}$
Coppo et al.¹³	Multi-fluid	$S_{LV} = ψ s γ_{LV} - ψ (1 - s) (1 - γ_{LV})$
Matamoros and Bruggemann¹⁴	Multi-fluid	${\overset{\cdot}{n}}_{{GLH}_{2} O} = K_{GL} M_{H_{2} O} \frac{A_{LG}}{V_{L}} \frac{P_{H_{2} O (G)} - P_{H_{2} O}^{Sat}}{RT}$
Al-Baghdadi and Al-Janabi¹⁵	Multi-fluid	${\overset{\cdot}{m}}_{evap} = M_{H_{2} O} ϖ N_{D} k_{xm} π D_{drop} \frac{x_{w 0} - x_{w \infty}}{1 - x_{w 0}}, RH \leq 100 %$ ${\overset{\cdot}{m}}_{cond} = ϖ C \frac{x_{w 0} - x_{w \infty}}{1 - x_{w 0}}, RH > 100 %$
Schwarz and Djilali¹⁶	Multi-fluid	$\nabla \cdot j_{w} = Cmax {(1 - s) \frac{x_{H_{2} O} P - P_{sat}}{RT} M_{H_{2} O}, - s ρ_{w}}$
Shimpalee et al.¹⁷	Liquid phase is a component of gas phase	$S_{{wv}_{p}} = \frac{M_{H_{2} O} \sum_{n of v} {mass}_{n of v} / M_{n of v}}{1 - (P_{wv}^{sat} / P)} [\frac{P_{wp}^{sat} - P_{wv}}{P}] * C$
Jang et al.¹⁸	Porosity correction	$S_{L} = {{Mk}_{c} \frac{ɛ_{eff} C_{H_{2} O}}{ρ RT} (P_{H_{2} O} - P_{sat}) if P_{H_{2} O} > P_{sat} k_{e} ɛ_{eff} s (P_{H_{2} O} - P_{sat}) if P_{H_{2} O} < P_{sat}$
Ren et al.¹⁹	Liquid phase is a component of gas phase	$S_{wvp} = - \frac{(p_{wv} - p_{wv}^{sat})}{RT} \times M_{H_{2} O} \times k_{c}$
Wang²⁰	Multiphase mixture	None
Wang et al.²¹	Multiphase mixture	${\overset{\cdot}{m}}_{fg} = ρ_{l} \frac{\partial s}{\partial t} + \nabla \cdot (ρ_{l} {\vec{u}}_{l})$
Berning et al.²²	Multi-fluid	None
Schwarz and Beale²³	Multi-fluid	$S_{w} = Cmax {(1 - s) \frac{P_{H_{2} O} - P_{sat}}{RT} M_{H_{2} O}, - s ρ_{w}}$
Schwarz and Djilali²⁴	Multi-fluid	$\nabla \cdot j_{w} = Cmax {(1 - s) \frac{x_{H_{2} O} P - P_{sat}}{RT} M_{H_{2} O}, - s ρ_{w}}$
Berning et al.²⁵	Multi-fluid	${\overset{\cdot}{m}}_{H_{2} O, evap} = k_{xm} π D^{2} M_{H_{2} O} \frac{x_{A 0} - x_{A \infty}}{1 - x_{A 0}}, RH < 100 %$ ${\overset{\cdot}{m}}_{H_{2} O, cond} ″ \frac{= 1}{4} {\bar{u}}_{m} α_{m} M_{H_{2} O} \frac{p_{g, H_{2} O} - p_{H_{2} O} (T)}{R_{g} T}, RH > 100 %$
Yuan et al.²⁶	Multiphase mixture	$r_{w} = C_{r} \max [(1 - s) \frac{p_{WV} - p_{sat}}{RT} M_{W, H_{2} O}, - s ρ_{l}]$
Wu et al.²⁷	Multi-fluid	$S_{vl} = {A_{pore} \frac{{Sh}_{c} D_{w}}{d} (1 - s) (ρ_{v} - ρ_{sat}), if ρ_{v} \geq ρ_{sat} A_{pore} \frac{{Sh}_{e} D_{w}}{d} (ρ_{v} - ρ_{sat}), if ρ_{v} < ρ_{sat}$
Cordiner et al.²⁸	Multi-fluid	${\overset{\cdot}{m}}_{gas - to - liquid} = γ M_{W, H_{2} O} \frac{(P_{V, H_{2} O} - P_{sat, H_{2} O} (T))}{R_{g} T}$
Fink and Fouquet²⁹	Multi-fluid	${\overset{\cdot}{M}}_{pc, H_{2} O} = \frac{A_{\lg}}{V} {\bar{ρ}}_{g} {\bar{D}}_{H_{2} O, g} \frac{{\tilde{Sh}}_{H_{2} O, g}}{d} \ln (1 - Bm)$
Kang et al.³⁰	Multiphase mixture and volume of fluid	None
Chiu et al.³¹	Porosity correction	$S_{L} = {M_{w} k_{c} \frac{ɛ_{eff} C_{W}}{ρ RT} (P_{w} - P_{sat}) if P_{W} > P_{sat} k_{e} ɛ_{eff} s (P_{W} - P_{sat}) if P_{W} < P_{sat}$
Sun et al.³²	Multiphase mixture	None
Cao et al.³³	Moisture diffusion	$S_{l} = A_{pore} \frac{{Sh}_{c} D_{w}}{\bar{d}} (1 - s) (ρ_{w} - ρ_{sat}) q +$ $A_{pore} \frac{{Sh}_{e} D_{w}}{\bar{d}} s (ρ_{w} - ρ_{sat}) (1 - q)$
Hossain et al.³⁴	Liquid phase is a component of gas phase	$S_{wlp} = - S_{wvp} = \frac{(p_{wv} - p_{wv}^{sat})}{RT} \times M_{H_{2} O} \times k_{c}$
Jiang and Wang³⁵	Multiphase mixture	${\overset{\cdot}{m}}_{gl} = ρ_{l} \frac{\partial s}{\partial t} + \nabla \cdot (ρ_{l} {\vec{u}}_{l})$
Mancusi et al.³⁶	Multi-fluid and volume of fluid	$r_{w} = C_{r} \max (((1 - s) \frac{P_{wv} - P_{sat}}{RT} M_{H_{2} O}), (- s ρ_{L}))$
Choopanya and Yang³⁷	Porosity correction	$r_{w} = C_{r} \max ([(1 - s) \frac{p_{WV} - p_{sat}}{RT} M_{W, H_{2} O}], [- s ρ_{l}])$
Ferreira et al.³⁸	Volume of fluid	None
Jo and Kim ³⁹	Volume of fluid	None
Khazaee and Sabadbafan ⁴⁰	Multiphase mixture	$S_{pc} = {\overset{\cdot}{m}}_{pc} = k_{cond} (1 - s_{l}) \frac{M_{H_{2} O}}{RT} P (X_{g}^{H_{2} O} - X_{sat}^{H_{2} O}) if X_{g}^{H_{2} O} \geq X_{sat}^{H_{2} O}$ $S_{pc} = {\overset{\cdot}{m}}_{pc} = k_{evap} s_{l} ρ_{l} P (X_{g}^{H_{2} O} - X_{sat}^{H_{2} O}) if X_{g}^{H_{2} O} \leq X_{sat}^{H_{2} O}$

Multi-fluid model

The multi-fluid approach of multiphase flow modelling consists of solving a set of conservation equations for mass and momentum for each phase separately. The two phases are usually coupled through the relative permeability and the phase change terms. The following are generic steady-state mass and momentum conservation equations in the porous media

Mass

\nabla \cdot (ɛ ρ_{g} {\vec{V}}_{g}) = S_{PC}

(10)

\nabla \cdot (ɛ ρ_{l} {\vec{V}}_{l}) = S_{PC}

(11)

Momentum

{\vec{V}}_{g} = - (1 - s) \frac{K_{g}}{μ_{g}} \nabla P_{g}

(12)

{\vec{V}}_{l} = - s \frac{K_{l}}{μ_{l}} \nabla P_{g} - D (s) \nabla s

(13)

The multi-fluid approach best suits high saturation conditions, as these demand a greater liquid resolution. The model can resolve complex liquid motion and accounts for the convection of liquid by the gas. As for the disadvantages, the multi-fluid model requires a high number of variables and the coupling of the phases can lead to unstable models.⁵⁵

Multiphase mixture model

As its name suggests, in the multiphase mixture model, the phases are considered to be mixed; and thus, a single set of conservation equations for the phase mixture is solved. Phase equilibrium is assumed and the mixture quantities such as density and so forth, as well as the relative velocity among different phases are evaluated subsequently. The mass conservation and momentum equations for the mixture read⁷

Mass

\nabla \cdot (ρ \vec{V}) = 0

(14)

Momentum

\vec{V} = - \frac{K}{ρ μ} \nabla P

(15)

The mixture density and velocity are defined as⁷

ρ = s ρ_{l} + (1 - s) ρ_{g}

(16)

and

ρ \vec{V} = {\vec{V}}_{l} ρ_{l} + {\vec{V}}_{g} ρ_{g}

(17)

The mixture species conservation equation is expressed as⁷

\nabla \cdot (γ_{k} \vec{V} c_{k}) = \nabla \cdot (D_{k g, eff} \nabla c_{k g}) - \nabla \cdot [(\frac{w_{k l}}{M_{k}} - \frac{c_{k g}}{ρ_{g}}) {\vec{j}}_{l}] + S_{k}

(18)

where the convection correction factor γ is a function of the liquid saturation.

The multiphase mixture approach is best used when the gas phase pressure dominates the liquid phase pressure or when high capillary numbers are encountered. This approach reduces the number of variables and can effectively model the influence of the gas pressure on the liquid. However, large number of mixture quantities need to be calculated and the model may have convergence issues at higher saturations.⁵⁵

Moisture diffusion model

The moisture diffusion model is based on the unsaturated flow theory. It is used to determine the transport of liquid water when the only driving force is the capillary pressure.

The transport equation of liquid water in the moisture diffusion model is written as⁵⁵

\nabla \cdot (D (s_{l}) \nabla s_{l}) + {\overset{\cdot}{S}}_{l} = 0

(19)

where the capillary diffusion coefficient D(s_l) and the mass source due to phase change

{\overset{\cdot}{S}}_{l}

are given by⁵⁵

D (s_{l}) = \frac{ρ_{l} k_{l}}{μ_{l}} \frac{\partial p_{c}}{\partial s_{l}}

(20)

{\overset{\cdot}{S}}_{l} = - {\overset{\cdot}{S}}_{g, H_{2} O}

(21)

The moisture diffusion model is well suited for low capillary numbers and when surface tension is the dominant force on the liquid. It only requires one additional equation to an existing single-phase model. A limitation of this approach is that the influence of the gas pressure on the liquid is not accounted for.⁵⁵

Porosity correction model

In the porosity correction model, the transport of liquid water is not modelled. Instead, it is assumed that liquid water fully occupy some pores of the porous media, meaning that no pore presents coexistence of gas and liquid water. Thus, the volume fraction open to the gas phase is determined by⁵⁵

ɛ_{g} = ɛ (1 - s)

(22)

This approach is best used in low relative humidity, very small pores, and low current density conditions. It does not require any additional transport equations over the single-phase model. A major drawback of this model is that it does not account for the motion of liquid water.⁵⁵

VOF model

VOF is a method for locating and advecting the interface in a gas–liquid two-phase flow. The method proceeds by reconstruction of the interface shape and then advection of the reconstructed interface in a given velocity field.⁵⁶ A single momentum equation such as the one given by Table 1, equation (2) is solved and the volume fraction of each of the fluids is tracked throughout the domain.

The VOF method can capture the effects of surface tension making it well suited for micro-channel flows in which surface tension is an important and sometimes dominant force. A limitation of the VOF technique for multiphase flow modelling is that it only deals with the interface of the liquid and gas. The entire fuel cell multiphase flow computation is still dealt with by a CFD package using the multi-fluid or the multiphase mixture model.

Phase change

In a thermodynamic system, the state of matter can change from one phase to another. This transformation process is known as phase change.

In PEM fuel cells, water may exist as three different phases: water vapour, liquid water, and dissolved water. Since this type of fuel cell operates at relatively low temperatures (50–100℃), phase change is thus frequently encountered.

In fact, liquid water present within the cell can evaporate into water vapour, and conversely, water vapour can condense into liquid water. Also, phase change occurs from water vapour to dissolved water (water uptake into the polymer phase of the CL) and from dissolved water to liquid water (water released from the polymer phase of the CL). In many phase change models, the occurrence of water phase change is determined by the gas-phase humidity.

In the selected papers, phase change kinetics is mainly implemented in the form of evaporation and condensation. The total amount of water undergoing evaporation depends on the level of undersaturation and the surface area of the liquid phase. Condensation is fundamentally different from evaporation as it can occur on any hydrophilic surface. Water vapour condenses into liquid water droplets when the local water vapour pressure exceeds the local saturation pressure.

In the majority of the articles reviewed, dissolved water resulting from the transport of protons in the CLs and the membrane is not considered as a separate water phase. Consequently, dissolved water related-phase change processes are neglected; dissolved water is simply treated as part of the liquid phase.

However, since the protons produced in the CL are assumed to be dissolved in water and are transported along with the water molecules through diffusion, due to a concentration gradient between the anode and the cathode, and electro-osmotic drag; dissolved water should be considered as a separate water phase, and thus dissolved water-related phase change processes should be treated separately. The governing equation of dissolved water transport is written as⁵²

D_{DW, eff} (\nabla \cdot \nabla c_{DW}) + \frac{n_{d}}{F} σ_{i, eff} (\nabla \cdot \nabla Φ_{i}) = S_{PC}

(23)

The total rate of water vapour transfer into dissolved water is calculated as follows⁴⁶

S_{vapor to dissolved} = κ_{adso} (c_{W}^{*} - c_{DW}) if c_{DW} < c_{W}^{*}

(24)

The total rate of dissolved water transfer into liquid water is calculated as follows⁴⁶

S_{dissolved to liquid} = κ_{deso} (c_{DW} - c_{LW}^{*}) if c_{DW} > c_{W}^{*}

(25)

where

c_{W}^{*}

represents the equilibrium water concentration and κ_adso and κ_deso are water adsorption and desorption mass transfer coefficient, respectively.

Numerical procedures

In this section, salient numerical features of the reviewed articles, such as system boundary, computational domain, and CFD solvers used are discussed. A summary of this information is given in Table 4.

Table 4.

Summarized numerical features of the reviewed papers.

Ref.	System boundary	Computational domain	CFD software
Shimpalee and Dutta¹⁰	Single flow channel with GDL and MEA	Single domain	Fluent
Berning and Djilali¹¹	Single flow channel with GDL and MEA	Multi-domain	CFX
Mazumder and Cole¹²	Single flow channel with GDL and MEA	Single domain	CFD-ACE+
Coppo et al.¹³	Single flow channel with GDL and MEA	Single domain	CFDesign
Matamoros and Bruggemann¹⁴	Single flow channel with GDL and MEA	Single domain	Own written codes
Al-Baghdadi and Al-Janabi¹⁵	Single flow channel with GDL and MEA	Multi-domain	Unspecified
Schwarz and Djilali¹⁶	Single flow channel with GDL and MEA	Single domain	Fluent
Shimpalee et al.¹⁷	Complete cell	Single domain	STAR-CD
Jang et al.¹⁸	Complete cell	Single domain	Unspecified
Ren et al.¹⁹	Double flow channels with GDL and MEA	Single domain	Unspecified
Wang²⁰	Single flow channel with GDL and MEA	Single domain	Own written codes
Wang et al.²¹	Single flow channel	Single domain	Fluent
Berning et al.²²	Single flow channel with GDL and CL	Multi-domain	unspecified
Schwarz and Beale²³	Complete cell	Single domain	Fluent
Schwarz and Djilali²⁴	Single flow channel with GDL and MEA	Single domain	Fluent
Berning et al.²⁵	Single flow channel with GDL and MEA	Multi-domain	CFX
Yuan et al.²⁶	Single flow channel with GDL and MEA	Single domain	Fluent
Wu et al.²⁷	Single flow channel with GDL and MEA	Single domain	Fluent
Cordiner et al.²⁸	Single flow channel with GDL and MEA	Multi-domain	Fluent
Fink and Fouquet²⁹	Complete cell	Multi-domain	AVL FIRE
Kang et al.³⁰	Complete cell	Single domain	Fluent
Chiu et al.³¹	Single flow channel with GDL and MEA	Single domain	CFD-ACE+
Sun et al.³²	Single flow channel with GDL and MEA	Multi-domain	Own written codes
Cao et al.³³	Single flow channel with GDL and MEA	Single domain	Own written codes
Hossain et al.³⁴	Single flow channel with GDL and MEA	Single domain	Fluent
Jiang and Wang³⁵	Single flow channel with GDL and MEA	Single domain	Fluent
Mancusi et al.³⁶	Single flow channel with GDL and MEA	Single domain	Fluent
Choopanya and Yang³⁷	Complete cell	Single domain	Fluent
Ferreira et al.³⁸	Single flow channel	Single domain	Fluent
Jo and Kim³⁹	Single flow channel	Single domain	Fluent
Khazaee and Sabadbafan⁴⁰	Complete cell	Single domain	Fluent

CFD: computational fluid dynamics; GDL: gas diffusion layer; MEA: membrane electrode assembly.

System boundary

In PEM fuel cell modelling, defining a manageable system boundary in order to reduce computational efforts without compromising accuracy is often necessary.

From Table 4, it is very apparent that the vast majority of the models analysed take advantage of the symmetric layout of the fuel cell geometry to only model a single GFC on both anode and cathode sides with the corresponding GDLs and CLs and the membrane as illustrated in Figure 4. This reduces the computational efforts considerably. A full 3D description of a complete cell or a fuel cell stack would require very large computing resources and overly long simulation time. Thus, a complete cell is only computed in five cases (in the models for geometry or operating condition effects mainly) and the simulation of a complete fuel cell stack is not performed.

Figure 4.

A single straight channel with corresponding GDL and MEA. GDL: gas diffusion layer; MEA: membrane electrode assembly.

Computational domain

In terms of modelling domain (computational domain), two distinct approaches are used in the reviewed papers for computing the numerical solutions of PEM fuel cell mathematical models according to the different regions within the cell. These are single-domain and multi-domain approaches, with a greater number of the models reviewed adopting the single-domain approach (see Table 4).

In the single-domain approach, the generic conservation law equations are used for the entire modelling domain and only the source terms vary according to each region or component, whereas in the multi-domain approach different forms of the same equations are solved in each region or component within the cell, and boundary or initial conditions may be required on internal interfaces.

The advantage to the single-domain approach is that it does not require the specification of interface conditions between the individual fuel cell components and aid in parallel computing. It may however involve a careful numerical manipulation of the equations to retain a model that is physically realistic since all of the equations are solved at once over the entire model which can also be wasteful in terms of memory and speed.

The multi-domain approach can efficiently deal with decoupled components. It is however inefficient when coupling dominates and it does not readily facilitate parallelization.

Solvers

As explained in ‘CFD’ section, for a CFD solver to solve a non-linear differential model, this must first be converted into a system of linear algebraic equations. In most CFD software, this conversion is performed by means of the finite volume method (FVM), which discretizes both the computational domain and the conservation equations using highly sophisticated algorithms.

For the coupling of velocity and pressure, the majority of the available CFD codes use the semi-implicit method for pressure-linked equations algorithm although the pressure-implicit split-operator algorithm is also used.

Based on the literature model overview, one can say that Fluent is well suited for 3D PEM fuel cell modelling. More details on the most frequently used CFD codes for PEM fuel cell simulation is provided in ‘CFD codes for PEM fuel cell modelling’ section.

CFD codes for PEM fuel cell modelling

There are several CFD codes that can be used for PEM fuel cell model simulation ranging from commercial software through open source free codes to own written programs.

Although most of these software use similar techniques for the discretization of both the computational domain and the partial differential equations (e.g. FVM) and have similar capabilities to deal with the multidimensional effects of the phenomena taking place in the fuel cell, they occasionally differ in their approach to modelling multiphase flow and the phase change processes.

This section gives a brief overview of the most commonly used CFD software for the numerical modelling and simulation of PEM fuel cells.

Commercial codes

AVL FIRE™

AVL FIRE is a multi-purpose thermo-fluid CFD software package that is capable of simulating fluid dynamics problems involving complex geometries and the interplay of advanced physics and chemistry.⁵⁷

The multiphase flow modules consist of Eulerian and Lagrangian multiphase modules. The Eulerian multiphase module uses the multi-fluid modelling approach of multiphase flow and thus allows the calculation of the volume fraction distribution for each phase in addition to all other flow variable. The Lagrangian multiphase module accounts for droplet break-up, turbulence dispersion, collision and coalescence, distortion and drag, evaporation as well as droplet/wall interaction. It thus allows a detailed study of the movements of the droplet.⁵⁷

CFD-ACE+

CFD-ACE+ is an advanced CFD and multi-physics software that enables coupled simulations of fluid, thermal, chemical, etc. phenomena. It is specifically designed for parallel computing on high performance workstations and clusters but it also works on normal computer systems.⁵⁸

The software provides built-in models for the modelling of electrochemistry and flow through porous media and small channels. The PEM fuel cell module accounts for the fundamental physics of PEM fuel cell. It includes a model for water transport through the membrane, a model for transport of liquid water saturation through the porous media, and the liquid saturation model for two-phase flow in the channels.⁵⁸

CFX

ANSYS CFX is a general purpose fluid dynamics program that offers an abundant choice of physical models, capturing any type of phenomena related to fluid flow. Its highly parallelized solver supports the simulation of flow in porous media, multiphase flow, phase change, and so forth.⁵⁹

Two distinct multiphase flow models are provided, Eulerian–Eulerian multiphase model and a Lagrangian particle tracking multiphase model. The Eulerian–Eulerian multiphase model allows the modelling of multiphase flow using either the multiphase mixture or the multi-fluid modelling approaches. The Lagrangian particle tracking multiphase model models the total flow of the particle phase by tracking a small number of particles through the continuum fluid. The particles could be solid particles, drops, or bubbles.⁶⁰

COMSOL Multiphysics®

COMSOL Multiphysics is a very versatile software platform for various coupled phenomena.⁶¹

A battery and fuel cell module that can be used to model the underlying electrochemical phenomena in the electrodes and electrolytes of batteries and fuel cells is provided. The CFD module allows the modelling of multiphase flows as well as flow in porous media. Multiphase flow is modelled using the bubbly flow, multiphase mixture, or Eulerian–Eulerian multiphase models. Phase changes are described using built-in step functions.⁶¹

Fluent

ANSYS Fluent is a general purpose fluid analysis software package with broad physical modelling capabilities for modelling flow, turbulence, heat transfer, and reactions for industrial applications.⁶²

It provides a fuel cell module that can be used for the modelling of electrochemistry, current and mass transport, heat source, liquid water formation, and transport in PEM fuel cells.⁶³

Two distinct approaches are used for the numerical calculation of multiphase flows: the Eulerian–Lagrangian approach and the Eulerian–Eulerian approach. In the Eulerian–Lagrangian approach the fluid phase is treated as a continuum by solving the Navier–Stokes equations, while the dispersed phase is solved by tracking a large number of particles, bubbles, or droplets through the calculated flow field. In the Eulerian–Eulerian approach, the different phases are treated mathematically as inter-penetrating continua. Three different Eulerian–Eulerian multiphase models are available in Fluent: VOF model, the mixture model, and the Eulerian model.⁶⁴

STAR-CD®

CD-adapco STAR-CD and STAR-CCM+ are well-established software platforms for industrial CFD simulation that are capable of solving complex multidimensional and multi-physics models.

STAR-CD has been used with a specialized PEM fuel cell module (es-pemfc) to gain a better understanding of water management as well as electrochemistry in PEM fuel cell, leading to optimized cell performance.⁶⁵

STAR-CCM+ can perform comprehensive simulations involving electrochemistry allowing fuel cell modelling. Multiphase flow modelling is completed using Eulerian multiphase model, with other models within the Eulerian framework, such as the Eulerian multiphase mixture model, VOF, fluid film, and the dispersed multiphase model.⁶⁶

Free codes

OpenFOAM®

OpenFOAM is an open source CFD code that can solve anything from complex fluid flows involving chemical reactions, turbulence, and heat transfer, to solid dynamics and electromagnetics.⁶⁷

The object-oriented design of the software which is written in C++ allows the implementation of own models and numerical algorithms. Thus, the user has complete freedom to customize and extend any existing functionality.⁶⁷

A wide range of solvers is provided for the simulation of flow in porous media and multiphase flows. The approaches available for the modelling of multiphase flows range from a system of two fluid phases model with one phase dispersed, through VOF phase fraction-based interface capturing approach, to multiphase mixture and multi-fluid models.⁶⁷

Other codes

Besides the conventional CFD codes outlined above, there exist other CFD codes which are also capable of performing the same numerical modelling and simulation of PEM fuel cells, and often with equal results. These are usually highly specialized own written programs developed by individual academic or industrial research groups. Multiphase flow is modelled using either the multiphase mixture or the multi-fluid modelling approaches mainly.

Summary and outlook

A comprehensive review of the open literature 3D multiphase CFD models for PEM fuel cell was conducted. This revealed the large volume of published work dealing with the modelling and simulation of transport phenomena in PEM fuel cells, as well as thermal management and the influence of geometry or operating conditions on the cell performance.

It is difficult to directly compare the published models because of the modelling of multiphase flow, phase change, and other modelling assumptions, as well as model validation methods used by the authors. Nevertheless, it appears that 1D and 2D models are only suitable when focusing on particular phenomena, whereas 3D models allow multi-physics simulations of a single component, a complete cell, or even a fuel cell stack and they can be used for design purposes.

Thus, the 3D PEM fuel cell models presented in a selection of journal papers were thoroughly analysed. This revealed that comprehensive 3D multiphase flow CFD models that elucidate the complex couplings of electrochemical kinetics and multiphase transport are still greatly in need though significant advances have been made in PEM fuel cell modelling.

Indeed, truly functional and predictive capabilities remain a challenge due to deficiencies in models for two-phase transport in the GFCs and the porous media. A further limitation of many of the existing comprehensive 3D models is that they do not consider the detailed composition and structure of the CLs, despite the important role of the CLs in determining the cell performance.

Therefore, more research is required to develop models that include the detailed composition and structure of the CLs. Such comprehensive models should also consider the effects of the movement of water droplet in the GFCs, the existence of phase change in both the anode and cathode sides of the fuel cell, as well as dissolved water transport and related phase change mechanisms. These matters have thus far received little attention.

As for the CFD codes used for PEM fuel cell modelling and simulation, it has been noted that the market is largely dominated by a few major brands, such as Fluent, CFD-ACE+, CFX, STAR-CD, COMSOL, AVL FIRE, and OpenFOAM. Although these software packages differ in appearance and capabilities, fundamentally they are all the numerical solvers of partial differential equations with attached multi-physics models and in some cases a specialized fuel cell module.

Finally, in view of all that has been mentioned so far, one may suppose that optimizing PEM fuel cell designs to provide better thermal and water managements, while reducing the costs associated with the testing and operation of PEM fuel cell remains a key technological challenge for this type of fuel cell. Thus, it is likely that in the future, with further advances in numerical procedures and computational techniques, CFD modelling and simulation will continue to play an important role in tackling this technological challenge for the PEM fuel cell technology.

Footnotes

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors acknowledge the financial support from the International Doctoral Innovation Centre, Ningbo Education Bureau, Ningbo Science and Technology Bureau, China’s MoST, and the University of Nottingham. This work was partially supported by Ningbo Natural Science Foundation Program (project code 2013A610107).

Appendix

References

Liu

. Fuel cell performance, New York: Nova Science Publishers, 2012, pp. 284.

Bird

Stewart

Lightfoot

. Transport phenomena, 2nd ed. New York; Chichester: John Wiley & Sons, 2007.

Barbir F and ScienceDirect. PEM fuel cells [electronic resource]: theory and practice. Sustainable world series. Amsterdam; London: Elsevier Academic. 1 online resource (xv, 433 pages), 2005.

Zikanov

. Essential computational fluid dynamics [electronic resource], Hoboken, NJ: Wiley, 2010.

Ramshaw

. Elements of computational fluid dynamics 2011; Vol. 2: London: Imperial College Press.

Hoogers

. Fuel cell technology handbook, Boca Raton, FL; London: CRC Press. 1 v. (various pagings), [16] p. of plates, 2003.

Wang

Chen

Cho

. PEM fuel cells: thermal and water management fundamentals, New York: Momentum Press, 2013.

Basu

. Recent trends in fuel cell science and technology, New York: Springer, 2007.

Bago͡skĭ

. Fuel cells [electronic resource]: problems and solutions, 2nd ed. The Electrochemical Society series. Hoboken, NJ: Wiley. 1 online resource (xviii, 385 p.), 2012.

10.

Shimpalee

Dutta

. Numerical prediction of temperature distribution in PEM fuel cells. Numer Heat Transfer Part A 2000; 38: 111–128.

11.

Berning

Djilali

. A 3D, multiphase, multicomponent model of the cathode and anode of a PEM fuel cell. J Electrochem Soc 2003; 150: A1589–A1598.

12.

Mazumder

Cole

. Rigorous 3-d mathematical modeling of PEM fuel cells – II. Model predictions with liquid water transport. J Electrochem Soc 2003; 150: A1510–A1517.

13.

Coppo

Siegel

von Spakovsky

. On the influence of temperature on PEM fuel cell operation. J Power Sources 2006; 159: 560–569.

14.

Matamoros

Bruggemann

. Simulation of the water and heat management in proton exchange membrane fuel cells. J Power Sources 2006; 161: 203–213.

15.

Al-Baghdadi

Al-Janabi

. Modeling optimizes PEM fuel cell performance using three-dimensional multi-phase computational fluid dynamics model. Energy Convers Manage 2007; 48: 3102–3119.

16.

Schwarz

Djilali

. 3D modeling of catalyst layers in PEM fuel cells. J Electrochem Soc 2007; 154: B1167–B1178.

17.

Shimpalee

Spuckler

Van Zee

. Prediction of transient response for a 25-cm(2) PEM fuel cell. J Power Sources 2007; 167: 130–138.

18.

Jang

. Three-dimensional numerical study on cell performance and transport phenomena of PEM fuel cells with conventional flow fields. Int J Hydrogen Energy 2008; 33: 156–164.

19.

Ren

. Transport mechanisms and performance simulation of a PEM fuel cell. Int J Energy Res 2008; 32: 514–530.

20.

Wang

. Modeling of two-phase transport in the diffusion media of polymer electrolyte fuel cells. J Power Sources 2008; 185: 261–271.

21.

Wang

Basu

Wang

. Modeling two-phase flow in PEM fuel cell channels. J Power Sources 2008; 179: 603–617.

22.

Berning

Odgaard

Kaer

. A computational analysis of multiphase flow through PEMFC cathode porous media using the multifluid approach. J Electrochem Soc 2009; 156: B1301–B1311.

23.

Schwarz

Beale

. Calculations of transport phenomena and reaction distribution in a polymer electrolyte membrane fuel cell. Int J Heat Mass Transfer 2009; 52: 4074–4081.

24.

Schwarz

Djilali

. Three-dimensional modelling of catalyst layers in PEM fuel cells: effects of non-uniform catalyst loading. Int J Energy Res 2009; 33: 631–644.

25.

Berning

Odgaard

Kaer

. A study of multi-phase flow through the cathode side of an interdigitated flow field using a multi-fluid model. J Power Sources 2010; 195: 4842–4852.

26.

Yuan

. Model prediction of effects of operating parameters on proton exchange membrane fuel cell performance. Renew Energy 2010; 35: 656–666.

27.

Berg

. Steady and unsteady 3D non-isothermal modeling of PEM fuel cells with the effect of non-equilibrium phase transfer. Appl Energy 2010; 87: 2778–2784.

28.

Cordiner

Lanzani

Mulone

. 3D effects of water-saturation distribution on polymeric electrolyte fuel cell (PEFC) performance. Int J Hydrogen Energy 2011; 36: 10366–10375.

29.

Fink

Fouquet

. Three-dimensional simulation of polymer electrolyte membrane fuel cells with experimental validation. Electrochim Acta 2011; 56: 10820–10831.

30.

Kang

. Liquid water flooding in a proton exchange membrane fuel cell cathode with an interdigitated design. Int J Energy Res 2011; 35: 1292–1311.

31.

Chiu

. A three-dimensional modeling of transport phenomena of proton exchange membrane fuel cells with various flow fields. Appl Energy 2012; 96: 359–370.

32.

Sun

. Numerical study of a 3D two-phase PEM fuel cell model via a novel automated finite element/finite volume program generator. Commun Comput Phys 2012; 11: 65–98.

33.

Cao

T-F

. Numerical investigation of the coupled water and thermal management in PEM fuel cell. Appl Energy 2013; 112: 1115–1125.

34.

Hossain

Islam

Pollard

. Investigation of species transport in a gas diffusion layer of a polymer electrolyte membrane fuel cell through two-phase modelling. Renew Energy 2013; 51: 404–418.

35.

Jiang

Wang

. Numerical modeling of liquid water motion in a polymer electrolyte fuel cell. Int J Hydrogen Energy 2014; 39: 942–950.

36.

Mancusi

. Numerical study of two-phase flow patterns in the gas channel of PEM fuel cells with tapered flow field design. Int J Hydrogen Energy 2014; 39: 2261–2273.

37.

Choopanya

Yang

. A CFD investigation of effects of flow-field geometry on transient performance of an automotive polymer electrolyte membrane fuel cell. Comput Therm Sci 2015; 7: 93–104.

38.

Ferreira

. Numerical simulations of two-phase flow in an anode gas channel of a proton exchange membrane fuel cell. Energy 2015; 82: 619–628.

39.

Kim

. Numerical simulation of water droplet dynamics in a right angle gas channel of a polymer electrolyte membrane fuel cell. Int J Hydrogen Energy 2015; 40: 8368–8383.

40.

Khazaee

Sabadbafan

. Effect of humidity content and direction of the flow of reactant gases on water management in the 4-serpentine and 1-serpentine flow channel in a PEM (proton exchange membrane) fuel cell. Energy 2016; 101: 252–265.

41.

Berning

Djilali

. Three-dimensional computational analysis of transport phenomena in a PEM fuel cell. J Power Sources 2002; 106: 284–294.

42.

Ticianelli

Beery

Srinivasan

. Dependence of performance of solid polymer electrolyte fuel-cells with low platinum loading on morphological-characteristics of the electrodes. J Appl Electrochem 1991; 21: 597–605.

43.

Uchida

. Investigation of the microstructure in the catalyst layer and effects of both perfluorosulfonate ionomer and PTFE-Loaded carbon on the catalyst layer of polymer electrolyte fuel cells. J Electrochem Soc 1995; 142: 4143–4149.

44.

Siegel

. Single domain PEMFC model based on agglomerate catalyst geometry. J Power Sources 2003; 115: 81–89.

45.

Rao

. A two-dimensional steady state model including the effect of liquid water for a PEM fuel cell cathode. J Power Sources 2007; 173: 375–393.

46.

Xing

. A two-phase flow and non-isothermal agglomerate model for a proton exchange membrane (PEM) fuel cell. Energy 2014; 73: 618–634.

47.

Wang

Cheng

. Multiphase flow and heat transfer in porous media. Adv Heat Transfer 1997; 30: 93–182.

48.

Baschuk

. Modelling of polymer electrolyte membrane fuel cells with variable degrees of water flooding. J Power Sources 2000; 86: 181–196.

49.

Cordiner

Mulone

Romanelli

. Thermal-fluid-dynamic simulation of a proton exchange membrane fuel cell using a hierarchical 3D-1D. J Fuel Cell Sci Technol 2007; 4: 317–327.

50.

Cordiner

. Polymer electrolyte fuel cell design based on three-dimensional computational fluid dynamics modeling. J Fuel Cell Sci Technol 2009; 6: 17–26.

51.

Dutta

Shimpalee

Van Zee

. Three-dimensional numerical simulation of straight channel PEM fuel cells. J Appl Electrochem 2000; 30: 135–146.

52.

Siegel

. A two-dimensional computational model of a PEMFC with liquid water transport. J Power Sources 2004; 128: 173–184.

53.

Shimpalee

. Predicting the transient response of a serpentine flow-field PEMFC I. Excess to normal fuel and air. J Power Sources 2006; 156: 355–368.

54.

Shimpalee

. Predicting the transient response of a serpentine flow-field PEMFC II: normal to minimal fuel and AIR. J Power Sources 2006; 156: 369–374.

55.

Sunden

Faghri

. Transport phenomena in fuel cells, International series on developments in heat transfer. Southampton: WIT Press, 2005.

56.

Tryggvason

Scardovelli

Zaleski

. Direct numerical simulations of gas-liquid multiphase flows, Cambridge University Press, 2011Cambridge.

57.

AVL, www.avl.com (2016, accessed 1 April 2016).

58.

ESI. CFD-ACE+, www.esi-cfd.com (2016, accessed 1 April 2016).

59.

ANSYS. ANSYS CFX, www.ansys.com (2016, accessed 1 April 2016).

60.

ANSYS. ANSYS CFX-Solver Theory Guide, www.ansys.com (2016, accessed 1 April 2016).

61.

COMSOL. COMSOL Multiphysics, www.comsol.com (2016, accessed 1 April 2016).

62.

ANSYS. ANSYS Fluent, www.ansys.com (2016, accessed 1 April 2016).

63.

ANSYS. Fluent Fuel Cells Module Manual, www.ansys.com (2016, accessed 1 April 2016).

64.

ANSYS. ANSYS Fluent Theory Guide, www.ansys.com (2016, accessed 1 April 2016).

65.

CD-adapco. STAR-CD, www.cd-adapco.com (2016, accessed 1 April 2016).

66.

CD-adapco. STAR-CCM+, www.cd-adapco.com (2016, accessed 1 April 2016).

67.

The OpenFOAM Foundation. OpenFOAM, www.openfoam.org (2016, accessed 1 April 2016).