**Optimization of Membrane Electrode Assembly of PEM Fuel Cell by Response Surface Method**

**Rohit K. S. S. Vuppala 1,2, Benitta A. Chaedir <sup>1</sup> , Lishuai Jiang <sup>3</sup> , Lianjun Chen <sup>3</sup> , Muhammad Aziz <sup>4</sup> and Agus P. Sasmito 1,\***


Received: 10 July 2019; Accepted: 21 August 2019; Published: 26 August 2019

**Abstract:** The membrane electrode assembly (MEA) plays an important role in the proton exchange membrane fuel cell (PEMFC) performance. Typically, the structure comprises of a polymer electrolyte membrane sandwiched by agglomerate catalyst layers at the anode and cathode. Optimization of various parameters in the design of MEA is, thus, essential for reducing cost and material usage, while improving cell performance. In this paper, optimization of MEA is performed using a validated two-phase PEMFC numerical model. Key MEA parameters affecting the performance of a single PEMFC are determined from sensitivity analysis and are optimized using the response surface method (RSM). The optimization is carried out at two different operating voltages. The results show that membrane thickness and membrane protonic conductivity coefficient are the most significant parameters influencing cell performance. Notably, at higher voltage (0.8 V per cell), the current density can be improved by up to 40% while, at a lower voltage (0.6 V per cell), the current density may be doubled. The results presented can be of importance for fuel cell engineers to improve the stack performance and expedite the commercialization.

**Keywords:** PEM fuel cell; membrane electrode assembly (MEA); response surface method; computational fuel cell dynamics

#### **1. Introduction**

As a clean energy device, a proton exchange membrane (PEM) fuel cell is a promising power-generating technology that has received increasing attention over the last decade. Fuel cell is an electrochemical device that converts chemical energy into electrical energy. Due to its high energy-conversion efficiency and zero-emission potential, the fuel cell is considered as the best power-generating device, especially in transportation applications. Among the different kinds of fuel cells, PEM fuel cell (PEMFC) offers desirable features, such as low operational temperature and high-power density, which makes it the most promising alternative technology for power production.

In order for the fuel cell technology to be competitive with conventional power systems, some challenges associated with it, including reliability, longevity, and cost, must be overcome. A better understanding of the system is required to achieve the ideal price-performance balance. Studies have been carried out to characterize the behavior of the PEMFC system as affected by different parameters. Wang et al. [1] conducted parametric experiments to study the effect of various operating parameters on the performance of a single PEMFC and used the results to validate the 3D model they developed. The parameters studied were pressure, fuel cell temperature, and anode and cathode relative humidity. It was observed that, generally, increasing the fuel cell temperature and pressure increases its performance, except when the temperature is higher than the gas stream humidification temperatures, especially at a low current density region. Cathode humidification temperature was found to have no significant impact on fuel cell performance, while increasing anode humidification temperature increases performance at the low current density region. These results are in accordance with the results obtained by Ferng et al. [2] and Amirinejad et al. [3] who concluded that operation at higher pressure and elevated temperature can improve the electrode kinetic performance and increase the ionic conductivity in membrane and electrodes, which results in high power density in the PEMFC system. Santarelli and Torchio [4] experimentally analyzed the performance of a PEMFC by varying cell temperature, anode and cathode flow temperatures in both saturation and dry conditions, and reactant pressure. The results showed that, in addition to cell temperature, an increase in reactant saturation temperature also leads to a better performance and the best improvements due to a pressure increase are observed when both anode and cathode are humidified. Yan et al. [5] extended the study of the effects of operating conditions to PEMFC with interdigitated flow field. Nafion-based PEM fuel cell performance analysis with various reactant humidification levels, which varied from no external humidification to a fully saturated anode and cathode, was carried out by Williams et al. [6]. Klika et al. [7] have developed a thermodynamically-consistent model based on polynomial functions to study the behavior of water sorption in Nafion membranes. A three-dimensional multiphase numerical model was developed by Fan et al. [8] to study the PEMFC performance at low external humidification. It was found that the dependency on external humidification of a PEMFC can be cut down at high current density, due to the water produced in the cathode catalyst layer that is sufficient to be employed to humidify both the cathode and anode polymer electrolyte.

In addition to operating conditions, there are other parameters affecting the fuel cell performance. Bayrakçeken et al. [9] found that membrane thickness, hot-pressing conditions of the gas diffusion layer (GDL), and the Teflon to carbon ratio in the GDL, are also significant parameters to provide good PEMFC performance. The study showed that thinner membrane thickness and higher Teflon:carbon ratio in the GDLs give better performances. Jiang et al. [10] implemented an effective "elementary effect" (EE) method based on Monte Carlo random experiments to analyze 22 uncertain parameters involved in their two-phase 1D analytical PEMFC model. Among all of the parameters, membrane thickness and volume fractions were found to be the most important factors influencing the cell performance. The effect of catalyst layer microstructure was recently investigated numerically by Carcadea et al. [11]. A CFD model was used to study the behavior of a PEMFC as a function of Pt loading, Pt particle radius, ionomer volume fraction, and carbon support, and to establish the optimum range of these parameters. It was observed that increasing the ionomer volume fraction in the catalyst layer (CL) leads to better performance due to the fact that the ionomer acts as a network for the mass and charge transport. Moreover, higher Pt loading and a lower particle radius are recommended to achieve better PEMFC performance. Lee et al. [12] investigated the performance improvement of a PEMFC as a function of gas diffusion layer porosity and impregnation of the Nafion solution.

The previously mentioned studies confirm the significance of various parameters on the operation of the PEMFC. It is, therefore, crucial to select the optimum values in order to achieve a high-performance fuel cell. Efforts have been made by researchers toward the optimization of critical parameters influencing the PEMFC operation using different approaches. Salva et al. [13] developed a one-dimensional analytical model and used it to obtain the operating conditions under which a single PEMFC provides the maximum power output for different current intensities. The optimization was carried out for every value of current intensity by solving the parametric table consisting of all possible combinations obtained from modifying the stoichiometry in the cathode and anode, relative humidity in the anode and cathode, and the operating temperature, while keeping the pressure constant. Wu et al. [14] employed a multi-resolution approach and the radial basis function (RBF) surrogate model for simulation and optimization of operating conditions for hydrogen polymer electrolyte fuel cells. Zervas et al. [15] performed a phosphoric acid fuel cell (PAFC) system optimization study based on meta-models that were derived by applying the linear regression and the RBF neural network methodology on the results produced by a CFD model. The optimization of different operating and design parameters on PEMFC using the Taguchi method was performed by Karthikeyan et al. [16], Solehati et al. [17], and Sasmito et al. [18]. Grujicic and Chittajallu [19] utilized a single-phase two-dimensional electrochemical model, coupled with a nonlinear constrained optimization algorithm, which was solved using sequential quadratic programming (SQP) to obtain the operational and geometric parameters for achieving the maximum electric current in a PEMFC. The parameters investigated include air inlet pressures and cathode thickness, cathode length for each shoulder segment of flow channel, and a fraction of cathode length associated with the flow channel. Similarly, Na and Gou [20] used SQP to optimize the fuel cell system efficiency and cost. Guo et al. [21] proposed an optimization algorithm that combines the teaching–learning based optimization (TLBO) with a differential evolution (DE) algorithm, known as the TLBO-DE method, to promote the efficiency of PEMFC. Behrou et al. [22] demonstrated the use of density-based topology optimization for the practical design of flow fields for PEMFCs, with goals to maximize both the output power and homogeneity of the current density distribution, as well as permit reduced costs and higher durability. The response surface methodology (RSM) has been employed by Kanani et al. [23] and Xuan et al. [24] to maximize the performance of a PEMFC system. Recently, a comprehensive evaluation of different optimization scenarios for a PEMFC is provided by Sohani et al. [25].

