Multi-Variate Optimization of Polymer Electrolyte Membrane Fuel Cells in Consideration of Effects of GDL Compression and Intrusion

This numerical study is based on the following specific assumptions

• 1.  
  The operation pressure is low, and hence, ideal gas mixtures are considered in the gas phase.
• 2.  
  The flow velocity is low and laminar.
• 3.  
  The gravitation effect is negligible.
• 4.  
  The effect of immobile liquid saturation in the porous region is negligible.

The PEMFC model considered herein is governed by five principles of conservation: mass, species, charge, momentum, and thermal energy. The above-stated equations are coupled with source terms related to electrochemical reactions in the anode (hydrogen oxidation reaction; HOR) and cathode (oxygen reduction reaction; ORR). The governing equations and source terms are listed in Table I and II, respectively, for reference.

Table I. PEMFC model: governing equations.

Governing equations
Mass	${\rm{\nabla }}\cdot \left(\rho \vec{u}\right)=0$	(37)
Momentum	$\left(\displaystyle \frac{1}{{\varepsilon }^{2}}\right){\rm{\nabla }}\cdot \left(\rho \vec{u}\vec{u}\right)=-{\rm{\nabla }}P+{\rm{\nabla }}\cdot \tau +{S}_{u}$	(38)
Species	Flow channels and porous media:
	$\begin{array}{l}{\rm{\nabla }}\cdot \left({\gamma }_{i}\rho {m}_{i}\vec{u}\right)={\rm{\nabla }}\cdot \left[{\rho }^{g}{D}_{i}^{g,eff}{\rm{\nabla }}\left({m}_{i}^{g}\right)\right]\\ \,+{\rm{\nabla }}\cdot \left[\left({m}_{i}^{g}-{m}_{i}^{l}\right){\vec{j}}^{l}\right]+{S}_{i}\end{array}$	(39)
	Water transport in membrane:
	$\begin{array}{l}{\rm{\nabla }}\cdot \left(\left(\displaystyle \frac{{\rho }^{{\rm{mem}}}}{EW}\right){{\rm{D}}}_{{\rm{w}}}^{{\rm{mem}}}{\rm{\nabla }}\lambda \right){M}_{w}\\ \,-{\rm{\nabla }}\cdot \left({n}_{d}\left(\displaystyle \frac{I}{F}\right)\right){M}_{w}+{\rm{\nabla }}\cdot \left(\left(\displaystyle \frac{{\kappa }^{mem}}{{\nu }^{l}}\right){\rm{\nabla }}{P}^{l}=0\right.\end{array}$	(40)
Charge	Proton transport:
	${\rm{\nabla }}\cdot \left({\kappa }^{{\rm{eff}}}{\rm{\nabla }}{\phi }_{e}\right)+{S}_{\phi }=0$	(41)
	Electron transport:
	${\rm{\nabla }}\cdot \left({\sigma }^{{\rm{eff}}}{\rm{\nabla }}{\phi }_{s}\right)-{S}_{\phi }=0$	(42)
Energy	${\rm{\nabla }}\cdot \left(\rho \vec{u}{C}_{p}^{g}T\right)={\rm{\nabla }}\cdot \left({k}^{eff}{\rm{\nabla }}T\right)+{S}_{T}$	(43)

Table II. PEMFC model: Source/sink terms.

Description	Expressions
Momentum	Porous media	${{\rm{S}}}_{{\rm{u}}}=-\displaystyle \frac{\mu }{K}\vec{u}$
Species	H2 in anode CL	${{\rm{S}}}_{{H}_{2},a}=\left[-\displaystyle \frac{j}{2F}\right]{{\rm{M}}}_{{H}_{2}}$
	O2 in cathode CL	${{\rm{S}}}_{{O}_{2},c}=\left[\displaystyle \frac{j}{4F}\right]{{\rm{M}}}_{{O}_{2}}$
	Water in anode CL	${{\rm{S}}}_{w,a}=\left[-{\rm{\nabla }}\cdot \left(\displaystyle \frac{{n}_{d}}{F}\vec{I}\right)+\displaystyle \frac{j}{4F}\right]{{\rm{M}}}_{w}$
	Water in cathode CL	${{\rm{S}}}_{w,c}=\left[-{\rm{\nabla }}\cdot \left(\displaystyle \frac{{n}_{d}}{F}\vec{I}\right)-\displaystyle \frac{j}{2F}\right]{{\rm{M}}}_{w}$
Energy	In anode CL	${{\rm{S}}}_{{\rm{T}},{\rm{a}}}=j\cdot \eta +\displaystyle \frac{{I}^{2}}{{k}^{eff}}$
	In cathode CL	${{\rm{S}}}_{{\rm{T}},{\rm{c}}}=j\left(\eta +T\displaystyle \frac{d{U}_{0}}{dT}\right)+\displaystyle \frac{{I}^{2}}{{k}^{eff}}$
	In membrane	${{\rm{S}}}_{{\rm{T}}}=\displaystyle \frac{{I}^{2}}{{k}^{eff}}$
Charge	In CLs:	${{\rm{S}}}_{\phi }=j$
Electrochemical reactions $\displaystyle \sum _{k}{s}_{i}{M}_{i}^{z}=n{e}^{-},{\rm{where}}\,\left\{\begin{array}{c}{M}_{i}=chemical\,formula\,of\,species\,i\\ {s}_{i}=stoichiometric\,coefficient\,\\ \,n=number\,of\,electrons\,transferred\end{array}\right.$
HOR on the anode side: ${{\rm{H}}}_{2}-2{{\rm{H}}}^{+}=2{{\rm{e}}}^{-}$
ORR on the cathode side: $2{{\rm{H}}}_{2}{\rm{O}}-{{\rm{O}}}_{2}-4{{\rm{H}}}^{+}=4{{\rm{e}}}^{-}$
Transfer current	HOR in anode CL:
density, $\left[A/{m}^{3}\right]$	$\begin{array}{l}j=\left(1-s\right)a{i}_{0,a}^{ref}{\left(\displaystyle \frac{{C}_{{H}_{2}}}{{C}_{{H}_{2},ref}}\right)}^{\displaystyle \frac{1}{2}}\,\exp \left[-\displaystyle \frac{{E}_{a}}{R}\left(\displaystyle \frac{1}{T}-\displaystyle \frac{1}{353.15}\right)\right]\\ \,\times \,\left(\exp \left(\displaystyle \frac{{\alpha }_{a}F}{RT}\eta \right)-\exp \left(-\displaystyle \frac{{\alpha }_{c}F}{RT}\eta \right)\right)\end{array}$	(44)
	ORR in cathode CL:
	$\begin{array}{l}j=-\displaystyle \frac{3{L}_{Pt}}{{r}_{Pt}\cdot {\rho }_{Pt}\cdot {\delta }_{CL}}{i}_{0,c}^{ref}{\left(\displaystyle \frac{{C}_{{O}_{2}}^{Pt}}{{C}_{{O}_{2},ref}}\right)}^{3/4}\exp \left(-\displaystyle \frac{{E}_{c}}{R}\left(\displaystyle \frac{1}{T}-\displaystyle \frac{1}{353.15}\right)\right)\exp \left(-\displaystyle \frac{{\alpha }_{c}}{{R}_{u}T}F\eta \right)\end{array}$	(45)
Overpotential	$\eta =\phi {\,}_{s}-{\phi }_{e}-U$	(46)
	where $U=\,{U}_{0}-\,\displaystyle \frac{RT}{nF}\,\mathrm{ln}\,\displaystyle \frac{{C}_{{O}_{2}}}{{C}_{{O}_{2\,ref}}}$

Additionally, the kinetic, transport, and physiochemical properties of the PEMFC components are listed in Table III. The M² model proposed by Wang and Cheng ³⁰ is employed herein, and the mixture properties are listed in Table IV. Table V lists a set of the species' transport properties, and they are correlated to the water content $\lambda ,$ which subsequently is a function of the water activity $a.$ ³¹

Table III. Kinetic, transport, and physiochemical properties.

Description	Value/Expression	References
Exchange current density of HOR × ECSA per unit CL volume, $a{I}_{0,{\rm{a}}}^{ref}$	1.2 × 1010 A m−3	24
Exchange current density for ORR, ${I}_{0,{\rm{c}}}^{ref}$	2.0 × 10−4 A cm−2-P−1t−1	11
Activation energy of anode, ${E}_{a}$	10.0 kJ mol−1	24
Activation energy of cathode, ${E}_{c}$	70.0 kJ mol−1	24
Transfer coefficient of HOR, ${\alpha }_{a}={\alpha }_{c}$	1	25
Transfer coefficient of ORR, ${\alpha }_{c}$	1	25
Reference H2/O2 molar concentration, ${C}^{ref}$	40.88 mol m−3	25
Permeability of GDL/CL, ${K}_{GDL}/{K}_{CL}$	1.0 × 10−12/1.0 × 10−13 m2	24
Equivalent weight of electrolyte in the membrane, $EW$	1.1 kg mol−1	24
Young's modulus of GDL	6.16 MPa	42
Poisson ratio of GDL	0.09	29
Faraday's constant, ${\rm{F}}$	96,485 C mol−1	25
Universal gas constant, ${R}_{u}$	8.314 $J/\left(mol\cdot K\right)$
H2 diffusivity in the anode gas channel, ${D}_{0,{H}_{2},a}^{g}$	1.1028 × 10−4 m2 s−1	25
H2O diffusivity in the anode gas channel, ${D}_{0{H}_{2}O,a}^{g}$	1.1028 × 10−4 m2 s−1	25
O2 diffusivity in the cathode gas channel, ${D}_{0,{O}_{2},c}^{g}$	3.2348 × 10−4 m2 s−1	25
H2O diffusivity in the cathode gas channel, ${D}_{0,{H}_{2}O,c}^{g}$	7.35 × 10−5 m2 s−1	25
Binary gas diffusivity (${D}_{i}^{g}$)	For nonporous regions
	${D}_{i}^{g}={D}_{o}^{g}{\left(\displaystyle \frac{T}{{T}_{0}}\right)}^{3/2}\left(\displaystyle \frac{{P}_{0}}{P}\right)$	(47)
Effective diffusivity (${D}_{i}^{g,eff}$)	For porous regions
	${D}_{i}^{g,eff}=\displaystyle \frac{\varepsilon }{{\tau }_{k}}{D}_{i}^{g}$	(48)

Table IV. Expressions used in the two-phase mixture model.

Description	Expression
Mixture density ($\rho$)	${\rm{\rho }}={{\rm{\rho }}}^{l}s+{\rho }^{g}\left(1-s\right)$	(49)
Gas mixture density $\left({\rho }^{g}\right)$	${\rho }^{g}=\left(\displaystyle \frac{P}{{R}_{u}T}\right)\displaystyle \frac{1}{{\sum }_{i}\displaystyle \frac{{m}_{i}^{g}}{{M}_{i}}}$	(50)
Mixture velocity ($\rho \vec{u}$)	$\rho \vec{u}={\rho }^{l}{\vec{u}}^{l}+{\rho }^{g}{\vec{u}}^{g}$	(51)
Mixture mass fraction	${m}_{i}=\displaystyle \frac{{\rho }^{l}s{m}_{i}^{l}+{\rho }^{g}\left(1-s\right){m}_{i}^{g}}{\rho }$	(52)
Relative permeability	${k}_{r}^{l}={s}^{3}$	(53)
	${k}_{r}^{g}={(1-s)}^{3}$	(54)
Kinematic viscosity of the two-phase mixture	$v={\left(\displaystyle \frac{{k}_{r}^{l}}{{v}^{l}}+\displaystyle \frac{{k}_{r}^{g}}{{v}^{g}}\right)}^{-1}$	(55)
Kinematic viscosity of the gas mixture	${v}^{g}=\displaystyle \frac{{\mu }^{g}}{{\rho }^{g}}=\displaystyle \frac{1}{{\rho }^{g}}\displaystyle \sum _{n}^{i=1}\displaystyle \frac{{x}_{i}{\mu }_{i}}{{\sum }_{j=1}^{n}{x}_{j}{\phi }_{ij}}$	(56)
	where ${\phi }_{ij}=\displaystyle \frac{1}{\sqrt{8}}{\left(1+\displaystyle \frac{{M}_{i}}{{M}_{j}}\right)}^{-1/2}{\left[1+{\left(\displaystyle \frac{{\mu }_{i}}{{\mu }_{j}}\right)}^{1/2}{\left(\displaystyle \frac{{M}_{j}}{{M}_{i}}\right)}^{1/4}\right]}^{2}$	(57)
	and
	${\mu }_{i}\left[N.s.{m}^{-2}\right]=\left\{\begin{array}{c}{\mu }_{{H}_{2}}=0.21\times {10}^{-6}{T}^{0.66}\\ {\mu }_{w}=0.00584\times {10}^{-6}{T}^{1.29}\\ {\mu }_{{N}_{2}}=0.237\times {10}^{-6}{T}^{0.76}\\ {\mu }_{{O}_{2}}=0.246\times {10}^{-6}{T}^{0.78}\end{array}\right.,$ T in kelvin	(58)
Relative mobility	${\lambda }^{l}=\displaystyle \frac{{k}_{r}^{l}}{{v}^{l}}v$	(59)
	${\lambda }^{g}=1-{\lambda }^{l}$	(60)
Diffusive mass flux	$\mathop{{j}^{l}}\limits^{\longrightarrow}={\rho }^{l}{\vec{u}}^{l}-{\lambda }^{l}\rho \vec{u}=\displaystyle \frac{K}{v}{\lambda }^{l}{\lambda }^{g}{\rm{\nabla }}{P}_{c}$	(61)
Capillary pressure Pc	${P}_{c}={P}^{g}-{P}^{l}=\sigma \,\cos \,\theta {\left(\displaystyle \frac{\varepsilon }{K}\right)}^{1/2}J\left(s\right)$	(62)
Leverett function J(s)	$J=\left\{\begin{array}{c}1.417(1-s)-2.120{(1-s)}^{2}+1.263{(1-s)}^{3}\\ 1.417s-2.120{s}^{2}+1.263{s}^{3}\end{array}\right.\begin{array}{c}{\rm{if}}\,{\theta }_{c}\lt {90}^{\circ }\\ {\rm{if}}\,{\theta }_{c}\gt {90}^{\circ }\end{array}$	(63)

Table V. Transport properties in the electrolyte.

Description	Expression
Membrane water content (λ)	$\lambda =\left\{\begin{array}{c}{\lambda }_{g}=0.043+17.81a-39.85{a}^{2}+36.0{a}^{3}\,for\,0\lt a\leqslant 1\,\\ {\lambda }_{l}=22\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\end{array}\right.$	(64)
	Water activity, $a=\displaystyle \frac{{C}_{w}^{g}{R}_{u}T}{{P}_{sat}}$ [00]	(65)
Electro-osmotic drag (EOD) coefficient of water $\left({n}_{d}\right)$	${n}_{d}=\displaystyle \frac{2.5\lambda }{22}$	(66)
Proton conductivity ($\kappa$)	$\kappa =\left(0.5139\lambda -0.326\right)\exp \left.\left[1268\left(\displaystyle \frac{1}{303}-\displaystyle \frac{1}{T}\right)\right.\right]$	(67)
Water diffusion coefficient (${D}_{w}^{mem}$)	${D}_{w}^{mem}=\left\{\begin{array}{c}2.692661843\cdot {10}^{-10}for\,\lambda \leqslant 2\\ \left\{\left.0.87\left(3-\lambda \right)+2.95\left(\lambda -2\right)\right\}\cdot {10}^{-10}\cdot {e}^{\left(7.9728-2416/T\right)}\,for\,2\lt \lambda \leqslant 3\right.\\ \left\{\,\{2.95\left(4-\lambda \right)+1.642454\left(\lambda -3\right)\right\}\cdot {10}^{-10}\cdot {e}^{\left(7.9728-2416/T\right)}for3\lt \lambda \leqslant 4\\ \left(\,(2.563-0.33\lambda +0.0264{\lambda }^{2}-0.000671{\lambda }^{3}\right)\cdot {10}^{-10}\cdot {e}^{\left(79728-2416/T\right)}for\,4\lt \lambda \leqslant {\lambda }_{a=1}^{g}\end{array}\right.$	(68)
Interfacial resistance of the water film	${{\rm{\Omega }}}_{w,int}={z}_{w}\displaystyle \frac{{\delta }_{w}}{{D}_{{O}_{2},w}}$	(69)

Figure 1a shows the computational domain of the single-cell PEMFC considered herein. Moreover, Table VI presents the details of the operating conditions and thermal properties. The anode and cathode inlet velocities are given by Eqs. 70 and 71 and are determined using the stoichiometric ratios (${\xi }_{a}\,\mathrm{and}\,{\xi }_{c},$ respectively). In terms of the mass flow, the no-slip and impermeability conditions have been considered at all the external surfaces except the inlet and outlet regions of the anode and cathode gas channels. While considering the thermal transport in the computational domain, an isothermal boundary condition is applied to the sidewalls of the anode and cathode and an adiabatic boundary condition is applied to the top and bottom surfaces (Eqs. 72 and 73). The PEMFC operation can be performed under a galvanostatic or potentiostatic mode. Hence, either a constant voltage or current density is applied at the outer sidewall of the cathode (Eq. 75) while the electric potential ${\phi }_{s}$ is fixed to zero; this is shown in Eq. 74. To determine the appropriate number of meshes, a grid-independent study is performed, and the results are presented in Fig. 1b. Considering the analysis accuracy and time required, the number of mesh elements is determined to be about 200,000. The PEMFC model herein is numerically implemented into a commercially available CFD program ANSYS Fluent ver. 19 (ANSYS, Inc., US) by employing user-defined functions. For the equation residuals, the convergence criterion is set to 10⁻⁸.

Zoom In Zoom Out Reset image size

Table VI. Boundary conditions.

Description	Expression
Anode inlet velocity	${u}_{in,a}=\displaystyle \frac{{\xi }_{a}\displaystyle \frac{I}{2F}{A}_{mem}}{{C}_{{H}_{2}}{A}_{a,ch}}$	(70)
Cathode inlet velocity	${u}_{in,c}=\displaystyle \frac{{\xi }_{c}\displaystyle \frac{I}{4F}{A}_{mem}}{{C}_{{O}_{2}}{A}_{c,ch}}$	(71)
Constant temperature on side walls	${T}_{side}=60^\circ C$	(72)
Heat flux on the top and bottom surfaces	${k}^{eff}{\rm{\nabla }}T=0$	(73)
Electric potential on anode outer BP surface	${\phi }_{s}=0$	(74)
Electric potential on cathode outer BP surface	${\phi }_{s}={V}_{cell}\,or\,\sigma {\rm{\nabla }}{\phi }_{s}=I$	(75)

RSA is an imperative technique based on mathematical and statistical methods and is used to build global approximations of an objective via functional evaluations at various points in the design space. In this study, RSA is used to denote a polynomial approximation wherein the training data are fitted by least square regression.

Generally, when a process with an observed response $y,$ input variables (${x}_{1},\ldots \ldots .,{x}_{d}$), and number of design input variables d are considered, the response of the process in terms of the variables is given as follows. ³³

$$\begin{eqnarray}&&y=f\left({x}_{1},\ldots \ldots .,{x}_{d}\right)+\varepsilon \end{eqnarray} \tag{ 1 }$$

where $f$ is the actual response function and $\varepsilon$ is the source of errors that are not accounted for in $f.$ Since the form of the actual response function is unknown, RSA is used to develop an approximate function for the actual function $f.$ The approximation is usually performed using a polynomial function and it aims to find a predicted response within the region of interest. The polynomial that is used to represent the predicted response $\hat{y}$ is usually a first- or second-order equation represented in Eqs. 2 and 3, respectively. ³³

$$\begin{eqnarray}&&{\hat{y}}_{i}={\beta }_{0}+\displaystyle \sum _{d}^{k=1}{\beta }_{k}{x}_{k}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\hat{y}}_{i}={\beta }_{0}+\displaystyle \sum _{k=1}^{d}{\beta }_{k}{x}_{k}+\displaystyle \sum _{k=1}^{d}{\beta }_{kk}{x}_{k}^{2}+\,\displaystyle \sum _{d}^{j=1}\displaystyle \sum _{k=1,k\lt j}^{d-1}{\beta }_{kj}{x}_{k}{x}_{j}\end{eqnarray} \tag{ 3 }$$

where ${x}_{k\,}\mathrm{and}\,{x}_{j}$ represent the set of input variables and the $\beta \,$parameters are called the weight coefficients.

Furthermore, we demonstrate the response surface generation using a second-order function approximation. The advantage of using the second-order equation (Eq. 3) is that it can smooth out the various levels of errors in the data and can capture the global trend in the variations; this makes it more robust and well suited for engineering design optimization. Based on the training data, the weight coefficients are determined by minimizing the mean square error (MSE) between the observed data and the prediction values, which is given as

$$\begin{eqnarray}&&MSE=\,\displaystyle \sum _{N}^{i=1}{\displaystyle \frac{\left({\hat{y}}_{i}-{y}_{i}\right)}{N}}^{2}\end{eqnarray} \tag{ 4 }$$

where ${y}_{i}$ denotes the observed data at the set of input variables, while ${\hat{y}}_{i}$ denotes the prediction values obtained using Eq. 3 within the range of N training data. To simply solve Eq 4, a matrix system is constructed (Eqs. 5 and 6).

$$\begin{eqnarray}&&{\boldsymbol{Y}}={{\boldsymbol{X}}}_{{\rm{RSA}}}{\boldsymbol{B}}+\varepsilon \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}{{\boldsymbol{X}}}_{{\rm{RSA}}}=\left[\begin{array}{cccc}1 & {x}_{11} & \cdots & {x}_{1p-1}\\ 1 & {x}_{21} & \cdots & {x}_{2p-1}\\ \vdots & \vdots & & \vdots \\ 1 & {x}_{N1} & \cdots & {x}_{Np-1}\end{array}\right],\,{\boldsymbol{Y}}=\left[\begin{array}{c}\begin{array}{c}{y}_{1}\\ {y}_{2}\\ \vdots \end{array}\\ {y}_{N}\end{array}\right],{\boldsymbol{B}}=\left[\begin{array}{c}\begin{array}{c}{\beta }_{0}\\ {\beta }_{1}\\ \vdots \\ {\beta }_{p-1}\end{array}\end{array}\right],{\boldsymbol{\varepsilon }}=\left[\begin{array}{c}\begin{array}{c}{\varepsilon }_{1}\\ {\varepsilon }_{2}\\ \vdots \end{array}\\ {\varepsilon }_{N}\end{array}\right]\end{eqnarray} \tag{ 6 }$$

where p represents the number of weight coefficients in each second-order equation that can be expressed in terms of the number of design variables d, $p=\tfrac{\left(d+2\right)\left(d+1\right)}{2}.$ Since N training data points yield N number of responses (${\boldsymbol{Y}}$) and errors (${\boldsymbol{\varepsilon }}$), the dimension of ${{\boldsymbol{X}}}_{{\rm{RSA}}}$ should be $N\times p.$ Equation 4 is further extended and expressed in the form of a matrix system as

$$\begin{eqnarray}&&\begin{array}{l}L=\displaystyle \sum _{n}^{i=1}{\varepsilon }_{i}^{2}={{\boldsymbol{\varepsilon }}}^{T}{\boldsymbol{\varepsilon }}={\left({\boldsymbol{Y}}-{{\boldsymbol{X}}}_{{\rm{RSA}}}{\boldsymbol{B}}\right)}^{T}\left({\boldsymbol{Y}}-{{\boldsymbol{X}}}_{{\rm{RSA}}}{\boldsymbol{B}}\right)\\ \,=\,{{\boldsymbol{Y}}}^{T}{\boldsymbol{Y}}-{{\boldsymbol{B}}}^{T}{{{\boldsymbol{X}}}_{{\rm{RSA}}}}^{T}{\boldsymbol{Y}}-{{\boldsymbol{Y}}}^{T}{{\boldsymbol{X}}}_{{\rm{RSA}}}{\boldsymbol{B}}-{{\boldsymbol{B}}}^{T}{{{\boldsymbol{X}}}_{{\rm{RSA}}}}^{T}{{\boldsymbol{X}}}_{{\rm{RSA}}}{\boldsymbol{B}}\end{array}\end{eqnarray} \tag{ 7 }$$

where L denotes the least square function. To determine the weight coefficient vector, ${\boldsymbol{B}},$ L is minimized by setting its derivative, $\displaystyle \frac{\partial L}{\partial {\boldsymbol{B}}},$ as zero.

$$\begin{eqnarray}&&\displaystyle \frac{\partial L}{\partial {\boldsymbol{B}}}=-2{{{\boldsymbol{X}}}_{{\rm{RSA}}}}^{T}{\boldsymbol{Y}}+2{{{\boldsymbol{X}}}_{{\rm{RSA}}}}^{T}{{\boldsymbol{X}}}_{{\rm{RSA}}}{\boldsymbol{B}}=0\end{eqnarray} \tag{ 8 }$$

By simplifying Eq. 8, the expression of B is given by

$$\begin{eqnarray}&&{\boldsymbol{B}}={\left[{{{\boldsymbol{X}}}_{{\rm{RSA}}}}^{T}{{\boldsymbol{X}}}_{{\rm{RSA}}}\right]}^{-1}{{{\boldsymbol{X}}}_{{\rm{RSA}}}}^{T}{\boldsymbol{Y}}\end{eqnarray} \tag{ 9 }$$

Neural networks are surrogate models that mimic the neuron structure. Neurons, the basic building blocks of nerves, are composed of three parts: receiving, processing, and judgement parts. As shown in Fig. 3, these are named input, hidden, and output layers, respectively. After the input layer receives information, the hidden layer processes the information, and then, the output layer makes a prediction. The function that processes information in the hidden layer is called the activation function. The activation function type varies depending on the neural network used. In the RBNN method, the Gaussian-radial function is used as the activation function, i.e., expressed as a function of input vector ${\boldsymbol{X}}$ comprising design variables (${x}_{1},{x}_{2},\ldots \mathrm{..}{x}_{d}$) as follows:

$$\begin{eqnarray}&&{h}_{i}\left({\boldsymbol{X}}\right)=\exp \left(-\displaystyle \frac{{\boldsymbol{X}}-{{{\boldsymbol{c}}}_{i}}^{2}}{2{{\sigma }_{i}}^{2}}\right)\end{eqnarray} \tag{ 10 }$$

where ${h}_{i}$ is the ${i}^{th}$ activation function. Moreover, ${{\boldsymbol{c}}}_{i}$ and ${\sigma }_{i}$ are the center and variance of input training data, respectively. ∥∥ denotes the Euclidean norm.

Zoom In Zoom Out Reset image size

The output layer that is the predicted response of RBNN is a linear combination of the activation functions and the corresponding weights, and it is given as

$$\begin{eqnarray}&&\hat{y}=\displaystyle \sum _{i=1}^{m}{\beta }_{i}{h}_{i}\left({\boldsymbol{X}}\right)\end{eqnarray} \tag{ 11 }$$

where $\hat{y}$ denotes the predicted value for an input vector ${\boldsymbol{X}}$ and ${\beta }_{i}$ is the $i$ ^th weight coefficient. Further, m is the number of activation functions and weight coefficients that should be less than the number of training data points ($m\leqslant N$). Since the accuracy of the RBNN model depends on the ${\beta }_{i}$ values, the learning process of the model denotes the calculation process of the best ${\beta }_{i}$ values. Based on the training data in the given domain space ${\boldsymbol{X}},$ ${\beta }_{i}$ is determined by minimizing the sum-squared error ($SSE$) between the observed data and the prediction values:

$$\begin{eqnarray}&&SSE=\,\displaystyle \sum _{N}^{i=1}{\left({\hat{y}}_{i}-{y}_{i}\right)}^{2}\end{eqnarray} \tag{ 12 }$$

where ${y}_{i}$ represents the observed response at individual ${\boldsymbol{X}},$ while ${\hat{y}}_{i}$ is a prediction value shown in Eq. 11. The weight coefficients ${\beta }_{i}$ are estimated using the least square solution of Eq. 12 in the matrix form of B :

$$\begin{eqnarray}&&{\boldsymbol{B}}={{\boldsymbol{H}}}^{T}{\boldsymbol{H}}{{\boldsymbol{H}}}^{T}{\boldsymbol{Y}}\end{eqnarray} \tag{ 13 }$$

where ${\boldsymbol{H}}$ is a design matrix comprising the individual activation functions. Here, the dimension of ${\boldsymbol{H}}$ should be $N\times m,$ where N training data points of the input layer are assigned to m activation functions of the hidden layer. ${\boldsymbol{Y}}\,$is a $N\times 1$ vector containing observed responses corresponding to the training data points, ${\boldsymbol{X}}.$

$$\begin{eqnarray}\begin{array}{l}{\boldsymbol{H}}=\left[\begin{array}{ccccc}{h}_{1}\left({{\boldsymbol{X}}}_{1}\right) & {h}_{2}\left({{\boldsymbol{X}}}_{1}\right) & {h}_{3}\left({{\boldsymbol{X}}}_{1}\right) & \cdots & {h}_{m}\left({{\boldsymbol{X}}}_{1}\right)\\ {h}_{1}\left({{\boldsymbol{X}}}_{2}\right) & {h}_{2}\left({{\boldsymbol{X}}}_{2}\right) & {h}_{3}\left({{\boldsymbol{X}}}_{2}\right) & \cdots & {h}_{m}\left({{\boldsymbol{X}}}_{2}\right)\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ {h}_{1}\left({{\boldsymbol{X}}}_{N}\right) & {h}_{2}\left({{\boldsymbol{X}}}_{N}\right) & {h}_{3}\left({{\boldsymbol{X}}}_{N}\right) & \cdots & {h}_{m}\left({{\boldsymbol{X}}}_{N}\right)\end{array}\right],\\ {\boldsymbol{B}}=\left[\begin{array}{c}\begin{array}{c}{\beta }_{1}\\ {\beta }_{2}\\ \vdots \end{array}\\ {\beta }_{m}\end{array}\right],\,{\boldsymbol{Y}}=\left[\begin{array}{c}\begin{array}{c}{y}_{1}\\ {y}_{2}\\ \vdots \end{array}\\ {y}_{N}\end{array}\right]\end{array}\end{eqnarray} \tag{ 14 }$$

KRG ³⁴ is a well-known surrogate model used to predict an accurate approximation in the design space. It can be used to predict both linear and nonlinear trends in the design space. The model herein is constructed using Eqs. 15–21.

The ith input vector of d-dimension is represented as ${{\boldsymbol{X}}}_{i}:$

$$\begin{eqnarray}&&{{\boldsymbol{X}}}_{i}={\left[{x}_{1},{x}_{2},\ldots ,{x}_{d}\right]}^{T}\,\end{eqnarray} \tag{ 15 }$$

As the observed response at ${{\boldsymbol{X}}}_{i}$ is taken as ${y}_{i},$ the vector for all N observed responses is

$$\begin{eqnarray}&&{\boldsymbol{Y}}={\left[{y}_{1},\,{y}_{2},\,\ldots ,\,{y}_{N}\,\right]}^{T}\end{eqnarray} \tag{ 16 }$$

Suppose that the predicted response from the KRG model at ${{\boldsymbol{X}}}_{i}$ is $\hat{y}\left({{\boldsymbol{X}}}_{i}\right),$ then the distance between any two chosen design points ${{\boldsymbol{X}}}_{i}$ and ${{\boldsymbol{X}}}_{l}$ is correlated with the prediction values. i.e., the difference between the predicted values reduces with the interval between the two points. The relationship between the correlated values is

$$\begin{eqnarray}&&cor\left[\hat{y}\left({{\boldsymbol{X}}}_{i}\right),\hat{y}\left({{\boldsymbol{X}}}_{l}\right)\right]=\exp \left(-\displaystyle \sum _{{\rm{j}}=1}^{d}{\beta }_{{\rm{j}}}{x}_{ij}-{{x}_{lj}}^{2}\right)\end{eqnarray} \tag{ 17 }$$

where ${x}_{ij}$ and ${x}_{lj}$ represent the $j$ ^th components of ${{\boldsymbol{X}}}_{i}$ and ${{\boldsymbol{X}}}_{l},$ respectively, and ${\beta }_{j}$ is their regression parameter. The correlation values at all the observed points can be written in terms of a $N\times N$ matrix ${\boldsymbol{\Psi }},$ i.e.,

$$\begin{eqnarray}{\rm{\Psi }}=\left[\begin{array}{ccc}cor\left[\hat{y}\left({X}_{1}\right),\hat{y}\left({X}_{1}\right)\right] & \cdots & cor\left[\hat{y}\left({X}_{1}\right),\hat{y}\left({X}_{N}\right)\right]\\ \vdots & \ddots & \vdots \\ cor\left[\hat{y}\left({X}_{N}\right),\hat{y}\left({X}_{1}\right)\right] & \cdots & cor\left[\hat{y}\left({X}_{N}\right),\hat{y}\left({X}_{N}\right)\right]\end{array}\right]\end{eqnarray} \tag{ 18 }$$

The KRG predictor $\hat{y}\left({\boldsymbol{X}}\right)$ is

$$\begin{eqnarray}&&\hat{y}\left({\boldsymbol{X}}\right)=\hat{\mu }+{{\boldsymbol{\psi }}}^{T}{{\boldsymbol{\Psi }}}^{-1}\left({\boldsymbol{Y}}-1\hat{\mu }\,\,\right)\end{eqnarray} \tag{ 19 }$$

where $1$ is a singular matrix of size $N\times 1$ and $\hat{\mu }$ is the mean value of the training dataset distribution that is calculated through maximum likelihood estimation:

$$\begin{eqnarray}&&\hat{\mu }=\displaystyle \frac{{1}^{{\boldsymbol{T}}}{{\boldsymbol{\Psi }}}^{-1}{\boldsymbol{Y}}}{{1}^{{\boldsymbol{T}}}{{\boldsymbol{\Psi }}}^{-1}1}\end{eqnarray} \tag{ 20 }$$

Here, ${\boldsymbol{\psi }}$ is the vector comprising correlation values between a given prediction point and all training data. The structure of ${\boldsymbol{\psi }}$ is expressed as

$$\begin{eqnarray}&&{\boldsymbol{\psi }}=\left[\begin{array}{c}cor\left[Y\left({\boldsymbol{X}}\right),Y\left({{\boldsymbol{X}}}_{1}\right)\right]\\ \begin{array}{c}\vdots \\ cor\left[Y\left({\boldsymbol{X}}\right),Y\left({{\boldsymbol{X}}}_{N}\right)\right]\end{array}\end{array}\right]\end{eqnarray} \tag{ 21 }$$

PSO is one of the most widely used optimization technique developed by Kennedy and Eberhart, ³⁵ and it is based on the communal behavior of particles in a population, and is also known as a swarm. The algorithm is a metaheuristic that is used to optimize nonlinear continuous functions and is based on an evolutionary search to obtain the global best. Figure 4 describes the functioning of the algorithm. In PSO, S represents the population of N number of particles, and the particles represent points that are uniformly distributed within the dimensional space $d,$ where each particle can flow through $d$ based on their individual and swarm behaviors.

$$\begin{eqnarray}&&{\boldsymbol{S}}={{\boldsymbol{X}}}_{1},{{\boldsymbol{X}}}_{2},{{\boldsymbol{X}}}_{3}\ldots {{\boldsymbol{X}}}_{N}\end{eqnarray} \tag{ 22 }$$

Zoom In Zoom Out Reset image size

Accordingly, the position, velocity, and fitness values of the $k$ ^th particles ($k=1,2,\cdots ,N$) can be described in the form of a vector, as given in Eqs. 23–25

Position vector,

$$\begin{eqnarray}&&{{\boldsymbol{X}}}_{k}=\left[{x}_{k1},{x}_{k2},\ldots ,{x}_{kd}\right]\end{eqnarray} \tag{ 23 }$$

Velocity vector,

$$\begin{eqnarray}&&{{\boldsymbol{V}}}_{k}=\left[{v}_{k1},{v}_{k2},\ldots ,{v}_{kd}\right]\end{eqnarray} \tag{ 24 }$$

Fitness values,

$$\begin{eqnarray}&&{f}_{k}=f\left({{\boldsymbol{X}}}_{k}\right)\end{eqnarray} \tag{ 25 }$$

In Eq. 25, ${f}_{k}$ is the objective function at each ${{\boldsymbol{X}}}_{k}.$ While searching for the optimal solution for the given problem, the particles update their positions and velocities according to each iteration as follows:

Particle velocity update,

$$\begin{eqnarray}&&{{\boldsymbol{V}}}_{k}\left(t+1\right)={{\boldsymbol{V}}}_{k}\left(t\right)+{C}_{1}{{\boldsymbol{R}}}_{1}\left({{\boldsymbol{P}}}_{k}-{{\boldsymbol{X}}}_{k}\left(t\right)\right)+{C}_{2}{{\boldsymbol{R}}}_{2}\left({{\boldsymbol{G}}}_{k}-{{\boldsymbol{X}}}_{k}\left(t\right)\right)\end{eqnarray} \tag{ 26 }$$

Particle position update,

$$\begin{eqnarray}&&{{\boldsymbol{X}}}_{k}\left(t+1\right)={{\boldsymbol{X}}}_{k}\left(t\right)+{{\boldsymbol{V}}}_{k}\left(t+1\right)\end{eqnarray} \tag{ 27 }$$

where $t\,$and $\left(t+1\right)$ indicate two consecutive iterations of the algorithm. ${{\boldsymbol{V}}}_{k}\left(t\right)$ and ${{\boldsymbol{X}}}_{k}\left(t\right)$ are the velocity and position vectors of the particle k at the ${t}^{{\rm{th}}}$ iteration, respectively. Further, ${{\boldsymbol{P}}}_{k}$ denotes the "personal best" of the ${k}^{{\rm{th}}}$ particle, i.e., the position of the best solution obtained during the entire iteration. ${{\boldsymbol{G}}}_{k}$ denotes the "global best," representing the position of the overall best solution obtained in the population group (S).${C}_{1}\mathrm{and}\,{C}_{2}$ are the acceleration factors and ${{\boldsymbol{R}}}_{1}$ and ${{\boldsymbol{R}}}_{2}$ are random numbers between 0 and 1. As shown in Eq. 26, ${{\boldsymbol{V}}}_{k}\left(t+1\right)$ is depends on three terms: the inertia term, ${{\boldsymbol{V}}}_{k}\left(t\right);$ cognitive term, ${C}_{1}{{\boldsymbol{R}}}_{1}\left({{\boldsymbol{P}}}_{k}-{{\boldsymbol{X}}}_{k}\left(t\right)\right);$ and social term, ${C}_{2}{{\boldsymbol{R}}}_{2}\left({{\boldsymbol{G}}}_{k}-{{\boldsymbol{X}}}_{k}\left(t\right)\right).$ While the inertia term prevents the particle from making drastic changes in direction by keeping a trail of the previous flow direction, the cognitive term helps the particle return to its previous best solution and the social term directs the particles to move toward the best solution in the entire population. While updating ${{\boldsymbol{V}}}_{k},$ the acceleration factors ${C}_{1}\,{\rm{and}}\,{C}_{2},$ which modulate the magnitudes of the steps in the direction of its personal best and global best, respectively, are usually chosen within the range of 0–4. ³⁶

Figure 4 briefly shows the working procedures of the PSO algorithm that comprises the following basics steps:

• 1.  
  Initialization step.
• 2.  
  Iteration step.
• 3.  
  Judging the best solution.

In the initialization step, the size of S is defined (k = 1, 2, ..., N) and the initial positions of the particles are considered as ${{\boldsymbol{X}}}_{k}\left(0\right).$ At every ${{\boldsymbol{X}}}_{k}\left(0\right),$ ${{\boldsymbol{P}}}_{k}$ is initialized to ${{\boldsymbol{X}}}_{k}\left(0\right).$ The fitness values of all particles at ${{\boldsymbol{X}}}_{k}\left(0\right)$ are evaluated as $f\left({{\boldsymbol{X}}}_{k}\left(0\right)\right),$ and the particle with the highest fitness value is assigned as the global best, ${{\boldsymbol{G}}}_{k}.$ In the second step, the values of ${{\boldsymbol{V}}}_{k}\left(t\right)$ and ${{\boldsymbol{X}}}_{k}\left(t\right)$ are updated to ${{\boldsymbol{V}}}_{k}\left(t+1\right)$ and ${{\boldsymbol{X}}}_{k}\left(t+1\right)$ according to Eqs. 26 and 27, respectively, and the fitness values $f\left({{\boldsymbol{X}}}_{k}\left(t+1\right)\right)\,$ are evaluated. If the updated fitness values $f\left({{\boldsymbol{X}}}_{k}\left(t+1\right)\right)$ are greater than the previous fitness values $f\left({{\boldsymbol{X}}}_{k}\left(t+1\right)\right)\geqslant f\left({{\boldsymbol{P}}}_{{\boldsymbol{k}}}\right),$ the ${{\boldsymbol{X}}}_{k}\left(t+1\right)$ values are assigned to new ${{\boldsymbol{P}}}_{k}$ (${{\boldsymbol{P}}}_{k}={{\boldsymbol{X}}}_{k}\left(t+1\right)$). Similarly, the global best ${{\boldsymbol{G}}}_{k}$ is updated by evaluating whether $f\left({{\boldsymbol{X}}}_{k}\left(t+1\right)\right)\,$ is greater than the desired fitness $f\left({{\boldsymbol{G}}}_{{\boldsymbol{k}}}\right).$ The algorithm terminates after all the particles reach the same global best value ${{\boldsymbol{G}}}_{k}$ or the maximum iteration has been reached.

Similar to PSO, GA employs metaheuristic and iterative processes for optimization inspired by natural selection. The optimization process mimics biological processes wherein genes are selected, mated, and mutated. ³⁷ GA comprises the following four steps (Fig. 5): generation of parent chromosomes, mating, selection, and mutation.

Zoom In Zoom Out Reset image size

In the first step, the process begins by randomly generating 1^st generation parent chromosomes to form population P, comprising N number of design vectors called chromosomes, which is denoted by ${{\boldsymbol{X}}}_{k}$ ($k=1,\,2,\,\,\ldots ,\,N$). Each chromosome is a d-dimensional vector comprising d number of genes represented by ${x}_{kd};$ random values between 0 and 1 are assigned to it. ${\boldsymbol{S}}$ and ${{\boldsymbol{X}}}_{k}$ are expressed as follows:

$$\begin{eqnarray}&&{\boldsymbol{S}}={{\boldsymbol{X}}}_{1},{{\boldsymbol{X}}}_{2},\cdots ,{{\boldsymbol{X}}}_{N}\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&{{\boldsymbol{X}}}_{k}=\left[{x}_{k1},{x}_{k2},\cdots ,{x}_{kd}\right]\end{eqnarray} \tag{ 29 }$$

In GA, the individual chromosomes in population P are evaluated at the end of each generation, where the fitness value (i.e., objective function value ${f}_{k}\,$) is represented as

$$\begin{eqnarray}&&{f}_{k}=f\left({{\boldsymbol{X}}}_{k}\right)\end{eqnarray} \tag{ 30 }$$

Once the ${1}^{st}$ generation of N chromosomes have been released, the second step, i.e., mating, begins to produce new offsprings N. Two offsprings ${{\boldsymbol{X}}}_{w}=\left[{x}_{w1},{x}_{w2},\cdots ,{x}_{wd}\right]$ and ${{\boldsymbol{X}}}_{z}=\left[{x}_{z1},{x}_{z2},\cdots ,{x}_{zd}\right]$ are afforded from two random parents in P, namely, ${{\boldsymbol{X}}}_{f}=\left[{x}_{f1},{x}_{f2},\cdots ,{x}_{fd}\right]$ and ${{\boldsymbol{X}}}_{m}=\left[{x}_{m1},{x}_{m2},\cdots ,{x}_{md}\right].$ The process of generating offsprings is controlled by the mating rate q that varies between 0 and 1:

Offspring 1,

$$\begin{eqnarray}&&{x}_{zi}=\left(1-q\right){x}_{fi}+q{x}_{mi}\end{eqnarray} \tag{ 31 }$$

Offspring 2,

$$\begin{eqnarray}&&{x}_{wi}=\left(1-q\right){x}_{mi}+q{x}_{fi}\end{eqnarray} \tag{ 32 }$$

where ${x}_{zi}$ and ${x}_{wi}$ represent the ${i}^{th}$ genes ($i=1,\,2,\,\,\ldots ,\,\,d$) of the offsprings 1 and 2, respectively. At the end of the mating process, a total population of 2N is obtained.

In the third step, based on the fitness values of individual chromosomes, N chromosomes with high fitness values are selected among 2N chromosomes. In the final step, a mutation that induces changes in genes of random chromosomes is performed with a very small predetermined probability of 0.1 to secure genetic diversity. ³⁸ Then, the newly selected N chromosomes become parent chromosomes in the next iteration. Thereafter, the iterative process continues from step 2 until the stopping criteria are met, i.e., a satisfactory fitness level has been obtained for the population or a maximum number of allowed iterations has been reached.

Once the surrogate models are developed using the generated training data, the accuracy of these models is tested at 16 test points, i.e. selected from 85 sample points. Numerical error estimators (${\text{adjusted}\,R}^{2}$ and RMSE) are used to evaluate the accuracy of the surrogate models, i.e., using Eqs. 33 and 34, respectively.

$$\begin{eqnarray}&&{\mathrm{adjusted}\,R}^{2}=1-\displaystyle \frac{r-1}{r-d-1}\left(1-\left(1-\displaystyle \frac{{\sum }_{i=1}^{r}{\left({y}_{i}-{\hat{y}}_{i}\right)}^{2}}{{\sum }_{i=1}^{r}{\left({y}_{i}-{\bar{y}}_{i}\right)}^{2}}\right)\right)\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&{\rm{RMSE}}=\sqrt{\displaystyle \frac{1}{r}\displaystyle \sum _{r}^{i=1}{\left({\hat{y}}_{i}-{y}_{i}\right)}^{2}}\end{eqnarray} \tag{ 34 }$$

where r is the number of test data points. Moreover, ${y}_{i}$ and ${\hat{y}}_{i}$ denote the observed responses that can be obtained using CFD simulations and the predicted values from the surrogate models, respectively; ${\bar{y}}_{i}$ denotes the mean value of the observed data at r test points and $d$ is the number of design parameters. The obtained ${\text{adjusted}\,R}^{2}$ and RMSE values demonstrate the fit of the model with the data considered herein. The ${\text{adjusted}\,R}^{2}$ values range from zero to one, where zero indicates a poor model accuracy and a value close to one indicates perfect predictions. For RMSE, smaller values, i.e., close to zero, indicate better model accuracy and higher values, i.e., close to one, indicate poor accuracy.