Reliability-Based Design Optimization of the PEMFC Flow Field with Consideration of Statistical Uncertainty of Design Variables

Abstract

Recently, with the fourth industrial revolution, the research cases that search for optimal design points based on neural networks or machine learning have rapidly increased. In addition, research on optimization is continuously reported in the field of fuel cell research using hydrogen as fuel. However, in the case of optimization research, it often requires a large amount of training data, which means that it is more suitable for numerical research such as CFD simulation rather than time-consuming research such as actual experiments. As is well known, the design range of fuel cell flow channels is extremely small, ranging from hundreds of microns to several millimeters, which means the small tolerance could cause fatal performance loss. In this study, the general optimization study was further improved in terms of reliability by considering stochastic tolerances that may occur in actual industry. The optimization problem was defined to maximize stack power, which is employed as objective function, under the constraints such as pressure drop and current density standard deviation; the performance of the optimal point through general optimization was about 3.252 kW/L. In the reliability-based optimization problem, the boundary condition for tolerance was set to 0.1 mm and tolerance was assumed to occur along a normal distribution. The optimal point to secure 99% reliability for the given constraints was 2.918 kW/L, showing significantly lower performance than the general optimal point.

1. Introduction

In the context of emphasizing the environmental friendliness of energy, hydrogen energy plays a crucial role [1]. While there are various ways to harness it, one example is fuel cells that utilize precisely designed electrochemical reactions to directly convert the chemical energy of hydrogen into electrical energy [2]. Depending on the type, fuel cells can operate in both high-temperature and low-temperature environments, and possess rapid responsiveness and portability [3,4].

Fuel cells are composed of a bipolar plate with engraved flow channels, a gas diffusion layer for uniform reactant diffusion, and a membrane electrolyte assembly where electrochemical reactions and ion transfer occur. Material transfer occurs through the flow channels engraved on the bipolar plate and the type and shape of these channels are crucial for determining the concentration overpotential which is one of the three types of overpotentials (activation overpotential, ohmic overpotential, and concentration overpotential) determining the performance of a fuel cell. Particularly, in high-power fuel cells, these channels play an extremely decisive role in the performance. In addition, there is a wide range of parameters in fuel cells, leading to active research in this area [5,6,7,8]. Especially in recent times, research on commercialization has begun in earnest, providing ample opportunities for optimization in various areas of fuel cell design, including materials, operating conditions, and shapes. There are various optimization methods and one intuitively accessible approach is empirical optimization, where optimization is carried out based on understanding the influence of parameters on performance [9,10]. Although this approach is easy to implement, it comes with high human and material costs, and the derived values rely on experience, making it challenging to ensure global optimal values or even confidence in local optimal values across the entire variable space. To address the drawbacks of such empirical methods, research has been conducted and applied in utilizing computer engineering-based optimization algorithms [11]. Particle swarm optimization (PSO) and genetic algorithms (GA) serve as examples of optimization algorithms implemented computationally, utilizing behavioral and evolutionary characteristics inspired by biology. They have been widely employed in optimizing various engineering problems [12,13].

Optimization techniques, through various methodologies, often require a substantial amount of responses to track the minimum or maximum values of the objective function. Therefore, to enhance the computational efficiency of optimization techniques, there is active utilization of regression models through artificial intelligence (AI) [14,15]. For fuel cells with various parameters, it is possible to model them diversely depending on the phenomena to be described [16]. Particularly, understanding the impact of numerical values and shapes at the cell and stack levels on performance involves a complex combination of physical phenomena, making it challenging to intuitively grasp. Consequently, there have been diverse attempts to optimize this through simulation and AI in recent years [13,17,18].

On the other hand, in some cases, there may be a need for more constraints rather than optimizing the entire variable space. For example, the operating temperature of a fuel cell is a crucial factor determining its performance. The electrochemical reaction rate increases with temperature, so under conditions with adequate moisture supply, the performance of a fuel cell increases with temperature. However, exceeding the temperature at which water vaporizes leads to a sharp decline in performance due to the mechanical deformation of the Nafion membrane [19]. In other words, due to design constraints, there is a restriction on the operating temperature, which ultimately influences the objective function, namely, the performance. Engineering objectives can lead to various constraints on variables, and in such cases, constrained optimization can be employed to address these issues. Linear programming (LP) and sequential quadratic programming (SQP) are representative methodologies for constrained optimization [20,21]. Constraints can be defined in the form of equalities or inequalities and can take the shape of linear or nonlinear functions. In constrained optimization, the goal is to obtain the optimal solution under conditions that satisfy these nonlinear (or linear) inequality (or equality) constraints.

In most engineering problems, the ultimate goal of optimization is to find the design point that represents the maximum or minimum value of the objective function within the valid design space. Typically, due to the nature of optimal design which requires a vast amount of data, computer simulation, as mentioned earlier, is actively utilized. However, computer simulation using deterministic variables focuses on identifying physical phenomena and behaviors, and indeed, computer simulation and traditional deterministic optimization may have the probabilistic disadvantage of not reflecting the uncertainties that can occur in actual industrial processes. Reliability-based design optimization (RBDO), which has been developed by K.K. Choi [22] and others [23,24], is designed to narrow the gap between such real-world industrial uncertainties and deterministic optimization based on probabilistic theory. It has been observed that uncertainties such as tolerances that can occur in actual industrial processes can be reflected by considering the probability distribution of each design variable, and through this, it is possible to derive an optimal design solution that guarantees reliability [25,26,27].

In this study, we introduce research that has optimized key design variables—channel width, land width, and channel depth—of parallel flow channels in PEMFC bipolar plates, focusing not only on performance but also on various objective functions such as pressure drop and current density uniformity. First, the multi-scale, multi-phase PEMFC model was validated against experimental data from PEMFCs with parallel flow channels and used to generate comprehensive data on key design variables via simulation. A neural network-based surrogate model was then trained with the acquired database to ascertain the optimal values for the parallel flow field design. A novel aspect of this study, which differentiates it from other optimization studies, was the application of a probability distribution to the optimal design values to account for the uncertainties and tolerance effects in the bipolar plate flow path production process, thus verifying the values’ reliability. Typically, the optimal design values obtained from general multi-objective optimization can partially fall into the infeasible region when uncertainties from the production process are accounted for with a probability distribution applied to these values. This study successfully demonstrated that it was necessary to revise the optimal design values, ensuring all output data fall within the feasible region through a probability statistics-based reliability optimization process. It also demonstrated the essential role of a novel reliability-based optimal design technology, which merged a three-dimensional physical model, artificial intelligence algorithms, and probabilistic statistical theory, in deriving practical optimal design values.

2. Multi-Scale, Multi-Phase 3D PEMFC Model

A 3D, multi-scale, two-phase PEMFC model based on the multi-phase mixture (M2) model proposed by Wang and Cheng [28] was considered in the current study. Modeling the PEMFC was carried out considering the various components such as the membrane, CLs, GDLs, and BPs. The 3D PEMFC model has been previously validated against experimental polarization curves measured under various cell designs and operating conditions [29]. The model assumptions and governing equations considered in this study are similar to the ones described in our previous study. Therefore, a very brief summary of the same has been reported in the Section 2. Additionally, the boundary conditions and numerical implementations of the PEMFC model featured using commercial computational fluid dynamics (CFD) software; ANSYS Fluent v23(ANSYS Inc., 2600 Ansys Dr, Canonsburg, PA, USA) is summarized in the Section 2. Furthermore, the multi-scale, two-phase PEMFC model employed in this study is shown in Figure 1.