Optimization of controlling parameters at the fuel cell system level, like membrane thickness, size of cathode catalyst particle, and protonic conductivity coefficient of the membrane, however, have been very limited. This can be attributed to the fact that they cannot be changed during the cell utilization [4], which makes it tedious and uneconomical to perform these studies experimentally. Moreover, the complex structure of MEA, comprising of polymer electrolyte membrane sandwiched by agglomerate catalyst layers at both the anode and cathode, adds complexity for the design of experiments due to a multiscale nature of the system. In addition, none of the research studies had focused on the membrane electrode assembly coupled with agglomerate catalyst layer parameters. This paper aims to develop a numerical model to simulate a polymer electrolyte membrane fuel cell with detailed multiscale MEA with an agglomerate catalyst layer model, and determine the optimum values of the previously mentioned parameters that provide maximum current density for various voltage values in the ideal range of operation. The study focuses on sensitivity analysis of design parameters of the membrane electrode assembly, including membrane thickness, membrane equivalent weight, the radius of the cathode catalyst particle, catalyst ionomer resistance, cathode catalyst porosity, a membrane protonic conductivity coefficient, a cathode catalyst hydrophobic angle, ionomer tortuosity, and a cathode catalyst volume fraction. The parameters which have significant impact on the current density magnitude are considered. Once the model is validated, it is used to carry out parameter optimization for better cell performance. In the optimization, the response surface methodology is employed for meta-modelling. RSM is a collection of statistical and mathematical methods for optimizing and predicting responses with limited experimental data at various input factors, as well as performing sensitivity analysis [23]. It is extensively used in the industrial world, particularly in situations where the output is swayed by several input parameters. This method has been widely used in different fields and applications such as metals removal [26], chemical extraction [27,28], and the chemical and environmental engineering field [29–31]. As compared to other methods, RSM being a collection of mathematical and statistical techniques, gives a better understanding into the role of different parameters at play and generates a continuous model for visualizing the effect of parameters on the entire range as opposed to the average value of the response [32]. This study pioneers the sensitivity analysis of parameters of MEA coupled with agglomerate catalyst layers using the design of the experiment RSM method, along with the validated three-dimensional numerical model. The total numbers of simulations for combinations of 10 parameters would have been computationally

**2. Methodology**

very expensive and could not have been done using traditional parametric studies. The RSM method reduces the number of simulations significantly. *2.1. Mathematical Formulation*

#### **2. Methodology** The schematic figure of a typical PEMFC and its functional layers is illustrated in Figure 1. The system consists of a proton exchange membrane (m), sandwiched between two catalyst layers (cl),

#### *2.1. Mathematical Formulation* two gas diffusion layers (gdl), two porous metal foam flowfields (ff), and two terminal plates (tp).

The schematic figure of a typical PEMFC and its functional layers is illustrated in Figure 1. The system consists of a proton exchange membrane (m), sandwiched between two catalyst layers (cl), two gas diffusion layers (gdl), two porous metal foam flowfields (ff), and two terminal plates (tp). The main assumptions/approximations adopted in the model are: The main assumptions/approximations adopted in the model are: • Thermal equilibrium: Local thermal equilibrium between all the phases. • Membrane: The membrane model takes into account the water flux due to electro-osmatic drag and diffusion.


**Figure 1.** Schematic view of a single PEMFC: (**a**) components of PEMFC, (**b**) computational domains with boundaries: I—insulation/wall, II—anode inlet, III—cathode inlet, and IV—outlets. **Figure 1.**Schematic view of a single PEMFC: (**a**) components of PEMFC, (**b**) computational domains with boundaries: I—insulation/wall, II—anode inlet, III—cathode inlet, and IV—outlets.

**Table 1.** Geometrical and operating parameters for two PEMFC experimental data sets. **Case a Segmented Cell Case Single-Gas Diffusion Layer [37]** Physical parameters , 7.3 × 10<sup>−</sup><sup>13</sup> m<sup>2</sup> 6.1 × 10<sup>−</sup><sup>11</sup> m<sup>2</sup> The mathematical model is comprised of governing equations for the conservation of mass, momentum, species, energy, charge, and water transport in the membrane. The physical parameters, geometry, and operating conditions for two different PEMFC experimental data sets that are used later for validation purposes can be found in Tables 1 and 2.

 0.9 0.635 0.4 0.77 **Table 1.** Geometrical and operating parameters for two PEMFC experimental data sets.


3 × 10<sup>−</sup><sup>5</sup> m 5.1 × 10<sup>−</sup><sup>5</sup> m


**Table 1.** *Cont.*

**Table 2.** Additional parameters for all cases.


2.1.1. Governing Equations

Conservation of Mass [34]:

$$\frac{\partial \rho}{\partial t} + \nabla \cdot \left(\rho \begin{array}{c} \overrightarrow{v} \end{array} \right) = \mathcal{S}\_m \tag{1}$$

where *S<sup>m</sup>* is the mass source added from the continuous phase to the dispersed second phase and any other user-defined sources.

Conservation of Momentum [34]:

$$\frac{\partial}{\partial t}(\rho \stackrel{\scriptstyle}{\boldsymbol{\upsilon}} \overline{\boldsymbol{\upsilon}}) + \nabla \cdot \left(\rho \stackrel{\scriptstyle}{\boldsymbol{\upsilon}} \stackrel{\scriptstyle}{\boldsymbol{\upsilon}} \overline{\boldsymbol{\upsilon}}\right) = -\nabla p + \nabla \cdot \left(\overline{\overline{\boldsymbol{\tau}}}\right) + \rho \stackrel{\scriptstyle}{\boldsymbol{\varsigma}} + \mathcal{S}\_{\text{mom}} \tag{2}$$

where *<sup>p</sup>* denotes pressure, <sup>=</sup> τ is the stress tensor, ρ → *g* denotes the gravitational body forces, and *Smom* is the momentum source term for porous media, which includes the gas diffusion layer, catalyst layer, and membrane. The stress tensor <sup>=</sup> τ is given by the equation below [34].

$$\overline{\overline{\pi}} = \mu \left[ \left( \nabla \ \overrightarrow{\overline{v}} + \nabla \ \overrightarrow{\overline{v}}^{T} \right) - \frac{2}{3} \ \nabla \cdot \overrightarrow{v} \ I \right] \tag{3}$$

where µ is the molecular viscosity and *I* is the unit tensor. The second term in the right-hand side of the equation represents the effects of volume dilation.

Species Transport [34]:

$$\frac{\partial}{\partial t}(\rho \,\, Y\_i) + \,\,\nabla \cdot \left(\rho \,\, \overrightarrow{\boldsymbol{v}} \,\, \, Y\_i\right) = \,\, -\nabla \cdot \boldsymbol{D}\_{i\varepsilon f\bar{f}} \,\, \mathrm{V} \mathbf{Y}\_i + \mathcal{S}\_m \tag{4}$$

where *Y<sup>i</sup>* denotes the local mass fraction of species *i* and *Di*,*e f f* is the effective diffusivity of the species. Note that the total species mixture should conserve the total mass, and, thus, the source terms in the conservation of mass should be equal to the source terms in the conservation of species [48].

Electric Potential [34]:

$$\nabla \cdot (\sigma \,\, \nabla \psi) + \mathbf{S} = \mathbf{0} \tag{5}$$

where ψ is the electric potential, σ is the electric conductivity in a solid zone or ionic conductivity in a fluid zone, and *S* is a source term.

Conservation of Energy [34]:

$$\frac{\partial}{\partial t}(\rho \to) + \nabla \cdot \left(\rho \mathbb{C}\_p \vec{\upsilon} \, T\right) = \nabla \cdot \left(\mathbb{k}\_{eff} \, \nabla T\right) + \mathbb{S}\_h \tag{6}$$

where *C<sup>p</sup>* is the specific heat capacity and *ke f f* is the effective thermal conductivity. The first term on the right-hand side of the equation represents the energy transfer due to conduction.

Volumetric source terms (*Sm*) for H<sup>2</sup> and O<sup>2</sup> and the dissolved water content λ in the triple-phase boundaries (catalyst layers) due to electrochemical reactions are shown below [34].

$$S\_{H\_2} = -\frac{M\_{wH\_2}}{2F} R\_{an} < 0\tag{7}$$

$$S\_{O\_2} = -\frac{M\_{\rm wO\_2}}{2F} R\_{\rm cat} < 0\tag{8}$$

