1. Introduction
Hydrogen will play an important role as an energy carrier in the future. Renewable electrical energy is used to produce it, carbon dioxide-free, by electrolysis from water. Fuel cells convert the hydrogen back into electricity efficiently, and here, polymer electrolyte membrane fuel cells are showing promising performance [
1,
2]. Not only in stationary applications but also in mobility applications such as buses or truck propulsion systems, this design has been and is being successfully deployed [
3,
4]. The main advantages of PEMFCs are high efficiency, low operating noise, high startup availability, and high dynamic range [
5,
6]. The multiphase transport processes in membrane electrode assemblies (MEAs) and the durability and, hence, degradation of PEMFCs in cyclic operation are still areas of research. Transient 3D-CFD simulations with high temporal and spatial resolution and their ability to provide deep insights into internal processes and states accelerate development and reduce the need for prototypes [
7,
8,
9].
The development of simulation models of fuel cells started with the pioneering work of Springer [
10] and Bernardi and Verbrugge [
11]. During the past years, many models with different dimensionality and physical depth were presented [
8,
12,
13,
14,
15,
16,
17,
18]. With the increasing computational power, more sophisticated models such as unsteady, non-isothermal, two-phase 3D-CFD models were developed. Many authors reported the use of 3D-CFD PEMFC models with different commercial CFD codes, such as OpenFoam [
19], Ansys Fluent/CFX [
20,
21,
22,
23], StarCCM+ [
24,
25], and AVL Fire™ [
26].
One of the first 3D multiphase and multicomponent models of PEMFCs was developed by Berning and Djilali [
23] and implemented into Ansys CFX (former AEA Technology). Catalyst layers were not fully modeled but were represented by thin interfaces. Their main focus was on the formation of liquid water inside the gas channels and GDLs. Fink et al. [
27] presented an unsteady, multiphase, multicomponent and non-isothermal 3D-CFD model which is also used in this work to carry out the CFD simulations.
Some of the authors focused on single-channel models [
28,
29,
30,
31,
32], but with the increase of computational power, discretizations of the entire geometry [
19,
21,
33,
34,
35,
36,
37,
38,
39] or modeling degradation [
27,
40,
41] have become more common.
Iranzo et al. [
21] used the PEMFC model implemented in Ansys Fluent to give a guideline on how to cope with the parameterization challenge of such models. They had experimental data from one polarization curve and tried to fit the simulation results to this data by adjusting the reference exchange current density (RECD) of the cathode catalyst layer (CCL). The remaining parameters were taken from data sheets, literature review, and measurements. The simulation results and measured data showed good agreement with the polarization curve from experimental data. In [
35], the simulated liquid water distribution with the CFD model was validated with neutron imaging.
Macedo-Valencia et al. [
36] proposed an extension of the Ansys Fluent model by modeling the catalyst layers and membrane with computational cells as well. They simulated a short stack with five cells and an active area of
cm
. Their model was validated against available experimental data of a polarization curve but showed quite high deviations. The parameters of the 3D-CFD model were not perfectly fitted.
With increasing computational power, articles featuring transient 3D-CFD simulations of PEMFCs have been published in recent years. Iranzo and Boillat [
42] simulated the unsteady gas transport in the cathode side of a PEMFC during an AC impedance test. Khajeh-Hosseini-Dalasm et al. [
43] developed a 3D transient two-phase isothermal model to investigate the time variation of liquid water formation and the gas phase transport under startup conditions. Similarly, Peng et al. [
31] developed a transient heat and mass transfer 3D-CFD model to investigate a PEMFC with a dead-ended anode channel. They conducted simulations with and without purging of hydrogen at the anode outlet to show the differences in fuel cell performance. Variations in pressure and relative humidities of fuel and oxidant were investigated. Yan et al. [
38] conducted an unsteady 3D-CFD investigation of a serpentine flow field PEMFC to evaluate the effects of serpentine flow field design on the transient characteristics.
Verma and Pitchumani [
44] created a transient single-channel model of a PEMFC to study water transport, especially the influence of membrane properties on the behavior during dynamic load changes. They found that current density steps lead to drying of the anode due to increased electroosmotic drag. In [
45], the same model was used to conduct a model-based parametric study on the effects of operating parameters during a step change in load. Low inlet gas humidity significantly affects the transient behavior of the fuel cell until a steady state is reached. In [
46], they investigated the behavior of the single-channel PEMFC under transient changes in voltage, pressure, and stoichiometry in relation to the power demand of automotive cycles. They also simulated 180 s of the US06 drive cycle and optimized cell potential, inlet pressures, and mass flows to deliver the required performance.
Besides the numerical effort of fully resolved CFD analysis, the question of the parameterization of the transport and Butler–Volmer equations arises. The scatter of experimental data for material parameters of PEMFC models was highlighted by Vetter and Schumacher [
47]. Even with accessible data from data sheets of the materials, there is still an uncertainty due to compression of the parts, contact resistances, or its dependence on, e.g., water content and temperature [
48]. Moreover, new materials such as ultra-thin reinforced short-side-chain membranes may lack reliable experimental material parameters [
49].
Iranzo et al. [
21] proposed to use the reference exchange current density as a calibration parameter to fit the simulated data to experimental data. A similar approach was made by Fink et al. [
40]. Recently, many authors have developed parameterization strategies for models with reduced dimensionality [
50,
51,
52,
53,
54]. However, to the authors’ knowledge, there is no work that addresses the efficient parameterization of 3D-CFD models that predict a wider range of steady-state and transient operating conditions.
This work is intended to be a guide for efficient parameterization of an unsteady 3D-CFD model for industrial, low temperature PEMFCs with large active area. Using the same approach as in [
21,
40], the authors of this work had great difficulty in accurately predicting other available experimental data with the fitted RECD. Based on the previous work of the authors [
55], numerical optimization methods with the use of numerically reduced models were used to parameterize a 3D-CFD model. For parameterization and validation, measurement results of three load steps and three polarization curves were available, which were measured at the institute’s stack test stand. These data had to be reduced from the stack level to the cell level in order to obtain suitable boundary conditions for the single fuel cell simulations. Therefore, the presented results represented an average large active area fuel cell of a commercial 30 kW PEMFC stack. The experimental data included species mass flows at the inlets, pressures at the outlets, cooling temperatures, and gas humidities at the inlets. The high number of necessary model evaluations of numerical optimization methods were enabled by the use of numerically reduced models, but with very similar physical modeling approaches. With the results from the parameterization process, the transition to CFD was carried out with a single-channel model of the fuel cell. This model provided acceptable accuracy and low simulation times and it enabled the further refinement of parameters for the CFD models. The results show excellent agreement with the measured data, not only with the quasi-3D model but also with the 3D-CFD models. The obtained parameter set was able to predict measured load steps as well as measured polarization curves with minimum deviation.
Despite the reduction from stack to cell level, an alternative meshing strategy is presented to enable the transient 3D-CFD simulation of large active area fuel cells with reasonable computational effort. This approach uses porosity sections instead of the original flow field and steers the flow into the desired direction with anisotropic permeabilites. The number of computational cells could be reduced by a factor of 87% compared to a computational mesh with fully resolved gas channels.
Using this innovative meshing strategy and a suitable parameterization, the transient simulation of a load change was performed and the results are presented. The objective was to show the internal processes and states during transient operation, in particular the effects of liquid water formation and transport on the dissolved water content and electric current uniformity in the cell. Finally, an outlook on further possible investigations on the parameterization of transient CFD models is given.
2. Materials and Methods
2.1. Governing Equations of the 3D-CFD-Model
The 3D-CFD model for low-temperature PEMFCs used in this work is part of the AVL Fire™ M software package. The governing equations of the fuel cell model and the associated source terms and interface conditions can be found in Fink et al. [
27]. The model accounts for laminar two-phase flow based on an Euler multi-phase approach with multi-component species diffusion and heat transport in flow channels and porous layers. Electron and ion transport in solid phases, electrochemical reactions in catalyst layers and dissolved water, and species transport in ionomer material are considered. The electrochemical reactions are modeled with the well-known Butler–Volmer equation. Parasitic reactions at opposite electrode due to gas crossover in the membrane are considered. Gas diffusion layers (GDLs) and catalyst layers (CLs) are defined as porous material and are modeled using volume-averaging methods such as the Darcy equation, including capillary effects. Because the mass transport limitations at the cathode were not relevant, the more sophisticated agglomerate model was not used. The macrohomogeneous approach was activated, which further accelerated the simulation and to achieve better compatibility with the quasi-3D model used for the parameterization.
2.2. Geometry of the Investigated Fuel Cell
The investigated fuel cell was part of a 30 PEMFC stack. The bipolar plates were made out of graphite and had imprinted gas channels. Liquid cooling was used to remove excess heat, and the cooling channels were located inside the bipolar plates. GDLs were made from carbon paper with hydrophobic treatment and microporous layer (MPL). The total thickness of the MEA was 450 m in compressed state with 210 for each GDL, 10 for the CCL, 5 for the anode catalyst layer (ACL), and 15 for the membrane. Hydrogen and air were in counter-flow setup. The y-direction was the through-plane direction. The cooling inlet was on the side of the air inlet and hydrogen outlet.
2.3. Computational Meshes
The generation of suitable finite-volume meshes for 3D-CFD simulation (preprocessing) is one of the most important tasks in the course of setting up such simulations. The computational time, accuracy, and required computational resources are strongly influenced by the size and quality of the resulting computational meshes. The unique geometry of fuel cells, thin in the main current flow direction but up to two orders of magnitude larger in the other two dimensions, presents a challenge for proper discretization. In particular, the small dimensions of the gas channels and the minimum requirements for capturing all flow phenomena can drive the number of computational cells into regions as large as 100 M. The presented computational meshes were created with the AVL FAME™ M mesher, which is part of the software package used.
Three different meshing approaches are presented in the following sections. The first approach is the most simplified one, as it represents one gas channel of the fuel cell. In this approach, the behavior of the FC is averaged on a gas channel. As a result, the effects of non-symmetry in the in-plane direction (x-z plane) of the MEA are not captured. However, the resulting small number of computational cells keeps the simulation time low.
The most detailed approach is to discretize the whole fuel cell geometry as is. Local effects and processes can be captured with high detail, e.g., the influence of the land area between the channels on the species concentration is visible. The main drawback is the large number of computational cells for fuel cells with large active area.
A reduction of computational cells could be achieved by an innovative approach, the homogenized channels mesh. The imprinted flow field of the bipolar plate was replaced with porous sections and these sections steered the flow, with the help of anisotropic permeabilities, into the direction of the original flow field. This approach is not as detailed as the discretization of the entire FC, e.g., the effects of land areas between channels cannot be captured due to the averaging approach. However, as the results later show, the local distributions are predicted very similarly. The much smaller number of computational cells allows simulation of longer load profiles with moderate computational effort.
The MEA was extruded with 5 layers of each GDL, 6 layers for each CL, and 10 layers for the membrane (five at the anode and five at the cathode side). The MPL was not modeled in order to keep the number of unknown material parameters low.
2.3.1. Single-Channel Model (SCM)
The channel model mesh represents a simplification of the geometry by reducing the original fuel cell geometry to a channel. Compared to the later-presented discretization approaches, this leads to a much smaller cell count and thus faster computation times. Furthermore, the channel model represents the same part of the fuel cell as the quasi-3D model used for parameterization. Due to the serpentine flow field, the channel model is longer than the actual width of the fuel cell. Hexahedral cells were used to discretize the channel geometry. Symmetry conditions were applied to further reduce the number of computational cells (see
Figure 1). The terminal boundary faces were also used to apply cooling boundary conditions (BCs). The in/outlet sections of the gas channels were extruded to develop the flow before entering the fuel cell.
Table 1 shows the resulting number of computational cells.
2.3.2. Fully Resolved Channels Model (FRCM)
For the fully resolved mesh, anode and cathode bipolar plate were meshed independently with polyhedral cells. Since anode and cathode flow field design were different, no symmetry could be used. Inside the gas channels, one boundary layer was created to account for the wall treatment. The minimum cell size was chosen to have at least three computational cells (excluding the boundary layer cells) at the height of the gas channels to achieve a compromise between accuracy and computational time. GDL, CL, and the membrane were extruded starting from the created polyhedral mesh. The membrane is divided in the middle into an anode and a cathode part, connected with an arbitrary connection (non-conformal mesh interface) due to non-symmetry of the cell geometry.
The resulting number of generated computational polyhedral cells is shown in
Table 2. In
Figure 2a, the created mesh with the labeling of the parts is shown.
Figure 2b shows the meshed geometry of the fuel cell with highlighted boundary faces. Since the gas channels were very small compared to the outline of the fuel cell, this resulted in a large number of computational cells. Therefore, high-performance computers with a sufficient number of cores and RAM were required to perform the 3D-CFD simulations in a reasonable amount of time.
2.3.3. Homogenized Channels Model (HCM)
The computational mesh with fully resolved channels results in a high number of cells, which is not suitable for transient simulations due to high demands of computational power and time. The goal is to significantly reduce the number of computational cells, but maintain the general properties of the flow field. Therefore, the imprinted gas channels were modeled with porous media sections that steer the flow into the original flow direction of the gas channels via anisotropic permeabilities, e.g., if the main direction of the flow field is along the
x-axis, the permeability in
x-direction is orders of magnitude lower compared to the
y- and
z-direction. The goal of this approach is to maintain the pressure drop and local distributions of the other variables, as in the FRCM. The application of the porous media approach requires a hexahedral-structured computational mesh of the flow fields.
Figure 3 shows the basic principle. The dark sections inside the hexahedral mesh belong to the bipolar plate (=solid region) and were necessary to ensure numerical stability and avoid artificial overshoots of the velocity. Most of the calculated quantities were not visibly affected by the 90° change in direction from one flow steering section to the other. Only the velocity field and liquid water volume fraction show unexpected behavior in the corners of the turns where the skewness and growth ratios of the computational cells are high, but that does not have any big influence on the overall solution. Isotropic permeabilities are defined in the input region (green) to ensure that the flow is not forced to flow in a particular direction.
The hexahedral-structured mesh is embedded into the mesh with polyhedral cells of the bipolar plate and connected with non-conformal mesh interfaces. Based on the isometric views of the created meshes, no difference between the HCM and the FRCM can be seen (see
Figure 2b).
The number of generated computational cells of the homogenized channel mesh is shown in
Table 3. The homogenized channel approach reduced the mesh count by a factor of eight compared to the FRCM (see
Table 4). The porosity of the gas channels was defined as the ratio between the area where the gas channels have contact with the GDL and the active area. The permeability
K in the main flow direction of the respective flow section was calculated based on Darcy’s law [
56] in 1D:
where
v is the average velocity in the channel,
is the dynamic viscosity,
is the length of the flow steering section, and
is the pressure drop of the considered flow section. The values for the velocity
v and pressure drop
over the length
were extracted from a simulation with the fully resolved channels model. The design of the analyzed flow field required the definition of five flow steering sections: (1) Flow All Directions, (2) Flow X IO, (3) Flow Z Main, (4) Flow Z2, and (5) Flow X Main (see
Figure 4).
In the through-plane direction (y-direction), the permeability was set to a lower value to avoid too-high species concentration gradients over the depth of the channels.
Table 3.
Number of computational cells of homogenized channel approach.
Table 3.
Number of computational cells of homogenized channel approach.
| Cathode | Anode |
---|
Bipolar Plate | 1,158,408 | 1,301,646 |
Channels | 474,327 | 286,403 |
GDL | 377,240 | 377,240 |
CL | 452,688 | 452,688 |
Membrane | 754,480 |
Sum | 5,632,120 |
A comparison of the resulting cell count of the created computational meshes can be found in
Table 4.
Table 4.
Comparison of total number of computational cells of the resulting models.
Table 4.
Comparison of total number of computational cells of the resulting models.
| Number of Cells | Reduction |
---|
Fully Resolved Channels Model
(FRCM) | 43,367,985 | |
Homogenized Channels Model (HCM) | 5,632,120 | −87% to FRCM |
Single-Channel Model (SCM) | 680,880 | −88% to HCM |
2.4. Experimental Data
The used experimental data are described in detail in Haslinger et al. [
55]. Three polarization curve measurements at different anode inlet pressures and cooling temperatures, as well as three load steps with approximately 120
duration, were available for simulation and validation.
Figure 5 shows the measured curves.
2.5. Numerical Settings
The numerical settings can influence a CFD simulation significantly and have to be customized for every simulation. In this paper, results from steady-state and transient simulations are presented. The maximum steady-state iterations for the FRCM and HCM were set to 7500 to obtain a converged solution. The initialization of the transient simulations was performed with quasi-steady-state time steps to have a converged solution as initial condition for the transient time stepping. The time step sizes were set according to the gradients of voltage, current density, and air mass flow and varied from 50 ms (high gradients) to 200 ms (low gradients). The iterations per time step ranged from 20 for SCM to 50 for HCM. With that approach, sufficient low residuals were achieved.
2.6. Boundary Conditions
In
Table 5, the boundary conditions are listed. The electric potential at the anode terminal surfaces was 0
and
at the cathode side for voltage-driven simulations (load steps). For the simulation of the polarization curves, the current density was prescribed at the cathode terminal faces. Excess heat was removed by heat convection, and a heat transfer coefficient and average cooling temperature were specified at the cooling faces. The remaining boundary faces were of the wall type, with no heat flux and a current density of zero. In the SCM, the electrical connection and cooling surfaces coincide. In addition, the SCM model had symmetric boundary conditions to reduce the number of computational cells.
2.7. Model Parameterization
A transient, single-phase, non-isothermal, quasi-3D model of a PEMFC based on the previous work of the authors [
55] was used to conduct the numerical optimization of parameters. The identification of the parameters with the greatest influence on performance was based on the results of Karpenko-Jereb et al. [
28]. They used the same CFD code as in this work to perform a theoretical study of the influence of parameters on the performance of a PEMFC. Their result is consistent with other sensitivity analyses of fuel cell models found in the literature [
18,
57].
It turned out that some of the parameters obtained from the numerical optimization process need further adjustments when used with the CFD models to predict the measured current density and potential accurately. The reasons for the discrepancies were, first, a different membrane water model, and, second, that the quasi-3D model did not account for the ohmic losses in the CLs [
58]. Low simulation times of the SCM allowed fast parameter variations with that model due to the reduced cell count compared to the FRCM and HCM. The following methodology was developed to achieve the best agreement with the measured data:
Quasi-3D Model:
Perform numerical optimization of parameters with a high number of function evaluations and measurement data as reference.
Transfer to CFD:
The reference ionic conductivity from the numerical optimization with the quasi-3D model had to be multiplied by 3. This was necessary in order to account for the ionic losses in membrane, ACL, and CCL because the used quasi-3D model did not account for ionic losses in the CLs. The factor of 3 was derived from the thicknesses of the parts and the assumption that the proton concentration in- and decreases with a constant gradient in the CLs. When ultra-thin membranes are used, the ion losses of the CLs may exceed the ion losses in the membrane. Therefore, modeling of the ion losses in the CLs is required to account for the actual ohmic losses. The SCM, as the simplest and fastest of the CFD models, was used to obtain a first glance at how the derived parameters behave with the CFD model. It turned out that some minor adjustments to the reference exchange current density and ionic conductivity were necessary to obtain a better agreement of the results of the CFD models with the measured data.
HCM/FRCM:
Use the updated parameter set to simulate more sophisticated models such as the FRCM or HCM.
For the optimization, only load step L3 was used as a reference; the other measured curves were used for validation purposes of the results.
4. Conclusions
The parameterization and the vast amount of resulting computational cells of 3D-CFD models of large active area FCs are challenging. In this work, both issues were addressed. The optimal parameter set was obtained using a quasi-3D model and numerical optimization methods with experimental data as a reference. A set of eight parameters (design variables), characterizing different parts of the polarization curve, were selected. Experimental data for calibration and validation from a stack test bench were available. Three measured polarizaton curves with different anode inlet pressures and cooling temperatures, as well as three load steps in different current regions, were recorded at the test bench. The load step with the biggest change in load was used as a reference for numerical optimization of parameters.
For the 3D-CFD simulations, three different meshing approaches with different levels of detail were presented: the SCM, FRCM, and HCM. As a first step, the SCM was parameterized with the obtained parameter set by the numerical optimization procedure. With minor adjustments to the reference exchange current density and the reference ionic conductivity, very good agreement between simulated and measured load steps and polarization curves was obtained. The use of only one load step that covers the full range of possible fuel cell current is enough to obtain a set of parameters that are able to simulate not only dynamic load changes but also steady-state polarization curves with high accuracy.
The SCM has only minor capabilities of representing the actual distribution of quantities over the active area. Thus, more detailed 3D-CFD models were created. The FRCM resulted in about 43 M computational cells, which might be used for steady-state simulations but not applicable for unsteady simulations of up to several minutes of physical time. To enable transient 3D-CFD simulations of the entire geometry of the fuel cell and simulate more than a few seconds in a reasonable computation time, an alternative meshing approach, the homogenized channel model, had to be used. With that approach, it was possible to reduce the cell count by about 87%, enabling the unsteady simulation of a load step with the high temporal and spatial distribution. The final parameter set obtained and improved with the SCM was used for parameterization of the HCM.
The steady-state results between the fully resolved and homogenized channel approach revealed only minor differences in pressure drop and distributions of solution quantities. One of the three available load steps (L3) from experimental data was simulated with the HCM, and the simulation results showed an excellent agreement of voltage and current with the measured data. Further insights into the processes and especially the connection of water management and current density during a steep load step of a PEMFC were made visible and discussed. Possible future work may include the reduction of design variables for the numerical optimization process. In addition, the authors experienced uniqueness problems with the obtained parameter set, as discussed in [
60]. Therefore, more detailed measurements need to be performed on the materials used in order to obtain a parameterization of the fuel cell model that reflects the actual material parameters as well as possible.