2.1. Model Assumptions

The following four assumptions were applied to the proposed 3D DU hydride model:

(1)

The present study was guided by a set of specific assumptions outlined below.

(2)

Ideal gas mixtures in the gas phase were considered due to the low operating pressures involved.

(3)

Flow velocity was assumed to be low, maintaining a laminar condition.

(4)

Negligible effects of gravity were taken into account.

(5)

The influence of immobile liquid saturation in porous regions was considered to be negligible.

2.2. Governing Equations and Source Terms

The 3D PEMFC model under consideration in this study was subject to five conservation principles: mass, species, charge, momentum, and thermal energy. These conservation equations were intricately linked with source terms corresponding to electrochemical reactions taking place in both the anode (hydrogen oxidation reaction, HOR) and cathode (oxygen reduction reaction, ORR). A summary of the governing equations and their associated source terms is available in

Table 1. PEM fuel cell model: Governing equations [17].

Governing Equations
Mass	$\nabla \bullet \left(\mathsf{\rho }\overrightarrow{u}\right)=0$	(1)
Momentum	$\left(\frac{1}{{\mathsf{\epsilon }}^2}\right)\nabla \bullet \left(\rho \overrightarrow{u}\overrightarrow{u}\right)=-\nabla P+\nabla \bullet \tau +S_u$	(2)
Species	Flow channels and porous media: $\nabla \bullet \left({\mathsf{\gamma }}_{\mathrm{i}}\rho m_i\overrightarrow{u}\right)=\nabla \bullet \left[{\rho }^g{D}_i^{g,eff}\nabla \left({\mathrm{m}}_{\mathrm{i}}^{\mathrm{g}}\right)\right]+\nabla \bullet \left[\left(m_i^g-m_i^l\right){\overrightarrow{j}}^l\right]+S_i$	(3)
	Water transport in membrane: $\nabla \bullet \left({\left(\frac{{\mathsf{\rho }}^{\mathrm{m}\mathrm{e}\mathrm{m}}}{EW}\right)\mathrm{D}}_{\mathrm{w}}^{\mathrm{m}\mathrm{e}\mathrm{m}}\nabla \lambda \right)M_w-\nabla \bullet \left(n_d\left(\frac{I}{F}\right)\right)M_w+\nabla \bullet (\left(\frac{{\kappa }^{mem}}{{\nu }^l}\right)\nabla P^l=0$	(4)
Charge	Proton transport: $\nabla \bullet \left({\mathsf{\kappa }}^{\mathrm{e}\mathrm{f}\mathrm{f}}\nabla {\varphi }_e\right)+S_{\varphi }=0$	(5)
	Electron transport: $\nabla \bullet \left({\mathsf{\sigma }}^{\mathrm{e}\mathrm{f}\mathrm{f}}\nabla {\varphi }_s\right)-S_{\varphi }=0$	(6)
Energy	$\nabla \bullet \left(\mathsf{\rho }\overrightarrow{u}{C}_p^{g}T\right)=\nabla \bullet \left(k^{eff}\nabla T\right)+S_T$	(7)

and

Table 2. PEM fuel cell model: source/sink terms. [17].

Description	Expression
Momentum	Porous media	${\mathrm{S}}_{\mathrm{u}}=-\frac{\mu }{K}\overrightarrow{u}$
Species	H2 in anode CL	${\mathrm{S}}_{H_2,a}=\left[-\frac{j}{2F}\right]{\mathrm{M}}_{H_2}$
	O2 in cathode CL	${\mathrm{S}}_{O_2,c}=\left[\frac{j}{4F}\right]{\mathrm{M}}_{O_2}$
	Water in anode CL	${\mathrm{S}}_{w,a}=\left[-\nabla \bullet \left(\frac{n_d}{F}\overrightarrow{I}\right)+\frac{j}{4F}\right]{\mathrm{M}}_w$
	Water in cathode CL	${\mathrm{S}}_{w,c}=\left[-\nabla \bullet \left(\frac{n_d}{F}\overrightarrow{I}\right)-\frac{j}{2F}\right]{\mathrm{M}}_w$
Energy	In anode CL	${\mathrm{S}}_{\mathrm{T},\mathrm{a}}=j\bullet \eta +\frac{I^2}{k^{eff}}$
	In cathode CL	${\mathrm{S}}_{\mathrm{T},\mathrm{c}}=j\left(\eta +T\frac{d{U}_0}{dT}\right)+\frac{I^2}{k^{eff}}$
	In membrane	${\mathrm{S}}_{\mathrm{T}}=\frac{I^2}{k^{eff}}$
Charge	In CLs:	${\mathrm{S}}_{\varphi }=j$
Electrochemical reactions             ${\sum }_k{s}_i{M}_i^z=n{e}^{-},\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\left\{\begin{array}{c}{M}_i=\mathrm{c}\mathrm{h}\mathrm{e}\mathrm{m}\mathrm{i}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{u}\mathrm{l}\mathrm{a}\mathrm{o}\mathrm{f}\mathrm{s}\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{e}\mathrm{s}\mathrm{i}\\ s_i=\mathrm{s}\mathrm{t}\mathrm{o}\mathrm{i}\mathrm{c}\mathrm{h}\mathrm{i}\mathrm{o}\mathrm{m}\mathrm{e}\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{e}\mathrm{f}\mathrm{f}\mathrm{i}\mathrm{c}\mathrm{i}\mathrm{e}\mathrm{n}\mathrm{t}\\ n=\mathrm{n}\mathrm{u}\mathrm{m}\mathrm{b}\mathrm{e}\mathrm{r}\mathrm{o}\mathrm{f}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{r}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{r}\mathrm{e}\mathrm{d}\end{array}\right.$ HOR on the anode side: ${\mathrm{H}}_2-2{\mathrm{H}}^{+}=2{\mathrm{e}}^{-}$ ORR on the cathode side: $2{\mathrm{H}}_2\mathrm{O}-{\mathrm{O}}_2-4{\mathrm{H}}^{+}=4{\mathrm{e}}^{-}$
Transfer current density, $\left[\mathrm{A}/{\mathrm{m}}^3\right]$	HOR in anode CL: $j=\left(1-s\right)a{i}_{0,a}^{ref}{\left(\frac{C_{H_2}}{C_{H_2,ref}}\right)}^{\frac{1}{2}}\mathit{exp}\left[-\frac{E_a}{R}\left(\frac{1}{T}-\frac{1}{353.15}\right)\right]\left(\mathit{exp}\left(\frac{{\alpha }_{a}F}{RT}\eta \right)-\mathit{exp}\left(-\frac{{\alpha }_{c}F}{RT}\eta \right)\right)$		(8)
	ORR in cathode CL: $j=-\frac{3L_{Pt}}{{r_{Pt}{\bullet \rho }_{Pt}\bullet \delta }_{CL}}{i}_{0,c}^{ref}{\left(\frac{C_{O_2}^{Pt}}{C_{O_2,ref}}\right)}^{3/4}exp\left(-\frac{E_c}{R}\left(\frac{1}{T}-\frac{1}{353.15}\right)\right)\exp \left(-\frac{{\alpha }_c}{R_{u}T}F\eta \right)$		(9)
Overpotential	$\eta ={\varphi }_s-{\varphi }_e-U$ where $U=U_0-\frac{RT}{nF}\ln \frac{C_{O_2}}{C_{O_{2ref}}}$		(10)

, respectively. Additionally,

Table 3. Kinetic, physiochemical properties and transport of the fuel cell [17].

Description	Value/Expression	Ref.
Activation energy of anode, $E_a$	10.0 kJ/mol	[24]
Activation energy of cathode, $E_c$	70.0 kJ/mol	[24]
Transfer coefficient of HOR, ${\alpha }_a={\alpha }_c$	1	[28]
Transfer coefficient of ORR, ${\alpha }_c$	1	[28]
Reference H2/O2 molar concentration, $C^{ref}$	40.88 mol/m3	[28]
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	[24]
Faraday’s constant, $\mathrm{F}$	96,485 C/mol	[28]
Universal gas constant, $R_u$	8.314 $\mathrm{J}/\left(\mathrm{m}\mathrm{o}\mathrm{l}\bullet \mathrm{K}\right)$
H2 diffusivity in the anode gas channel, $D_{0,H_2,a}^g$	1.1028 × 10−4 m2/s	[28]
H2O diffusivity in the anode gas channel, $D_{0H_{2}O,a}^g$	1.1028 × 10−4 m2/s	[28]
O2 diffusivity in the cathode gas channel, $D_{0,O_2,c}^g$	3.2348 × 10−4 m2/s	[28]
H2O diffusivity in the cathode gas channel, $D_{0,H_{2}O,c}^g$	7.35 × 10−5 m2/s	[28]
Binary gas diffusivity ($D_i^g)$	For nonporous regions $D_i^g=D_o^g{\left(\frac{T}{T_0}\right)}^{3/2}\left(\frac{P_0}{P}\right)$	(11)
Effective diffusivity ($D_i^{g,eff})$	For porous regions $D_i^{g,eff}=\frac{\epsilon }{{\tau }_k}{D}_i^g$	(12)

compiles constitutive equations detailing the kinetic and physicochemical properties of various PEMFC components. For the M2 model,

Table 4. Expressions used in the two-phase mixture model [17].

Description	Expression
Mixture density ($\mathsf{\rho })$	$\mathsf{\rho }={\mathsf{\rho }}^{l}s+{\rho }^g\left(1-s\right)$	(13)
Gas mixture density ${(\rho }^g)$	${\rho }^g=\left(\frac{P}{R_{u}T}\right)\frac{1}{{\sum }_i\frac{m_i^g}{M_i}}$	(14)
Mixture velocity ($\rho \overrightarrow{u})$	$\rho \overrightarrow{u}={\rho }^l{\overrightarrow{u}}^l+{\rho }^g{\overrightarrow{u}}^g$	(15)
Mixture mass fraction	$m_i=\frac{{\rho }^{l}s{m}_i^l+{\rho }^g\left(1-s\right)m_i^g}{\rho }$	(16)
Relative permeability	$k_r^l=s^3$	(17)
	$k_r^g=(1-s{)}^3$	(18)
Kinematic viscosity of the two-phase mixture	$v={\left(\frac{k_r^l}{v^l}+\frac{k_r^g}{v^g}\right)}^{-1}$	(19)
Kinematic viscosity of the gas mixture	$v^g=\frac{{\mu }^g}{{\rho }^g}=\frac{1}{{\rho }^g}{\sum }_{i=1}^n\frac{x_i{\mu }_i}{{\sum }_{j=1}^n{x}_j{\varphi }_{ij}}$	(20)
	where ${\varphi }_{ij}=\frac{1}{\sqrt{8}}{\left(1+\frac{M_i}{M_j}\right)}^{-1/2}{\left[1+{\left(\frac{{\mu }_i}{{\mu }_j}\right)}^{1/2}{\left(\frac{M_j}{M_i}\right)}^{1/4}\right]}^2$ and	(21)
	${\mu }_i\left[\mathrm{N}\bullet \mathrm{s}\bullet {\mathrm{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	(22)
Relative mobility	${\lambda }^l=\frac{k_r^l}{v^l}v$	(23)
	${\lambda }^g=1-{\lambda }^l$	(24)
Diffusive mass flux	$\overrightarrow{j^l}={\rho }^l{\overrightarrow{u}}^l-{\lambda }^l\rho \overrightarrow{u}=\frac{K}{v}{\lambda }^l{\lambda }^g𝛻P_c$	(25)
Capillary pressure Pc	$P_c=P^g-P^l=\sigma \mathit{cos}\theta {\left(\frac{\epsilon }{K}\right)}^{1/2}J\left(s\right)$	(26)
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.120s^2+1.263s^3\end{array}\right.\begin{array}{c}\mathrm{if}{\theta }_c<90^{\circ }\\ \mathrm{if}{\theta }_c>90^{\circ }\end{array}$	(27)

provides an overview of the two-phase mixture properties, while

Table 5. Transport properties in the electrolyte.

Description	Expression
Membrane water content (λ)	Water $\lambda =\left\{\begin{array}{c}{\lambda }_g=0.043+17.81a-39.85a^2+36.0a^{3}for0<a\le 1\\ {\lambda }_l=22\end{array}\right.$	(28)
	Water activity, $a=\frac{C_w^g{R}_{u}T}{P_{sat}}$	(29)
Electro-osmotic drag (EOD) coefficient of water $\left(n_d\right)$	$n_d=\frac{2.5\lambda }{22}$	(30)
Proton conductivity ($\kappa$)	$\kappa =\left(0.5139\lambda -0.326\right)\mathit{exp}\left.\left[1268\left(\frac{1}{303}-\frac{1}{T}\right)\right.\right]$	(31)
Water diffusion coefficient ($D_w^{mem}$)	$D_w^{mem}=\left\{\begin{array}{c}2.692661843\bullet {10}^{-10}for\lambda \le 2\\ \left\{0.87\left(3-\lambda \right)+2.95(\lambda -2)\}\bullet {10}^{-10}\bullet e^{\left(7.9728-2416/T\right)}for2<\lambda \le 3\right.\\ \left\{2.95\left(4-\lambda \right)+1.642454\left(\lambda -3\right)\right\}\bullet {10}^{-10}\bullet e^{\left(7.9728-2416/T\right)}for3<\lambda \le 4\\ \left(2.563-0.33\lambda +0.0264{\lambda }^2-0.000671{\lambda }^3\right)\bullet {10}^{-10}\bullet e^{\left(79728-2416/T\right)}for4<\lambda \le {\lambda }_{a=1}^g\end{array}\right.$	(32)
Interfacial resistance of the water film	${\mathsf{\Omega }}_{w,int}=z_w\frac{{\delta }_w}{D_{O_2,w}}$	(33)

outlines information on species and transport properties.

2.3. Microscale CL Model

The schematic representation of the micro CL model is illustrated in Figure 1, while

Table 6. Correlations of the microscale CL model [17].

Description	Expression
Carbon volume fraction (${\epsilon }_{Pt}$)	${\epsilon }_{Pt}=\frac{1}{{\rho }_{Pt}}\frac{L_{Pt}}{{\delta }_{CL}}$	(34)
Pt particle volume fraction (${\epsilon }_C$)	${\epsilon }_c=\frac{L_{Pt}}{{\delta }_{CL}}\left(\frac{1}{{\rho }_{c}wt\%}\right)$	(35)
Ionomer volume fraction (${\epsilon }_e$)	${\epsilon }_e=\left(\frac{I}{C}\right)\frac{{\rho }_c}{{\rho }_e}{\epsilon }_c$	(36)
CL thickness (${\delta }_{CL}$)	${\delta }_{CL}=\left[\frac{\left(\frac{I}{C}\right)\frac{{\rho }_C}{{\rho }_e}\left(\frac{1}{{\rho }_C\times wt\%}\right)+\left(\frac{1}{{\rho }_{Pt}}+\frac{1}{{\rho }_C\times wt\%}\right)}{1-{\epsilon }_{CL}}\right]L_{Pt}$	(37)
Thickness of the ionomer film (${\delta }_e$)	${\delta }_e=\left[{\left(1+\frac{{\epsilon }_e}{{\epsilon }_c}\right)}^{\frac{1}{3}}-1\right]r_c$	(38)
Thickness of water film (${\delta }_w$)	${\delta }_w=\left[{\left(1+\frac{{\epsilon }_e}{{\epsilon }_c}+\frac{s}{{\epsilon }_c}\right)}^{\frac{1}{3}}-1\right]r_c-{\delta }_e$	(39)
Oxygen balance in a single Pt/C particle domain	$C_{O_2}^g-C_{O_2}^{Pt}=-{\Omega }_T\frac{j_c}{4F{\bullet a}_c}={\Omega }_T\frac{I}{4F{\bullet {\delta }_{CL}\bullet a}_c}=R_T\frac{I}{4F}$	(40)
Volumetric surface area of ionomer ($a_c)$	$a_c=\frac{3{\epsilon }_c}{r_c^3}{\left(r_c+{\delta }_e+{\delta }_w\right)}^2$	(41)
Volumetric surface area of Pt ($a_{Pt})$	$a_{Pt}=\frac{a_{ECSA}{L}_{Pt}}{{\delta }_{CL}}$	(42)
	${(a}_{ECSA}=\frac{3}{r_{Pt}{\bullet \rho }_{Pt}})$	(43)
Total oxygen transport resistance (${\mathsf{\Omega }}_T)$	${\Omega }_T={\Omega }_{w,int}+\frac{{\delta }_w^{eff}}{D_{O_2,w}}+{\Omega }_{e,int}+\frac{{\delta }_e^{eff}}{D_{O_2,e}}+{\Omega }_{Pt,int}$	(44)
Total transport resistance	$R_T=\frac{{\mathsf{\Omega }}_T}{{{\delta }_{CL}\bullet a}_c}$	(45a)
Interfacial resistance in ionomer	${\Omega }_{e,int}=z_e\frac{{\delta }_e}{D_{O_2,e}}$	(45b)
Interfacial resistance in water film	${\Omega }_{w,int}=z_w\frac{{\delta }_w}{D_{O_2,w}}$	(45c)
Interfacial resistance in platinum particle	${\Omega }_{Pt,int}=z_{Pt}\frac{{\left(r_c+{{\delta }_{Ti{O}_2}+\delta }_e\right)}^2}{r_{Pt}^2}\frac{{\rho }_{Pt}}{{\rho }_c}{\left(\frac{r_{Pt}}{r_c}\right)}^3\left(\frac{1-wt\%}{wt\%}\right)\frac{{\delta }_e}{D_{O_2,e}}$	(45d)
Effective ionomer thickness	${\delta }_e^{eff}=\frac{{\left(r_c+{\delta }_e+{\delta }_w\right)}^2}{{n_{Pt}\bullet r}_{Pt}^2}{\delta }_e$	(46)
Effective water thickness	${\delta }_w^{eff}=\frac{{\left(r_c+{\delta }_e+{\delta }_w\right)}^2}{{n_{Pt}\bullet r}_{Pt}^2}{\delta }_w$	(47)
Number of Pt particles on a single carbon particle	$n_{Pt}=\frac{{\rho }_c}{{\rho }_{Pt}}{\left(\frac{r_c}{r_{Pt}}\right)}^3\left(\frac{wt\%}{1-wt\%}\right)$	(48)
Oxygen concentration on the Pt particle surfaces $\left(C_{O_2}^{Pt}\right)$	$f\left(C_{O_2}^{Pt}\right)=B{\left(C_{O_2}^{Pt}\right)}^{\gamma }+\frac{4F}{{\Omega }_T}{\frac{3{\epsilon }_c}{r_c^3}{\left(r_c+{\delta }_e\right)}^{2}C}_{O_2}^{Pt}-\frac{4F}{{\Omega }_T}{\frac{3{\epsilon }_c}{r_c^3}{\left(r_c+{\delta }_e\right)}^{2}C}_{O_2}^g=0$             $where$ $B=\frac{3L_{Pt}}{{r_{Pt}{\bullet \rho }_{Pt}\bullet \delta }_{CL}}{i}_{0,c}^{ref}{\left(\frac{1}{C_{O_2,ref}}\right)}^{3/4}exp\left(-\frac{E_c}{R}\left(\frac{1}{T}-\frac{1}{353.15}\right)\right)exp\left(-\frac{{\alpha }_c}{RT}F{\eta }_c\right)$	(49)

provides a comprehensive list of parameters for a CL particle within a microscale CL domain. These parameters encompass volume fractions such as Pt (${\epsilon }_{Pt}$), carbon (${\epsilon }_C$), ionomer (${\epsilon }_e$), and ${\delta }_{CL}$. The values were determined based on factors like the ionomer to carbon weight ratio (I/C), weight percent of Pt/C, and Pt loading. Given the hydrophilic nature of the ionomer layer, a layer of water covers it, and the thickness of the ionomer and water films was determined using Equations (38) and (39), respectively. Furthermore, Equations (41)–(43) were employed to calculate parameters such as the oxygen flux across ionomer layers and the water, the volumetric surface area of the ionomer ($a_C$), $j_c$ for the oxygen reduction reaction (ORR), the active volumetric surface area of Pt ($a_{Pt}$), and the electrochemical active surface area (ECSA). The interfacial resistance in the water and ionomer films, along with the effective thicknesses of the water and ionomer films (${\delta }_w^{eff}$, ${\delta }_e^{eff}$), were utilized to determine the total oxygen transport resistance, denoted by the variable ${\Omega }_T$. For a more in-depth understanding of microscale modeling, interested readers are recommended to consult our previous works [17,29].

2.4. Boundary Conditions and Numerical Implementation

Figure 1 illustrates the computational domain for a single-channel PEMFC geometry, including the operating and boundary conditions. Additional details can be found in

Table 7. Boundary and operating conditions [17].

Description	Expression
Anode inlet velocity	$u_{in,a}=\frac{{\xi }_a\frac{I}{2F}{A}_{mem}}{C_{H_2}{A}_{a,ch}}$	(50)
Cathode inlet velocity	$u_{in,c}=\frac{{\xi }_c\frac{I}{4F}{A}_{mem}}{C_{O_2}{A}_{c,ch}}$	(51)
Constant temperature on side walls	$T_{side}=60°C$	(52)
Heat flux on the top and bottom surfaces	$k^{eff}\nabla T=0$	(53)
Electric potential on anode outer BP surface	${\varphi }_s=0$	(54)
Electric potential on cathode outer BP surface	${\varphi }_s=V_{cell}or\sigma 𝛻{\varphi }_s=I$	(55)
Pressure(anode/cathode), $P_a/P_c$	2/2 bar
Operating temperature, $T_{cell}$	333.15 K
Stoichiometry(anode/cathode), ${\xi }_a/{\xi }_c$	1.5/2.0

. The anode and cathode inlet velocities, formulated based on the stoichiometric ratios (${\xi }_a\mathrm{a}\mathrm{n}\mathrm{d}{\xi }_c$, respectively), are provided. Mass flow conditions, enforcing no-slip and impermeability criteria, were applied to all exterior surfaces, except for the anode and cathode inlet and outlet areas. Within the computational domain, isothermal boundary conditions were applied to the side walls, while adiabatic boundary conditions were imposed on the bottom and top surfaces of the anode and cathode, suitable for the stack environment. The PEMFC can operate under either galvanostatic or potentiostatic mode. Consequently, a constant voltage or current density was applied to the cathode outer sidewall, while the electric potential ${\varphi }_s$ was fixed at zero on the anode sidewall, as depicted in Figure 1. The number of mesh elements for the single channel geometry, covering an MEA surface area of 2 × 25 mm², was estimated to be around 2 hundred thousand based on the grid-independent study results represented in Figure 2. Employing user-defined functions, the PEMFC model was numerically integrated into the commercially available CFD program ANSYS Fluent ver. 23 (ANSYS, Inc., Canonsburg, PA, USA). Convergence criteria for equation residuals were set at 10⁻⁸.

3. Reliability Based Design Optimization

In most engineering problems, design optimization has been used to derive theoretically optimal design performance that can be manufactured in the real industry. In other words, the results of the optimization process through the CFD simulation which can be directly related to industrial manufacturing should include considering uncertainties that might appear in the manufacturing process such as tolerances. As design variables in the manufacturing process have uncertainties, engineers should design products that satisfy the constraints of engineering problems as well as reliability in terms of the quality of products. To obtain the optimal point which satisfied both constraints and reliability, RBDO process [30] was suggested in this paper and the general schematic is shown in Figure 2.

3.1. General Formulation of RBDO

In general, the definition of RBDO problem can be mathematically expressed as follows:

$$Minimizef\left(\mathrm{d}\right)$$

(56)

$$SubjecttoP\left\{G_j\left(\mathrm{X}\right)\le 0\right\}<P_{f_j}^{Tar},j=1,\dots ,nc$$

(57)

$${\mathrm{d}}_{\mathrm{L}}\le \mathrm{d}\le {\mathrm{d}}_{\mathrm{U}},\mathrm{d}\in {\mathrm{R}}^{\mathrm{n}\mathrm{d}},\mathrm{X}\in {\mathrm{R}}^{\mathrm{n}\mathrm{r}}$$

(58)

Here, $f$ denotes the objective function, $\mathrm{d}$ denotes the vector of design variables which means the average of the vector of random variables, $\mathrm{X}$. $G_j\left(\mathrm{X}\right)$ represents the jth deterministic constraint functions, thus if $G_j\left(X_i\right)$ is less then 0, the design point $X_i$ would be structural/performance failure point in terms of jth constraint $G_j$. Moreover, $P$ represents the probability of failure, thus $P_{f_j}^{Tar}$ means target probability of failure in terms of jth constraint. ${\mathrm{d}}_{\mathrm{L}}$ and ${\mathrm{d}}_{\mathrm{U}}$ represent the upper bound and lower bound of the vector of design variables, respectively. nc and nd denote the number of probabilistic constraints and number of design variables, respectively, and nr denotes the number of random variables. Equation (57) expresses the probabilistic constraint of RBDO which includes probabilistic conditions. In order to determine whether Equation (57) is satisfied, the probability of failure should be calculated. In general, the probability of failure can be calculated as follows:

$$P_F=P\left\{\mathrm{X}\in {\mathsf{\Omega }}_{\mathrm{F}}\right\}=P\left\{G_j\left(\mathrm{X}\right)\le 0\right\}=\underset{{\mathsf{\Omega }}_{\mathrm{F}}}{\overset{}{\int }}{h}_{\mathrm{X}}\left(\mathrm{X}\right)d\mathrm{X}=E\left[I_{{\mathsf{\Omega }}_{\mathrm{F}}}\left(\mathrm{X}\right)\right]$$

(59)

Here, $E[\bullet ]$ denotes the expectation operator and $h_{\mathrm{X}}\left(\mathrm{X}\right)$ denotes the joint probability densigy function (JPDF) of all random variables. ${\mathsf{\Omega }}_{\mathrm{F}}$ and $I_{{\mathsf{\Omega }}_{\mathrm{F}}}$ which means that the failure sets are expressed as follows:

$${\mathsf{\Omega }}_{\mathrm{F}}=\left\{\mathrm{X}:G_j\left(\mathrm{X}\right)\le 0\right\}$$

(60)

$$I_{{\mathsf{\Omega }}_{\mathrm{F}}}\left(\mathrm{X}\right)=\left\{\begin{array}{c}1,if{G}_j\left(\mathrm{X}\right)\le 0\\ 0,otherwise\end{array}\right.$$

(61)

Stochastic sensitivity can be obtained from partial differentiation about $P_F$ and is expressed as follows:

$$\frac{{\partial P}_F\left(\mathrm{d}\right)}{\partial d_i}=\frac{\partial }{\partial d_i}\underset{{\mathsf{\Omega }}_{\mathrm{F}}}{\overset{}{\int }}{h}_{\mathrm{X}}\left(\mathrm{X}\right)d\mathrm{X}$$

(62)

3.2. Sampling-Based RBDO

In the case of problems where physical phenomena were complexly overlapped, multiple integrations on JPDF must be performed in the failure region for reliability analysis, which cannot be calculated through mathematical methods [31]. As an alternative, sampling-based methods have been developed along with various approximate analysis methods such as the second-order reliability method (SORM) or the first-order reliability method (FORM) [32,33,34,35]. In particular, sampling-based methods such as Monte Carlo simulation (MCS) can be applied to many optimization problems regardless of the complexity of the problem [36]. When the probability of failure is obtained through sampling-based methods, Equations (59) and (60) can be modified as follows [37]:

$$P\left\{G_j\left(\mathrm{X}\right)\le 0\right\}\cong \frac{1}{M}\sum _{k=1}^M{I}_{{\mathsf{\Omega }}_{\mathrm{F}}}\left({\mathrm{X}}^{\left(k\right)}\right)$$

(63)

$$\frac{{\partial P}_F\left(\mathrm{d}\right)}{\partial d_i}\cong \frac{1}{M}\sum _{k=1}^M{I}_{{\mathsf{\Omega }}_{\mathrm{F}}}\left({\mathrm{X}}^{\left(k\right)}\right)s_{d_i}^{\left(1\right)}\left({\mathrm{X}}^{\left(k\right)}\right)$$

(64)

Here, ${\mathrm{X}}^{\left(k\right)}$ is kth realization of X and M denotes number of samples of MCS. $s_{d_i}^{\left(1\right)}\left(\mathrm{X}\right)$ represents the first-order score function, if design variables are statistically independent each other, it can be expressed as follows [37]:

$$s_{d_i}^{\left(1\right)}\left(\mathrm{X}\right)\equiv \frac{\partial \mathrm{l}\mathrm{n}{h}_{\mathrm{X}}\left(\mathrm{X}\right)}{\partial d_i}=\frac{\partial \mathrm{l}\mathrm{n}{h}_{X_i}\left(X_i\right)}{\partial d_i}$$

(65)

The marginal probability density function (Marginal PDF) and the cumulative distribution function (CDF) of normal distribution can be obtained analytically, and are shown in Equations (66) and (67) respectively.

$$h_{\mathrm{X}}\left(\mathrm{X}\right)=\frac{1}{\sqrt{2\pi \sigma }}{e}^{-0.5{\left[\frac{x-\mu }{\sigma }\right]}^2}$$

(66)

$$H_{\mathrm{X}}\left(\mathrm{X}\right)=\frac{1}{\sqrt{2\pi }}{\int }_{-\mathrm{\infty }}^{\frac{x-\mu }{\sigma }}\exp \left(-\frac{1}{2}{\zeta }^2\right)d\zeta$$

(67)

Therefore, using the sampling-based analysis method, the probability of failure or sensitivity of the probabilistic constraint can be directly obtained through Equations (63) and (64). More importantly, the actual values for M samples of MCS should be obtained, however it might be impossible to obtain the actual value through experiment. In the case of simulation, it is extremely inefficient to obtain simulation results at M points in every step. Hence, a surrogate model was employed for replacing CFD simulation or experiments [38,39]. The general flow of sampling-based RBDO using a surrogate model is shown in Figure 3.

3.3. Metamodeling (Surrogate)

MLP is considered one of the best surrogate models based on a neural network’s structure. Over the past few years, the model has been seen to be used exclusively for predictions in the field of product design and engineering. The MLP as seen in Figure 4a is generally a three-layered model which consists of the input layer, hidden layer, and output layer. In the present study, an MLP having a sigmoid activation function in the hidden layer was constructed, trained, and tested. The complete framework of the model and the set of equations that constitute them are explained as follows:

$$sigmoid=\frac{1}{1+e^{-x}}$$

(68)

Since the MLP works based on matrix calculations in general, the inputs which usually consist of training data are represented in the form of a matrix which is presented in Equation (69).

$$X=\left[\begin{array}{ccc}{x}_{11}& \cdots & x_{N1}\\ ⋮& \ddots & ⋮\\ x_{1d}& \cdots & x_{Nd}\end{array}\right]$$

(69)

where, $X$ represents the matrix of $N$ training data and $d$ represent the number of design variables considered in the optimization problem.

The second layer, also known as the hidden layer, consists of nodes that are connected to the nodes of the input layer and have weights between them. The weights are represented as a matrix and are given in Equation (70) as follows,

$$W_{input-hiddenlayer}=\left[\begin{array}{ccc}{w}_{11}& \cdots & w_{d1}\\ ⋮& \ddots & ⋮\\ w_{1p}& \cdots & w_{dp}\end{array}\right]$$

(70)

where, $p$ represents the number of weights in the hidden layer.

The movement of data from the input to the hidden layer is known as the feed-forward process. The input to the hidden layer is nothing but the product of weights and the training data and is given in Equation (71).

$$net=\left[\begin{array}{ccc}{w}_{11}& \cdots & w_{d1}\\ ⋮& \ddots & ⋮\\ w_{1p}& \cdots & w_{dp}\end{array}\right]\left[\begin{array}{c}{x}_{n1}\\ x_{n2}\\ ⋮\\ x_{nd}\end{array}\right]=\left[\begin{array}{c}\begin{array}{c}{x}_{n1}{w}_{11}+x_{n2}{w}_{21}+\dots x_{nd}{w}_{d1}\\ x_{n1}{w}_{12}+x_{n2}{w}_{22}+\dots x_{nd}{w}_{d2}\\ ⋮\end{array}\\ x_{n1}{w}_{1p}+x_{n2}{w}_{2p}+\dots x_{nd}{w}_{dp}\end{array}\right]$$

(71)

As discussed earlier, the nodes of the hidden layer are embedded with a sigmoid activation function as a result of which, the data ($net$) entering the nodes of the hidden layer become activated and can be represented in Equation (72) as follows:

$$h=sigmoid\left(net\right)=\left[\begin{array}{c}{h}_1\\ h_2\\ ⋮\\ h_p\end{array}\right]=\left[\begin{array}{c}\begin{array}{c}sigmoid(x_{n1}{w}_{11}+x_{n2}{w}_{21}+\dots x_{nd}{w}_{d1})\\ sigmoid(x_{n1}{w}_{12}+x_{n2}{w}_{22}+\dots x_{nd}{w}_{d2})\\ ⋮\end{array}\\ sigmoid(x_{n1}{w}_{1p}+x_{n2}{w}_{2p}+\dots x_{nd}{w}_{dp})\end{array}\right]$$

(72)

where, $h$ represents the output of the hidden layer upon activation. The output $h$ of the hidden layer is further processed to the output layer by linearly combining them with the weights between the output layer and the hidden layer. The weights between the output and hidden layer are given in Equation (73) as,

$$W_{hiddenlayer-output}=\left[\begin{array}{ccc}{w}_{11}& \cdots & w_{p1}\\ ⋮& \ddots & ⋮\\ w_{m1}& \cdots & w_{mp}\end{array}\right]$$

(73)

Furthermore, the linear summation of $W_{hiddenlayer-output}\times h$ to the output layer is activated by a sigmoid function and then presented as the outputs of the MLP and is given in Equation (74) as follows:

$$outpu{t}_k=sigmoid(w_{k1}{h}_1+w_{k2}{h}_2+\dots w_{kp}{h}_p)$$

(74)

where, $k$ denotes the $k^{th}$ objective function.

The learning process of MLP includes multiple feed-forward and back-propagation processes. Each pass of feed-forward and back-propagation is called an epoch that can be considered a hyperparameter. The model is able to minimize the error and predict approximate results by performing multiple epochs and optimizing the model’s hyperparameters like the learning rate and number of nodes in the hidden layers.

4. Results and Discussion

4.1. Experimental Validation of a 3D PEMFC Model and Construction of the MLP Models

The 3D, multi-scale, and multi-phase PEMFC model, delineated in Section 2, underwent validation against experimental data to establish its accuracy. Figure 5 provides visual depictions of the intricate experimental configuration and the individual cell hardware. This PEMFC single cell consisted of anode and cathode GDLs, a catalyst-coated membrane (CCM), gaskets, and BPs, which included a machined flow field. The anode and cathode GDLs employed SGL 25 BC with a thickness of 235 µm, provided by SGL Carbon, Germany. The CCM, supplied by Vinatech in South Korea, had an overall thickness of 38 µm. It comprised a PEM (Nafion NR212 with a thickness of 20 µm), and anode and cathode catalyst CLs with a Pt loading of 0.4 mg/cm². Teflon gaskets and sub-gaskets were applied to seal the cell (active area was 3 × 3 cm²), with thicknesses of 105 µm and 170 µm, respectively. These PEMFC components were assembled with a clamping pressure of 20 kgf/cm², taking into account the compression of the GDLs to prevent any leakage. The experimental conditions are represented in Table 7 and the 3D PEMFC simulation was conducted to secure the validity of the model. As depicted in Figure 6, the PEMFC model precisely replicated the polarization characteristics of the experimental cell, thereby substantiating the model’s accuracy and validity. This validated 3D PEMFC model was subsequently employed to produce sample data essential for the training of MLP-based surrogate models. The MLP was trained within a defined effective design space, which is characterized as follows.

$$Maximize{V}_{cell}=f(w_{land},w_{ch},d_{ch})$$

(75)

$$Minimize∆P=f(w_{land},w_{ch},d_{ch})$$

(76)

$$subjectedto\left\{\begin{array}{c}{w}_{land},0.3\mathrm{m}\mathrm{m}\le w_{land}\le 1.00\mathrm{m}\mathrm{m}\\ w_{ch},0.3\mathrm{m}\mathrm{m}\le w_{ch}\le 2.00\mathrm{m}\mathrm{m}\\ d_{ch},0.3\mathrm{m}\mathrm{m}\le d_{ch}\le 2.00\mathrm{m}\mathrm{m}\end{array}\right.$$

(77)

After establishing the optimization problem, our focus shifted to the construction of surrogate models based on MLP networks. Using Latin hypercube sampling (LHS) depicted in Figure 7, we generated 121 sample points. For each, the cell voltage ($V_{cell}$), pressure drop ($∆P$), and standard deviation of current density distribution (${\sigma }_I$) were computed via a comprehensive 3D, two-phase, multi-scale PEMFC model. In particular, ${\sigma }_I$ was estimated by using local current density values ($I_i$) within the membrane as follows:

$${\sigma }_I=\sqrt{\frac{{\sum }_{i=1}^N{\left(I_i-I_{avg}\right)}^2}{N-1}}$$

(78)

The resulting dataset, comprising V_cell, ∆P, and σ_I, was divided into training (100 samples) and testing (21 samples) subsets to train and validate three surrogate models, respectively, for each output variable. Model performance was evaluated using adjusted R-squared and root mean square error (RMSE) metrics, as shown in Figure 7. As expounded in Section 3.3, the MLP models demonstrated robust prediction capabilities for highly nonlinear datasets compared to alternative surrogate models. The training set achieved an adjusted R-squared value of 0.9897 and RMSE values of 25 mV for V_cell, 0.9921 and 0.8464 Pa/cm for ∆P, and 0.9324 and 0.0694 A/cm² for σ_I. Furthermore, the MLP’s efficacy was tested against 21 randomly generated LHS data points, with the error analysis presented in Figure 7. The test data exhibited an adjusted R-squared value of 0.9961 with RMSE values of 8 mV for V_cell, 0.9920 and 0.4979 Pa/cm for ∆P, and 0.9589 and 0.0462 A/cm² for σ_I. These error metrics decisively confirmed the MLP’s adeptness at processing nonlinear data.

4.2. Comparative Analysis of Single-Objective and Multi-Objective Optimization Methods

The procedure to identify the optimal design point on the response surface with the aid of a trained MLP model was executed. Through simulation, the response surface, which was informed by the single cell voltage objective function, was transformed to represent stack power density, ${\dot{P}}_{stack}$ in kW/L, incorporating geometric considerations. In order to optimize the algorithm, PSO was employed to find the optimal point. As an algorithm optimization, PSO provides efficient performance. In detail, in PSO, particles share information with each other and seek the optimal solution, thus even if a particle is trapped in the local optimal, it can escape through particles including the global optimal. Through PSO, a deterministic optimization algorithm, the optimal point is represented in Table 1. The projected performance at this design point was 3.252 kW/L at the operating current density of 1.43 A/cm², demonstrating a marked improvement of approximately 1.585 kW/L over the baseline stack power density. This enhancement stems from a substantial narrowing of the land width ($w_{land}$) within the parallel channel flow configuration, which improved the oxygen transfer beneath the land to the reactive zones, thereby diminishing concentration overpotentials. Notably, the channel depth ($d_{ch}$) was reduced to less than a third of that of the reference stack, yielding gains not only in oxygen supply efficiency to the GDL but also in reducing the stack’s overall volume. Nonetheless, this decrease in $d_{ch}$ led to an increased pressure drop ($∆P$), highlighting the necessity for multi-objective optimization to balance the enhancement of stack power density against pressure drop mitigation. A decrease in $d_{ch}$ is conducive to stack power density improvements due to volume reduction, yet it exacerbates the $∆P$, suggesting a potential trade-off between these objectives. This interplay is analyzable via the construction of a Pareto front using the NSGA-2 multi-objective algorithm. The fidelity of the MLP model for pressure drop, employed in the Pareto front construction, is corroborated by Figure 7b, and the resultant Pareto front is depicted in Figure 8.

As observed in Figure 8, the optimal solution obtained via single-objective optimization resided at the extremity among the Pareto solutions, implying neglect of pressure drop considerations. For instance, at a 4 Pa/cm pressure drop as depicted in the optimum point via multi-objective optimization, ${OP}_{multi}$, the stack power density marginally decreased from 3.252 kW/L to roughly 2.821 kW/L relative to the optimum point via single-objective optimization, ${OP}_{sgl}$, albeit with a substantial pressure drop advantage of approximately 26 Pa/cm. This observation indicates that an optimal solution predicated solely on stack power density might not always be suitable, thereby affirming the utility of multi-objective optimization that incorporates pressure drop as a concurrent objective function, offering a more viable optimal design alternative.

4.3. RBDO Accounting for Output Constraints and Production Tolerances of Design Variables

To identify a design point that was both optimal and realistic, we employed three MLP-based surrogate models for ${\dot{P}}_{stack}$, $∆P$, and ${\sigma }_I$, as illustrated in Figure 7. Furthermore, we established the following constraints for $∆P$ and ${\sigma }_I$.

$$Maximize{\dot{P}}_{stack}\left(w_{land},w_{ch},d_{ch}\right)$$

(79)

$$Subjectto\left\{\begin{array}{c}∆P<8\mathrm{P}\mathrm{a}/\mathrm{c}\mathrm{m}\\ {\sigma }_I<0.3\mathrm{A}/{\mathrm{c}\mathrm{m}}^2\end{array}\right.$$

(80)

Following Equations (79) and (80), we performed constraint optimization on the response surfaces developed using three MLP-based surrogate models wherein the constrained optimal solution was traced using the SQP algorithm. Within the design space, the constrained optimum point (COP) is detailed in Table 1, and the stack power density was predicted to be 3.008 kW/L. Although this was approximately a 7.51% decrease compared to the stack power density at the single-objective based deterministic optimum point ($DOP$), it can be considered a realistic design value when other performance indicators ($∆P$ and ${\sigma }_I$) are taken into account. Figure 9 represents the contour of pressure drop and current density distribution in DOP and COP. In terms of current density, both cases showed a similar trend along the channel, thus the standard deviation of current density obtained was 0.2217 in DOP and 0.2269 in COP, respectively. However, in terms of pressure drop, the COP showed it had significantly improved compared with the DOP case. Constraint in pressure drop was set around 8 Pa, the maximum pressure in COP represents under 8 Pa.

Figure 10 displays the 2D surface plots for each design variable. It was observed that the COP did not lie within the feasible region that satisfied the constraints for both $∆P$ and ${\sigma }_I$. Instead, it was located at the interface between the infeasible and feasible regions.

Considering the dimensional tolerances inherent in actual component production, the theoretically derived optimal values of the design variables, as identified at the COP in

Table 8. Comparative summary of simulation results from various optimization techniques at $I=1.43\mathrm{A}/{\mathrm{c}\mathrm{m}}^2$.

Design Cases	Design Variables			${\dot{\mathit{P}}}_{\mathit{s}\mathit{t}\mathit{a}\mathit{c}\mathit{k}}$ [kW/L]	$\mathit{∆}\mathit{P}$ [Pa/cm]	Reliability
	${\mathit{w}}_{\mathit{l}\mathit{a}\mathit{n}\mathit{d}}$ [mm]	${\mathit{w}}_{\mathit{c}\mathit{h}}$ [mm]	${\mathit{d}}_{\mathit{c}\mathit{h}}$ [mm]			PG1 [%]	PG2 [%]
Ref. point	1.00	1.00	1.00	1.667	2.832	100	0.00
Deterministic Optimum point (DOP)	0.30	1.48	0.30	3.252	30.209	0.00	100
Constraint Optimum point (COP)	0.30	1.51	0.50	3.008	8.000	49.79	100
RBDO $\left({OP}_{80\%}\right)$	0.30	1.46	0.53	2.974	6.918	80.04	99.99
RBDO $\left({OP}_{90\%}\right)$	0.30	1.51	0.54	2.958	6.453	90.21	99.99
RBDO $\left({OP}_{90\%}\right)$	0.30	1.55	0.58	2.918	5.519	99.13	99.99

, may not be practically reliable. This unreliability highlights the inherent limitations of multi-objective optimization. It is particularly critical to note that the design of PEMFC components, with dimensions varying from several micrometers to tens of millimeters, is sensitive to minor production variances; these can significantly influence PEMFC performance.

Consequently, this study implemented an additional probabilistic optimization process that accounted for potential production errors, thereby enhancing the reliability of the results. In mass production, it is assumed that the distribution of errors in design variables follows a normal distribution centered around the optimal values at the COP as depicted in Figure 11. Process capability quantifies the variability in the quality of products from a well-managed process, with the process capability index ($C_p$) serving as a method for this quantification [41,42]. This index is a statistical measure that assesses if a process is capable of consistently producing parts of high quality within predefined specification limits. It calculates the ratio between the permissible variation, namely the specified limits of quality, and the actual variation observed in the process. When specific upper and lower specification limits are established, the process capability can be mathematically expressed, typically illustrated in Equation (81),

$$C_p=\frac{USL-LSL}{6\sigma }$$

(81)

where USL and LSL denote the upper and lower specification limits for process quality. When $C_p$ exceeds 1, the output of the process is well within the acceptable range, surpassing the defined specification limits. A $C_p$ of exactly 1 signifies that the process output is at the threshold, precisely meeting the specification limits. Conversely, a $C_p$ below 1 indicates that the process output fails to achieve the established specification limits. In this study, we established the minimum acceptable criterion within the permissible range as a $C_p$ of 1. We assumed that the upper and lower specification limits for all design variables ($w_{land},w_{ch},d_{ch}$) were ±0.1 mm. Based on these limits, we calculated the standard deviation ($\sigma$) to determine the distribution of the input variables.

At this design point, variables adhered to numerical distributions inclusive of tolerances, consequently shaping the distribution of the objective function values. This distribution reflected the probability of constraint satisfaction at the optimal point, as illustrated in Figure 12. For the baseline case of an unconstrained optimum (DOP), constraints were ignored, resulting in a distribution within the infeasible region. In contrast, for COP, solutions complied with the constraint surface, yielding approximately half of the solutions within the feasible region, and the remainder distributed across the infeasible domain. Upon elevating the constraint satisfaction probability to 99% for all parameters (i.e., $∆P<20\mathrm{P}\mathrm{a}/\mathrm{c}\mathrm{m}$ and ${\sigma }_I<0.3\mathrm{A}/{\mathrm{c}\mathrm{m}}^2$), a notable shift from COP to reliability based optimum point with target reliability 99% (${OP}_{99\%}$) occurred, detailed in Figure 11 and Table 1. Subsequently, the solutions for ${OP}_{99\%}$ predominantly resided in the feasible region, as shown by the objective function outcomes in Figure 12, with a near-universal compliance (approximately 99%). The reliability calculations indicated a 99.1% confidence level for $∆P<8\mathrm{P}\mathrm{a}/\mathrm{c}\mathrm{m}$ and a propensity to fully meet ${\sigma }_I<0.3\mathrm{A}/{\mathrm{c}\mathrm{m}}^2$. Addressing the errors inherent in actual production led to a verified reduction in stack power density, necessitating a performance compromise from 3.008 kW/L at COP to 2.918 kW/L at ${OP}_{99\%}$, a reduction of about 3.0%. This finding elucidated the intrinsic trade-off between the reliability of objective function constraints and the system’s operational efficacy. Key simulation results at various optimum points are compared in the Section 2.

Moreover, Figure 13 compares and analyzes the performance under constraint conditions at reliabilities of 90% and 80%. The corresponding optimal solutions are detailed in Table 1, from which 50,000 samples were generated, taking into account the tolerances at these optimal points. The results, illustrated in Figure 13a and 13b for 90% and 80% reliability, respectively, showed an increase in the number of design points within the infeasible region when compared to the 99% reliability case. Conversely, as reliability decreased, a marginal improvement in stack power density was noted, with values rising to 2.958 kW/L and 2.974 kW/L for 90% and 80% reliability, respectively. To facilitate a more instinctive comprehension, Figure 14 depicts a 2D distribution of the 50,000 samples according to $∆P$ and ${\dot{P}}_{stack}$, elucidating the trade-off between the reliability of the constraint conditions and the objective function.

5. Conclusions

The study presents a comprehensive 3D, multi-scale, and multi-phase model of PEMFCs, validated against experimental data. The model accurately replicated polarization characteristics, establishing its validity. It was then used to generate sample data for training MLP surrogate models. LHS generated 121 sample points, computing cell voltage, pressure drop, and standard deviation of current density distribution for each. The MLP models showed robust prediction capabilities for nonlinear datasets, with high adjusted R-squared and low RMSE values. An optimal design point was identified using the MLP and PSO to maximize stack power density. This yielded a significant improvement over baseline density by narrowing channel width and reducing channel depth.

However, single-objective optimization prioritized stack power density, while multi-objective optimization considered pressure drop as well. Results showed that solely optimizing for power density may neglect pressure drop concerns, highlighting the importance of multi-objective optimization. Constraint optimization using surrogate models predicted a COP considering both power density and pressure drop, which was 7.51% lower in performance than the DOP. However, the COP may not always lie within feasible regions due to statistical production variances. To address this, a probabilistic optimization process was implemented to account for production errors, assessing process capability using the Cp. This ensured reliability by quantifying the variability in product quality within predefined specification limits, with a Cp of 1 indicating optimal process performance. The study established upper and lower specification limits for design variables and calculated standard deviations to determine the distribution of input variables, enhancing the reliability of results in realistic production.

In terms of production process, at the COP, solutions adhered to constraints, with approximately half within the feasible region. Increasing constraint satisfaction probability to 99% shifted to a reliability-based optimum point (${OP}_{99\%})$, where nearly all solutions resided in the feasible region. Reliability calculations showed a 99.1% confidence level for pressure drop and a tendency to meet current density distribution requirements. Accounting for production errors led to a slight reduction in stack power density from 3.008 kW/L at COP to 2.918 kW/L at ${OP}_{99\%}$. Additionally, the study compared performance at reliabilities of 90% and 80%, observing more design points within the infeasible region as reliability decreased. However, a marginal improvement in stack power density was noted at lower reliabilities. Overall, the analysis highlights the trade-off between constraint reliability and system performance.