$$S\_{\lambda} = -\frac{M\_{wH\_2O}}{2F} \mathcal{R}\_{\text{cat}} > 0\tag{9}$$

where *MwH*<sup>2</sup> , *MwH*2*O*, and *MwH*<sup>2</sup> are the molecular mass of hydrogen, oxygen, and water, respectively, and *F* is the Faraday constant.

Mass transfer and water transport occurring in the PEMFC model is considered to be in two different phases, which are discussed below.

#### 1. Dissolved phase

The dissolved phase exists in the catalyst layers (ionomers) and in the membrane. The generation and transport of the dissolved water is described by the equation below [49].

$$\frac{\partial}{\partial t} \Big( \varepsilon\_i M\_{w, H\_2O} \frac{\rho\_i}{EW} \,\lambda \Big) + \nabla \cdot \left( \stackrel{\rightharpoonup}{\mathcal{I}}\_{\mathcal{U}} \frac{\eta\_{\mathcal{U}}}{F} M\_{\mathcal{W}} \right) = \nabla \cdot \Big( M\_w D\_w^i \,\nabla \,\lambda \Big) + S\_{\lambda} + S\_{\mathcal{g}d} + S\_{\mathcal{l}d} \tag{10}$$

where *<sup>i</sup>* denotes the porous media porosity, λ denotes the dissolved water content, η is the osmotic drag coefficient, and <sup>→</sup> ı *<sup>m</sup>* is the ionic current density, calculated as <sup>→</sup> ı *<sup>m</sup>* = −β*m*σ*mem* ∇ φ*mem*. In the right hand-side of the equation, *D<sup>i</sup> <sup>w</sup>* represents the diffusion coefficient of the water content, *S*<sup>λ</sup> denotes the water generation rate due to the cathode side reaction in the catalyst layer, *Sgd* is the rate of mass change between gas and dissolved phases, and *Sld* is the rate of mass change between the liquid and the dissolved phases. The mass change rates are shown in the equations below [50].

$$\mathcal{S}\_{\mathcal{S}^d} = \left(1 - s^{\theta}\right) \gamma\_{\mathcal{S}^d} M\_{w, H\_2O} \, \frac{\rho\_{\mathcal{i}}}{EW} \left(\lambda\_{eq} - \lambda\right) \tag{11}$$

$$\mathcal{S}\_{ld} = \left(\mathbf{s}^{\Theta}\right) \gamma\_{ld} M\_{w, H\_2O} \,\, \frac{\rho\_l}{EW} \left(\lambda\_{\text{eq}} - \lambda\right) \tag{12}$$

where ρ<sup>i</sup> is the dry ionomer or membrane density, *EW* is the equivalent weight of the membrane, *s* denotes the liquid saturation, λ*eq* denotes the equilibrium water content, θ is the exponential liquid coverage, and γ*gd*, γ*ld* are the gas and liquid mass exchange rate constants and are user-specified parameters. The equilibrium water content λ*eq* can be calculated using the equation below [50].

$$\lambda\_{eq} = 0.3 + 6a(1 - \tanh(a - 0.5)) + 0.69(\lambda\_{a=1} - 3.52)a^{0.5} \left(1 + \tanh\left(\frac{a - 0.89}{0.23}\right)\right) \tag{13}$$
 
$$+ \text{s.} (\lambda\_{s=1} - \lambda\_{a=1})$$

where *a* is the water activity, defined as:

$$a = \frac{p\_{\text{uv}}}{p\_{\text{sat}}} \tag{14}$$

where *pwv* is the water vapour partial pressure and *psat* is the saturation pressure. Both λ*s*=<sup>1</sup> and λ*a*=<sup>1</sup> in Equation (13) are user-specified parameters.

#### 2. Liquid Phase

Liquid is present in all the porous electrodes and gas channels. The driving force of the liquid water transport is the liquid pressure gradient ∇ *pl* , as shown in the liquid water transport equation below [50].

$$\frac{\partial}{\partial t}(\varepsilon\_l \rho\_l \, s) = \,^\cdot \nabla \cdot \left(\frac{\rho\_l \, \text{KK}\_l}{\mu\_l} \, \text{V}p\_l\right) + \,^\cdot \text{S}\_{gl} - \,^\cdot \text{S}\_{ld} \tag{15}$$

In Equation (15), ρ*<sup>l</sup>* is the liquid water density, µ*<sup>l</sup>* is the liquid dynamic viscosity, *K* is the absolute permeability, *K<sup>r</sup>* is the relative permeability, *p<sup>l</sup>* is liquid pressure, and *Sgl* is the rate of mass change between the gas and liquid phases. Replacing the liquid pressure with the sum of capillary pressure *p<sup>c</sup>* and gas pressure *p*, Equation (15) can be written as the equation below.

$$\frac{\partial}{\partial t}(\varepsilon\_{l}\rho\_{l}\,\mathrm{s}) = \,\,\,\nabla \cdot \left(\frac{\rho\_{l}\mathrm{K}\mathrm{K}\_{\mathrm{r}}}{\mu\_{l}}\,\,\,\nabla(p\_{c}+p)\right) + \,\,\,\mathcal{S}\_{gl} - \,\,\mathrm{S}\_{ld} \tag{16}$$

The mass transfer rate between the gas and liquid phases can be calculated based on the unidirectional diffusion theory [50,51].

$$\mathcal{S}\_{gl} \begin{cases} & \mathcal{Y}\_{gl} \in \text{sD}\_{\mathcal{S}l} \frac{M\_{\text{w},H\_2O}}{RT} p \ln \left( \frac{p - p\_{\text{sat}}}{p - p\_{\text{uv}}} \right) & p\_{\text{uv}} \le p\_{\text{sat}}\\ & \mathcal{Y}\_{\mathcal{S}l} \in (1 - s) D\_{\mathcal{S}l} \frac{M\_{\text{w},H\_2O}}{RT} p \ln \left( \frac{p - p\_{\text{sat}}}{p - p\_{\text{uv}}} \right) & p\_{\text{uv}} > p\_{\text{sat}} \end{cases} \tag{17}$$

where is porosity, <sup>γ</sup>*gl* is the geometric factor of the droplet size, and <sup>D</sup>*gl* <sup>h</sup> m2/s i , in the function of temperature [K] and pressure [Pa], takes the following form.

$$D\_{gl} \begin{cases} \ 0.365 \cdot 10^{-4} \left( \frac{T}{343} \right)^{2.334} \cdot \left( \frac{10^5}{p} \right), & \text{cathode} \\\ \ 1.79 \cdot 10^{-4} \left( \frac{T}{343} \right)^{2.334} \cdot \left( \frac{10^5}{p} \right), & \text{anode} \end{cases} \tag{18}$$

Equation (16) is solved in all the regions from the anode GDL-channel interface to the cathode GDL-channel interface. At both interface boundaries, the liquid water flux is considered to leave the gas diffusion layers into the gas channel only. The flux is assumed to be driven by the capillary pressure [50].

$$f\_{\rm liq} = \max[\Theta \epsilon \, sp\_{\mathcal{L}}, 0] \tag{19}$$

where Θ is the coefficient of liquid water removal.

Liquid saturation in the channels is calculated from the Leverett function.

$$p\_c \begin{cases} \begin{array}{c} \sigma \cos \theta\_c \ \sqrt{\frac{\xi}{K}} \text{ J } (1-s) \text{.} \\ \sigma \cos \theta\_c \ \sqrt{\frac{\xi}{K}} \text{ J s} \end{cases} & \begin{array}{c} \theta\_c < 90^{\circ} \\ \theta\_c > 90^{\circ} \end{array} \end{cases} \tag{20}$$

where

$$J(\mathbf{x}) = 1.417\mathbf{x} - 2.12\mathbf{x}^2 + 1.263\mathbf{x}^3 \tag{21}$$

Liquid water transport in the gas channels is determined to predict the pressure drop increase using the following correlation.

$$\frac{\partial}{\partial t}(\rho\_l \,\mathrm{s}) \, \, \, \, \nabla \cdot \Big(\boldsymbol{p}\_l \,\overrightarrow{\boldsymbol{v}}\_l \,\mathrm{s}\,\Big) = \,\, \, \nabla \cdot \Big(\boldsymbol{D}\_{\mathrm{liq}} \,\, \mathrm{Vs}\,\Big) \tag{22}$$

where *<sup>D</sup>liq* is the liquid water diffusion coefficient in the gas channel and the liquid velocity <sup>→</sup> *v l* is assumed to be a fraction of the gas velocity <sup>→</sup> *v <sup>g</sup>*, i.e., <sup>→</sup> *v <sup>l</sup>* = χ → *v <sup>g</sup>*.

Since the flow-channels (ff) in our model are porous in nature, user-defined function (UDF) is used to add the corresponding source term to X, Y, and Z momentum for the inertial losses.

$$\mathcal{S}\_{\bar{i}} = -\left(\frac{\mu}{\alpha} \,\, v\_{\bar{i}} + \,\, \mathcal{C}\_2 \,\, \frac{1}{2} \,\, \rho \,\, |v| \, v\_{\bar{i}}\right) \tag{23}$$

where *S<sup>i</sup>* denotes the source for the *ith* (x, y, or z) momentum equation, |*v*| denotes the magnitude of the velocity, α is the permeability, and *C*<sup>2</sup> is the inertial resistance factor.

In laminar flows through porous media, the pressure drop is typically proportional to the velocity and the constant *C*<sup>2</sup> can be considered zero. Ignoring the convective acceleration and diffusion, the porous media model is reduced to Darcy's law.

$$S\_i = -\frac{\mu}{\alpha} \,\, v\_i \tag{24}$$

The volumetric heat sources in various zones can be found in Table 3, where the variables *i<sup>s</sup>* and *i<sup>m</sup>* represent the magnitudes of the solid phase and membrane phase current density, respectively, and *L* (<0) is the latent heat due to water condensation.


**Table 3.** Energy source terms for the governing equations of conservation energy [50].

#### 2.1.2. Electrochemistry and Cathode Particle/Agglomerate Model

The driving force of these reactions is the surface overpotential, which is the difference between the phase potential of the solid and of the electrolyte or membrane. The phenomenon is accounted for in two equations: one for the electron transport in the catalyst layer, solid grids of porous media, and the current collector, and the other for the protonic conduction or transport of H<sup>+</sup> at the catalyst and the membrane [52,53].

$$\nabla \cdot \left( \sigma\_{\text{sol}} \, \, \nabla \, \phi\_{\text{sol}} \right) + R\_{\text{sol}} = 0 \tag{25}$$

$$\nabla \cdot \left( \beta\_m \sigma\_{mem} \,\, \nabla \, \phi\_{mem} \right) + R\_{mem} = 0 \tag{26}$$

where σ denotes the electrical conductivity in ohm-m−<sup>1</sup> , φ denotes the electric potential in volts, and *R* denotes the volumetric transfer current in A·m−<sup>3</sup> , which is also known as exchange current density, expressed as:

$$R\_{an} = \begin{pmatrix} \zeta\_{an} \ j\_{an} \ (T) \end{pmatrix} \begin{pmatrix} [A] \\ \hline [A]\_{ref} \end{pmatrix}^{\prime \mu a} \left( e^{a\_{\rm int}^{an} F \eta\_{\rm int} / RT} - e^{-a\_{\rm ext}^{an} F \eta\_{\rm int} / RT} \right) \tag{27}$$

$$R\_{\rm cat} = \left(\zeta\_{\rm cat} \, j\_{\rm cat}\left(T\right)\right) \left(\frac{[\mathbb{C}]}{[\mathbb{C}]\_{ref}}\right)^{\mathbb{V}\_{\rm cat}} \left(-e^{\alpha\_{\rm an}^{\rm{at}}F\eta\_{\rm cat}/RT} + e^{-\alpha\_{\rm cat}^{\rm{at}}F\eta\_{\rm cat}/RT}\right) \tag{28}$$

In the above equations, *j* (*T*) is the reference exchange current density per active surface area [A·m−<sup>2</sup> ], ζ is the specific active surface area [m−<sup>1</sup> ], [ ] and [ ]*re f* are the species local concentration and its reference value [kmol·m−<sup>3</sup> ], γ is the concentration dependence, α *an an* and α *an cat* are anode and cathode dimensionless transfer coefficients of the anode electrode, respectively, α *cat an* and α *cat cat* are the anode and cathode dimensionless transfer coefficients of cathode electrode, η*an* is the surface overpotential, F is the Faraday constant (9.65 <sup>×</sup> <sup>10</sup><sup>7</sup> <sup>C</sup>·kmol−<sup>1</sup> ), R is the universal gas constant, and *T* is the temperature.

The reference exchange current densities *jan*(*T*) and *jcat*(*T*) are dependent on a local temperature, described as follows [52,53].

$$j\_{\rm an}(T) = j\_{\rm an}^{ref} \ e^{-E\_{\rm an}/RT(1 - T/T\_{\rm an}^{ref})} \tag{29}$$

$$j\_{\rm cat}(T) = j\_{\rm cat}^{\rm ref} e^{-E\_{\rm cat}/RT(1 - T/T\_{\rm cat}^{\rm ref})} \tag{30}$$

where *E* and *Tre f* are user-specified activation energy and reference temperature, respectively, and *j re f an* and *j re f cat* are the associated exchange current densities at the specified reference temperature.

The driving force for the kinetics is the local surface overpotential η, also known as the activation loss. It is defined as the difference between the solid and membrane potentials, φ*sol* and φ*m*.

$$
\eta\_{an} = \phi\_{sol} - \phi\_m - \mathcal{U}\_{an}^0 \tag{31}
$$

$$
\eta\_{\rm cat} = \phi\_{\rm sol} - \phi\_m - \mathcal{U}\_{\rm cat}^0 \tag{32}
$$

The half-cell potentials at the cathode and anode, *U*<sup>0</sup> *an* and *U*<sup>0</sup> *cat*, can be calculated using the Nernst equation [50].

$$\mathcal{U}\_{an}^{0} = E\_{an}^{0} - \frac{\Delta S\_{an}}{2F} \left( T - T^{0} \right) - \frac{RT}{2F} \ln \frac{p\_{H\_{2}}}{p^{0}} \tag{33}$$

$$\mathcal{U}\_{\rm cat}^{0} = \mathbb{E}\_{\rm cat}^{0} - \frac{\Delta S\_{\rm cat}}{2F} \Big(T - T^{0}\Big) - \frac{\mathcal{R}T}{2F} \ln \frac{p\_{\rm H\_{2}o}}{p\_{\rm sat} \sqrt{p\_{\rm o} / p^{0}}} \tag{34}$$

where *E* <sup>0</sup> denotes the reversible potential, ∆*S* denotes the reaction entropy, *psat* denotes the saturation pressure of water, *T* <sup>0</sup> and *p* <sup>0</sup> are standard temperature and pressure, and *pH*<sup>2</sup> , *po*<sup>2</sup> , and *pH*2*<sup>o</sup>* are the partial pressures of hydrogen, oxygen, and water vapor, respectively.

In computing the cathode transfer current using Equation (28), the mass transport resistance of the microstructure is not taken into account in the equation [50]. The resistance consists of two parts, which includes the resistance due to ionomer film *Rion* and the resistance due to liquid water film

surrounding the particles *Rliq*. The volumetric transfer current in the cathode layer is represented by the formula below.

$$R\_{\rm cat} = 4F \frac{\mathcal{L}\_{\rm O\_2}}{c\_{\rm O\_2} / \dot{f}\_{\rm O\_2}^{ideal} + R\_{\rm ion} + R\_{\rm liq}} \tag{35}$$

where *cO*<sup>2</sup> is the oxygen concentration at the wall, *Rion* is a user-specified value, and *Rliq* can be calculated as the equation below.

$$R\_{liq} = \frac{\zeta\_{cat} r\_p^2}{K\_w D\_w} \cdot \frac{\sqrt[3]{1 + \frac{\text{sc}}{1 - \epsilon}}}{3(1 - \epsilon)}\tag{36}$$

where ζ*cat* is the specific active surface area for the cathode catalyst [m−<sup>1</sup> ], *s* is the liquid saturation, is the porosity, *r<sup>p</sup>* is the particle diameter, and *KwD<sup>w</sup>* is the product of oxygen solubility and diffusivity in liquid water.

The parameter *j ideal O*2 in Equation (35) is defined as the equation below.

$$j\_{O\_2}^{ideal} = \frac{R\_{cat}^0}{4F} \tag{37}$$

In this case, *R* 0 *cat* is the ideal current transfer, computed using Equation (28) without considering the resistances.

#### 2.1.3. Constitutive Relations

The density of the gas is given by the equation below.

$$
\rho = \frac{p \, M}{R \, T}
$$

where the mixture molecular weight, expressed in terms of molar fraction of individual species *x<sup>i</sup>* , is given by the equation below.

$$M = M\_{\rm O\_2} \mathfrak{x}\_{\rm O\_2} + M\_{\rm H\_2} \mathfrak{x}\_{\rm H\_2} + M\_{\rm H\_2O} \mathfrak{x}\_{\rm H\_2O} + M\_{\rm N\_2} \mathfrak{x}\_{\rm N\_2}$$

The molar fractions are related to the mass fractions shown below.

$$w\_i = \frac{M\_i \mathfrak{x}\_i}{M}$$

Molar concentration of species *i* is defined as:

$$c\_i = \frac{P}{RT} \times x\_i$$

and can be calculated as:

$$c\_l = \mathfrak{x}\_l (c\_{O\_2} + c\_{O\_2} + c\_{H\_2O} + c\_{N\_2})$$

The relative humidity percentage η is a function of water saturation pressure *p sat H*2*O* .

$$\eta = \frac{p\text{x}\_{H\_2O}}{p\_{H\_2O}^{\text{sat}}} \times 100$$

$$p\_{H\_2O}^{\text{sat}} = 101325 \times 10^{c\_1 + c\_2(T - T\_0) + c\_3(T - T\_0)^2 + c\_4(T - T\_0)^3}$$

Assuming the ratio *<sup>x</sup>O*<sup>2</sup> *xN*<sup>2</sup> = <sup>21</sup> <sup>79</sup> , the oxygen molar fraction at the inlet can be determined from the equation below.

$$\mathfrak{x}\_{O\_2}^{in} = 1 - \frac{\mathfrak{x}\_{H\_2O}^{in}}{1 + 79/21}$$

while the mole fraction on the anode side is defined below.

$$\propto\_{H\_2\mathcal{A}} = 1 - \chi\_{H\_2O\mathcal{A}}$$

The inlet velocities at the cathode and anode side are given by the following equations, respectively.

$$v\_a^{in} = \xi\_a^{in} \frac{i\_{avg}}{2FA\_{inlet}} \times \frac{1}{c\_{H\_2}^{in}}$$

$$v\_c^{in} = \xi\_c^{in} \frac{i\_{avg}}{4FA\_{inlet}} \times \frac{1}{c\_{O\_2}^{in}}$$

The average current density *iavg* is given by the equation below.

$$i\_{avg} = \frac{1}{L} \int\_{o}^{L} i \cdot e\_y \, dx$$

where *L* is the length of the fuel cell.

#### 2.1.4. Boundary Conditions

The boundaries of the system as illustrated in Figure 1b are as follows.

I. At the side walls:

$$\mu^{(\mathcal{g})} = 0, \quad \frac{\partial \phi\_i^{(\mathcal{g})}}{\partial \mathbf{x}} = \frac{\partial \phi^{(\mathcal{s})}}{\partial \mathbf{x}} = \frac{\partial \phi^{(m)}}{\partial \mathbf{x}} = \frac{\partial \mathbf{s}}{\partial \mathbf{x}} = \frac{\partial T}{\partial \mathbf{x}} = \mathbf{0}$$

II. At the anode inlet:

$$\dot{m}\_{\mathfrak{a}}^{(\mathfrak{g})} = \dot{m}\_{\mathfrak{a}}^{\mathrm{in}} \; \omega\_{\mathfrak{H}\_{\mathfrak{L}}}^{(\mathfrak{g})} = \omega\_{\mathfrak{H}\_{\mathfrak{L}}\mathfrak{a} \; \prime}^{\mathrm{in}} \; \omega\_{\mathfrak{H}\_{\mathfrak{L}}\mathcal{O}}^{(\mathfrak{g})} = \left. \omega\_{\mathfrak{H}\_{\mathfrak{L}}\mathcal{O}\mathcal{A}}^{\mathrm{in}} \right. \; T = T^{\mathrm{in}} \; \mathrm{s} = 0$$

III. At the cathode inlet:

$$\dot{m}\_{\mathfrak{c}}^{(\mathfrak{g})} = \dot{m}\_{\mathfrak{c}}^{\dot{\mathfrak{m}}} \; \; \omega\_{H\_2}^{(\mathfrak{g})} = \omega\_{H\_2\mathfrak{c}\ \prime}^{\dot{\mathfrak{m}}} \; \; \omega\_{H\_2\mathcal{O}}^{(\mathfrak{g})} = \left. \omega\_{H\_2\mathcal{O}\mathcal{c}}^{\dot{\mathfrak{m}}} \right. \; T = T^{\dot{\mathfrak{m}}} \; \; \mathcal{S} = \mathbf{0}$$

IV. At the outlets:

$$p^{(\mathcal{g})} = p^{\text{ref}} \,, \quad \frac{\partial \omega^{(\mathcal{g})}\_{\text{i}}}{\partial \mathbf{x}} = \frac{\partial \phi^{(\mathbf{s})}}{\partial \mathbf{x}} = \frac{\partial \mathbf{s}}{\partial \mathbf{x}} = \frac{\partial T}{\partial \mathbf{x}} = \mathbf{0}$$

V. At the anode wall terminal:

$$
\phi^{(s)} = 0
$$

VI. At the cathode wall terminal:

$$-\sigma^{(s)}\frac{\partial\phi^{(s)}}{\partial n} = i^{\text{set}}$$

The governing equations together with the constitutive relations and appropriate boundary conditions are then solved numerically.

#### *2.2. Numerical Method*

The developed mathematical model is implemented and customized in the commercial computational software ANSYS Fluent and its PEMFC module together with user-defined functions. The latter allows for changes in constitutive relations, parameters, and—to some extent—the governing equations. Fluent solves all the equations throughout the domain, so variables in the layers that should not be solved are set to zero.

The computational domains (Figure 1b) are created in the commercial software ANSYS Design Modeller and ANSYS Meshing. The whole domain is defined porous, except for the current collectors (cc) and the wall terminals. The whole domain is partitioned into smaller domains for running it in a parallel mode. The computational model is partitioned using the Cartesian Z-direction method to prevent any floating-point exceptions or errors. With convergence criteria of 10−<sup>6</sup> for the residuals of all the conservation equations, iterations are performed, after the mesh independence test to ensure an accurate solution.

#### *2.3. Response Surface Methodology*

In general concern of a process or system involving a response *y* that depends on many other controllable input variables ξ1, ξ2, ξ3, . . . , ξ*<sup>k</sup>* , its relationship with *y* is given by the equation below [54].

$$y = f(\xi\_1, \xi\_2, \xi\_3, \dots, \xi\_k) + \epsilon$$

where represents other sources for variability unaccounted for in *f*. Treating as a statistical error, we assume it to have a normal distribution with a mean zero and variance σ 2 .

#### *2.4. Kriging-Based Response Surface Methodology*

Design and analysis of computer experiments (DACE) is also called the kriging method of response surface generation. It involves training points in estimating the unknown parameters α and predicting new response points. It interpolates the model at all the training points. This method is used to generate response surfaces for voltage values of 0.6 and 0.8 V. The response can be expressed by the equation below [54].

$$
\hat{y}(\mathbf{x}) = f(\mathbf{x}) + z(\mathbf{x})
$$

where *f*(*x*) is a low-order polynomial that interpolates the design points. Typically, a constant value is found in order to predict for modelling complex input-output relations. Hence, the output can be viewed as:

$$
\mathfrak{z}(\mathfrak{x}) = \mathfrak{z} + z(\mathfrak{x}).
$$

*z*(*x*) is a Gaussian stochastic function representing the realization of a random process with zero mean, variance σ 2 , and its covariance is given by the equation below.

$$\operatorname{Cov}(Z) = \sigma^2 \mathbf{R} \{ \mathbf{x}^i, \mathbf{x}^j \}$$

where *R x i* , *x j* is the correlation matrix, defined by the equation below.

$$\begin{aligned} \mathcal{R}\left(\mathbf{x}^i, \mathbf{x}^j\right) &= \exp\left[-d\left(\mathbf{x}^i, \mathbf{x}^j\right)\right] \\\\ d\left(\mathbf{x}^i, \mathbf{x}^j\right) &= \sum\_{l=1}^k \partial\_l \left(\left\|\mathbf{x}^i\_l - \mathbf{x}^j\_l\right\|^p\right) \end{aligned}$$

where *i*, *j* denote the two training points, *l* refers to the design parameter, θ is the positive weight factor related to each design parameter, and *k* denotes the number of design parameters. The mean parameter β is evaluated by the equation below.

$$\beta = \left[A^T \, \mathcal{R}^{-1} \, A\right]^{-1} A^T \, \mathcal{R}^{-1} y$$

where *A* is an *n* × *n* matrix of training set points depending on the choice of the function *f*(*x*). The parameters θ and *p* ensure best fit to the training data. They are evaluated by using a maximizing likelihood estimation (MLE) [55].

$$-\frac{1}{2} \left[ n \ln(2\pi) + n \ln \sigma^2 + \ln|R| + \frac{1}{2\sigma^2} (y - A\,\beta)^T R^{-1} (y - A\beta) \right]^2$$

where the maximum likelihood σ 2 is expressed by the equation below.

$$
\sigma^2 = \frac{1}{n} (y - A\beta)^T \mathbb{R}^{-1} (y - A\beta)
$$

The response at a new point *x* , *y*ˆ(*x*) is directly evaluated by applying the equation below.

$$\hat{y}(\overline{\mathbf{x}}) = \boldsymbol{\beta} + \boldsymbol{r}^T(\overline{\mathbf{x}}) \boldsymbol{R}^{-1} (\boldsymbol{y} - \boldsymbol{A}\boldsymbol{\beta})$$

where *r*(*x*) is a correlation vector between *x* and all the training points.

#### **3. Results and Discussion**

#### *3.1. Validation*

Due to the complexity of the model being solved, model validation with experimental data is imperative to prevent misleading conclusions in predicting the behavior of the fuel cell system. In this work, we aimed to validate our model with two experimental single cells. The first experimental cell was equipped with a Gore Primea 5510 membrane, which is a microscopically reinforced composite membrane. The expressions for various phenomenological membrane models are generally based on the Nafion membrane and so need to be adapted to account for the Gore membrane used in the experiment. Two parameters have been, therefore, adapted to validate the model using both the polarization curve and its iR-corrected counterpart.

The iR-corrected potential is given by the equation below.

$$E\_{IR} = E\_{rev} - \eta\_a - \eta\_c$$

where η*<sup>a</sup>* (> 0) and η*<sup>c</sup>* (< 0) are the corresponding overpotentials of the anode and the cathode catalyst layers, respectively. In this case, we focused on parameter adaption of the cathode reference exchange current density, *j re f cat* , and the cathode transfer coefficient, α*cat* , and retained the anode counterparts from the work of Wang and co-workers [39,56].

Two points from the experimentally determined iR-corrected polarization curve were chosen: one at a low current density and the other at a higher current density. Initially, the membrane protonic conductivity coefficient, β*m*, was set to one. Once a good agreement for the two points was obtained, we predicted the complete iR-corrected polarization curve. A good agreement for the whole range was achieved, as can be seen in Figure 2. Subsequently, β*<sup>m</sup>* was varied to finish the validation for the full polarization curve, where the cell voltage is defined below.

$$E\_{cell} = E\_{rev} - \eta\_a - \eta\_c - \sum i^{(s)} \left(\frac{1}{\sigma\_{eff}}\right) - \sum i^{(m)} \left(\frac{1}{\sigma\_{eff}}\right)^2$$

In this case, the last two terms on the right-hand side of the equation account for the various ohmic losses in the solid and membrane functional layers, respectively.

*Molecules* **2019**, *24*, x 14 of 35

is the iR-corrected experimental potential,

**Figure 2.** Polarization curves: [ ] is the experimental potential, [ ] is the iR-corrected experimental potential, [ ] is the experimental power density, [— · ] is the predicted potential, [ — ] is the predicted iR-corrected potential, and [···] is the predicted power density. **Figure 2.** Polarization curves: [ ] is the experimental potential, [ ] is the iR-corrected experimental potential, [ is the experimental power density, [— · ] is the predicted potential, [ — ] is the predicted iR-corrected ] is the experimental power density, [— · ] is the predicted potential, [ — ] is the predicted iR-corrected potential, and [··· ] is the predicted power density.

Furthermore, the predicted local current density distribution along the top of the anode terminal was compared with the experimental counterpart. Figure 3 illustrates this local current density distribution for each measured average current density in Figure 2. It can be observed that better prediction is achieved at lower currents, while the model is less accurate at higher current densities, with the most deviation observed near the inlet and outlet. This could be due to the inlet boundary condition in the simulation that does not represent correctly the position of the inlet manifolds location as in the experiments. potential, and [···] is the predicted power density. Furthermore, the predicted local current density distribution along the top of the anode terminal was compared with the experimental counterpart. Figure 3 illustrates this local current density distribution for each measured average current density in Figure 2. It can be observed that better prediction is achieved at lower currents, while the model is less accurate at higher current densities, with the most deviation observed near the inlet and outlet. This could be due to the inlet boundary condition in the simulation that does not represent correctly the position of the inlet manifolds location as in the experiments. Furthermore, the predicted local current density distribution along the top of the anode terminal was compared with the experimental counterpart. Figure 3 illustrates this local current density distribution for each measured average current density in Figure 2. It can be observed that better prediction is achieved at lower currents, while the model is less accurate at higher current densities, with the most deviation observed near the inlet and outlet. This could be due to the inlet boundary condition in the simulation that does not represent correctly the position of the inlet manifolds location as in the experiments. *Molecules* **2019**, *24*, x 15 of 28 is the experimental power density, [— · ] is the predicted potential, [ — ] is the predicted iR-corrected potential, and [···] is the predicted power density. Furthermore, the predicted local current density distribution along the top of the anode terminal was compared with the experimental counterpart. Figure 3 illustrates this local current density distribution for each measured average current density in Figure 2. It can be observed that better prediction is achieved at lower currents, while the model is less accurate at higher current densities, with the most deviation observed near the inlet and outlet. This could be due to the inlet boundary condition in the simulation that does not represent correctly the position of the inlet manifolds

location as in the experiments.

**Figure 3.** Local current density distribution. Predicted (lines) and experimental (points) for current densities: 1.3 Acm<sup>−</sup><sup>2</sup> ( ), 1.1 A cm<sup>−</sup><sup>2</sup> ( ), 0.9 A cm<sup>−</sup><sup>2</sup> ( ), 0.7 A cm<sup>−</sup><sup>2</sup> ( ), 0.5 A cm<sup>−</sup><sup>2</sup> ( ), 0.4 A cm<sup>−</sup><sup>2</sup> ( ), 0.3 A cm<sup>−</sup><sup>2</sup> ( ), 0.2 A cm<sup>−</sup><sup>2</sup> ( ), 0.1 A cm<sup>−</sup><sup>2</sup> ( ), and 0.05 A cm<sup>−</sup><sup>2</sup> ( ). **Figure 3.** Local current density distribution. Predicted (lines) and experimental (points) for current densities: 1.3 Acm−<sup>2</sup> ( densities: 1.3 Acm<sup>−</sup><sup>2</sup> ( ), 1.1 A cm−<sup>2</sup> ( ), 0.9 A cm−<sup>2</sup> ( ), 0.7 A cm−<sup>2</sup> ( ), 0.7 A cm<sup>−</sup><sup>2</sup> ( ), 0.5 A cm−<sup>2</sup> ( ), 0.5 A cm<sup>−</sup><sup>2</sup> ( ), 0.4 A cm−<sup>2</sup> ( ), 0.4 A cm<sup>−</sup><sup>2</sup> ( ), 0.3 A cm−<sup>2</sup> ( ), 0.2 A cm−<sup>2</sup> ( ), 0.2 A cm<sup>−</sup><sup>2</sup> ( ), 0.1 A cm−<sup>2</sup> ( ), and 0.05 A cm−<sup>2</sup> ( ).

), 1.1 A cm<sup>−</sup><sup>2</sup> (

phase liquid water.

The model is also validated against the second set of experimental polarization curve, especially

.

phase liquid water.

density of 1 A/cm<sup>2</sup>

2.3 (

The model is also validated against the second set of experimental polarization curve, especially at a limiting current density to justify the MEA model with an agglomerate catalyst layer. The model was validated with experimental fuel cells with a single-layered gas diffusion studied by Han et al. [37], as depicted in Figure 4. It shows that the model has good agreement with the curve and is able to predict the limiting current density due to mass transport limitations and the presence of two-

**Figure 4.** Experimental (points) and predicted (line) polarization curve for case b.

To determine the dominant parameters affecting the PEMFC performance and optimize them, response surface generation and local sensitivity analysis were carried out at medium (0.6 V) and high (0.8 V) voltages. The response surfaces were created using the design of experiments (DOE) method of central composite design (CCD) and the kriging method of response surface generation. The CCD was employed to capture the non-linear interactions that cannot otherwise be described by linear functions. Hence, experimental designs for quadratic response surfaces, like three-level factorial, Box-Behnken, central composite, and Doehlert designs [57], should be used instead. The list of design variables considered for the DOE and their upper and lower bounds are tabulated in Table 4. The base case parameters correspond to case (a) in Table 1, which has been validated for both global and local current densities (Figures 2 and 3). The inputs used for the base-case correspond to a current

 = 1.052 m/s), which are calculated from constitutive relations, as explained in the previous section. The response surface generated was then carefully evaluated. If the response surface was not within the desired limits of accuracy, determined by the "goodness of fit," it was modified by adding

= 0.173 m/s) and cathode stoichiometry of

*3.2. Response Surface Generation and Optimization*

, i.e., anode stoichiometry of 3.35 (

The model is also validated against the second set of experimental polarization curve, especially at a limiting current density to justify the MEA model with an agglomerate catalyst layer. The model was validated with experimental fuel cells with a single-layered gas diffusion studied by Han et al. [37], as depicted in Figure 4. It shows that the model has good agreement with the curve and is able to predict the limiting current density due to mass transport limitations and the presence of two-phase liquid water. *Molecules* **2019**, *24*, x 16 of 28

**Figure 4.** Experimental (points) and predicted (line) polarization curve for case b. **Figure 4.** Experimental (points) and predicted (line) polarization curve for case b.

#### *3.2. Response Surface Generation and Optimization 3.2. Response Surface Generation and Optimization*

To determine the dominant parameters affecting the PEMFC performance and optimize them, response surface generation and local sensitivity analysis were carried out at medium (0.6 V) and high (0.8 V) voltages. The response surfaces were created using the design of experiments (DOE) method of central composite design (CCD) and the kriging method of response surface generation. The CCD was employed to capture the non-linear interactions that cannot otherwise be described by linear functions. Hence, experimental designs for quadratic response surfaces, like three-level factorial, Box-Behnken, central composite, and Doehlert designs [57], should be used instead. The list of design variables considered for the DOE and their upper and lower bounds are tabulated in Table 4. The base case parameters correspond to case (a) in Table 1, which has been validated for both global and local current densities (Figures 2 and 3). The inputs used for the base-case correspond to a current density of 1 A/cm<sup>2</sup> , i.e., anode stoichiometry of 3.35 ( = 0.173 m/s) and cathode stoichiometry of 2.3 ( = 1.052 m/s), which are calculated from constitutive relations, as explained in the previous section. The response surface generated was then carefully evaluated. If the response surface was not within the desired limits of accuracy, determined by the "goodness of fit," it was modified by adding refinement points to the kriging method. Additional CFD simulations were run to generate the response surface data to improve the confidence of the response surface. This iterative method was done until a good fit and reliable response surface was achieved and is illustrated in Figure 5. To determine the dominant parameters affecting the PEMFC performance and optimize them, response surface generation and local sensitivity analysis were carried out at medium (0.6 V) and high (0.8 V) voltages. The response surfaces were created using the design of experiments (DOE) method of central composite design (CCD) and the kriging method of response surface generation. The CCD was employed to capture the non-linear interactions that cannot otherwise be described by linear functions. Hence, experimental designs for quadratic response surfaces, like three-level factorial, Box-Behnken, central composite, and Doehlert designs [57], should be used instead. The list of design variables considered for the DOE and their upper and lower bounds are tabulated in Table 4. The base case parameters correspond to case (a) in Table 1, which has been validated for both global and local current densities (Figures 2 and 3). The inputs used for the base-case correspond to a current density of 1 A/cm<sup>2</sup> , i.e., anode stoichiometry of 3.35 (*v in <sup>a</sup>* = 0.173 m/s) and cathode stoichiometry of 2.3 (*v in <sup>c</sup>* = 1.052 m/s), which are calculated from constitutive relations, as explained in the previous section. The response surface generated was then carefully evaluated. If the response surface was not within the desired limits of accuracy, determined by the "goodness of fit," it was modified by adding refinement points to the kriging method. Additional CFD simulations were run to generate the response surface data to improve the confidence of the response surface. This iterative method was done until a good fit and reliable response surface was achieved and is illustrated in Figure 5.

Protonic conduction coefficient of the membrane 0.9 0.5 1.5 Hydrophobic angle of the cathode catalyst layer 95° 90° 180° Ionomer tortuosity of the cathode catalyst layer 1 0.7 1.5 Ionomer volume fraction of the cathode catalyst layer 1 0.5 1.0

Proton Exchange Membrane thickness 0.03 mm 0.005 mm 0.05 mm Equivalent weight of membrane 1100 kg kmol<sup>−</sup><sup>1</sup> 700 kg kmol<sup>−</sup><sup>1</sup> 1500 kg kmol<sup>−</sup><sup>1</sup> Radius of cathode catalyst particle 10<sup>−</sup><sup>7</sup> m 10<sup>−</sup><sup>8</sup> m 10<sup>−</sup><sup>7</sup> m

**Table 4.** Design variable values for the DOE method.

**Design Variable Base Case Lower Bound Upper Bound**


**Table 4.** Design variable values for the DOE method.

**Figure 5.** Flow diagram of response surface methodology. **Figure 5.** Flow diagram of response surface methodology.

Goodness of fit was used to determine the reliability of the response surface predicted, i.e., how close the predicted values are to the observed values. The predicted values from the response surface were compared against the values observed from design points for both 0.6 and 0.8 V. The results demonstrate good agreement for both voltage values, which is shown in Figure 6. Goodness of fit was used to determine the reliability of the response surface predicted, i.e., how close the predicted values are to the observed values. The predicted values from the response surface were compared against the values observed from design points for both 0.6 and 0.8 V. The results demonstrate good agreement for both voltage values, which is shown in Figure 6.

**Figure 6.** Goodness of fit for response surface generated at: (**a**) 0.8 V and (**b**) 0.6 V. **Figure 6.** Goodness of fit for response surface generated at: (**a**) 0.8 V and (**b**) 0.6 V.

#### *3.3. Local Sensitivity 3.3. Local Sensitivity*

In determining which variables influence the output current density the most, the following relation was used. In determining which variables influence the output current density the most, the following relation was used.

$$local\ sensitivity = \frac{Output\_{\text{max}} - Output\_{\text{min}}}{Output\_{\text{average}}}$$

 The local sensitivity analysis was carried out using outputs obtained from the DOE. Figures 7a and 8a illustrate the relative impacts of the different input parameters on the local sensitivity for the 0.6 and 0.8 V, respectively. The corresponding sensitivity curves (Figures 7b and 8b) show the output variation with changes in one input parameter, while keeping the other parameters constant. The results show that, for both voltage values, parameters that have the most impact are the membrane protonic conductivity coefficient and the membrane thickness. The cathode catalyst ionomer volume fraction, cathode catalyst porosity, cathode catalyst ionomer tortuosity, volume fraction, cathode catalyst hydrophobic angle, and the radius of cathode agglomerate particles have a minor impact on The local sensitivity analysis was carried out using outputs obtained from the DOE. Figures 7a and 8a illustrate the relative impacts of the different input parameters on the local sensitivity for the 0.6 and 0.8 V, respectively. The corresponding sensitivity curves (Figures 7b and 8b) show the output variation with changes in one input parameter, while keeping the other parameters constant. The results show that, for both voltage values, parameters that have the most impact are the membrane protonic conductivity coefficient and the membrane thickness. The cathode catalyst ionomer volume fraction, cathode catalyst porosity, cathode catalyst ionomer tortuosity, volume fraction, cathode catalyst hydrophobic angle, and the radius of cathode agglomerate particles have a minor impact on the output that can be considered negligible when compared to these two parameters. *Molecules* **2019**, *24*, x 19 of 28

the output that can be considered negligible when compared to these two parameters.

**Figure 7. Figure 7.** Local sensitivity analysis at 0.6 V presented in: ( Local sensitivity analysis at 0.6 V presented in: **a**) bar plot and ( (**a**) bar plot and **b** ( ) sensitivity curve. **b**) sensitivity curve.

**Figure 8.** Local sensitivity analysis at 0.8 V presented in: (**a**) the bar plot and the (**b**) sensitivity **Figure 8.** Local sensitivity analysis at 0.8 V presented in: (**a**) the bar plot and the (**b**) sensitivity curve.

curve. Based on these results, the membrane thickness and the membrane protonic conductivity coefficient are chosen as the varying parameters to study the interrelated response, as well as to perform the optimization. The three-dimensional response charts illustrating how the two variables affect the current output are depicted in Figures 9 and 10 for both medium and high voltages, Based on these results, the membrane thickness and the membrane protonic conductivity coefficient are chosen as the varying parameters to study the interrelated response, as well as to perform the optimization. The three-dimensional response charts illustrating how the two variables affect the current output are depicted in Figures 9 and 10 for both medium and high voltages, respectively. These response surfaces are used for the system optimization. *Molecules* **2019**, *24*, x 21 of 28

respectively. These response surfaces are used for the system optimization.

**Figure 9. Figure 9.** Three Three-dimensional response surface plot at 0.6 V. -dimensional response surface plot at 0.6 V.

**Figure 10. Figure 9.** Three Three-dimensional response surface plot at 0.8 V. -dimensional response surface plot at 0.6 V.

#### *3.4. Optimization*

Optimization was carried out using the non-linear programming through a quadratic Lagrangian (NLPQL) approach based on the response surfaces generated previously. NLPQL is a gradient-based algorithm that provides a refined local optimization result. It supports a single constraint on the output parameter and is limited to continuous parameters. The NLPQL approximates derivatives by a central difference scheme and finds candidate points by iterations. This approach was used to determine the optimum values of the variables considered.

Table 5 shows the optimization result that provides maximum current output for three candidate points at 0.8 V. It can be observed that the response surface prediction agrees well with the CFD model, with a maximum error of approximately 0.17%, occurred at candidate point 1 (case i). The current density output increased from 0.18 A cm−<sup>2</sup> at the base to an average of 0.2472 A cm−<sup>2</sup> , which indicates almost a 40% increase. This increase is attributed to the change in the protonic conductivity coefficient from 0.9 to 1.45 and attributed to the decrease in membrane thickness from 0.03 mm to the lowest bound value of 0.005 mm.

Similarly, a good agreement between the response surface and the computational fluid dynamics (CFD) model is obtained for the medium-voltage case, as shown in Table 6. At 0.6 V, the error is slightly higher, with a maximum error of 0.57%, which occurred at candidate point 1. As in the case at a high voltage value, the protonic conduction coefficient increased from 0.9 to 1.5 m and the membrane thickness reduced to the lowest bound value. However, the current density output for the 0.6 V almost doubled, from 1.23 A cm−<sup>2</sup> at the base case to 2.4 A cm−<sup>2</sup> , which gives approximately a 96% increase. This large improvement demonstrates the significance of optimization at medium voltage levels, as compared to high voltage values.


*Molecules* **2019**, *24*, 3097

It can be observed that the extent of optimization required increased greatly when the voltage is reduced from a high to a medium level. This is due to the fact that an increase in current density requires better conductivity of the H<sup>+</sup> ions through the membrane to maintain system net neutrality. In addition, the reduction in membrane thickness decreases the resistance offered by the membrane for the hydrogen ions to move from the anode to the cathode side. The sensitivity analysis also shows that the influence of membrane thickness on the output is greater than that of the protonic conductivity coefficient and rises slightly when the voltage is reduced from 0.8 to 0.6 V. Furthermore, the different values of parameters, other than the two mentioned, in Tables 5 and 6, confirms the insignificance of these variables to the current output and, hence, varying them would not affect the output largely.

In short, it can be deduced that, in designing high performance PEMFC, one needs to aim for thinner MEA with higher membrane protonic conductivity, which can be achieved by using carbon-reinforced membrane or water absorbent materials including polytetrafluoroethylene (PTFE), polyvinylidene fluoride (PVDF), and fluorinated ethylene propylene (FEP) or silica gels.

#### **4. Conclusions**

A numerical study of the two-phase PEMFC with a detailed multiscale agglomerate catalyst layer model was developed and validated against two sets of experimental data for iR-corrected and full polarization curves, including at limiting current densities due to mass transport limitation and the local current density distributions. The model is then extended and coupled with response surface methodology to optimize the design of the membrane electrode assembly (MEA), i.e., membrane thickness, equivalent weight of membrane, radius of agglomerate catalyst particle, cathode catalyst ionomer resistance, porosity of the catalyst layer, a membrane protonic conductivity coefficient, a hydrophobic angle of the catalyst layer, ionomer tortuosity, and a catalyst layer, at high and medium voltages. From sensitivity analysis, it was found that the membrane thickness and membrane protonic conductivity coefficient yield the most significant factor. Reducing the membrane thickness by 40% and increasing protonic conductivity by 50% gives rise to a current density of up to 40% at a higher voltage and up to 100% at a medium voltage. This finding could help fuel cell engineers and designers to carefully manufacture MEA with optimum parameters for a high-performance fuel cell system.

Future work will focus on a combined optimization of design and operating parameters simultaneously for better MEA design, thermal, water, and gas management. A more advanced optimization algorithm including artificial intelligence and machine learning will be considered as well.

**Author Contributions:** Conceptualization, A.P.S. and M.A. Methodology, A.P.S. and R.K.S.S.V. Software, R.K.S.S.V. and B.A.C. Validation, A.P.S., R.K.S.S.V. and B.A.C. Formal analysis, R.K.S.S.V., B.A.C., A.P.S. and M.A. Investigation, R.K.S.S.V., B.A.C. and A.P.S. Resources, A.P.S., M.A., L.J. and L.C. Data curation, A.P.S. and M.A. Writing—original draft preparation, R.K.S.S.V. and B.A.C. Writing—review and editing, A.P.S., M.A., L.J. and L.C. Visualization, R.K.S.S.V. and B.A.C. Supervision, A.P.S. and M.A. Project administration, L.J. and L.C. Funding acquisition, A.P.S., L.J. and L.C.

**Funding:** This research received no external funding.

**Acknowledgments:** R.K.S.S.V. and B.A.C. gratefully acknowledge support from MITACS Globalink Research Award. A.P.S., L.J., and L.C. acknowledge the support from the SDUST Open Grant. We acknowledge anonymous reviewers for the valuable feedback and suggestions.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**




#### **References**


#### **Sample Availability:** Not Available.

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*
