In this study, a three-dimensional geometry of a four-serpentine channel PEMFC is built in COMSOL (Figure 1). The geometric parameters are presented in Table 1. The model includes seven components, i.e., cathode gas channel, anode gas channel, gas catalytic layer, GDL, and proton exchange membrane. In the anode catalytic layer, hydrogen is catalytically decomposed into protons, which reach the cathode through the proton exchange membrane. The electrons generated by hydrogen decomposition go to the cathode through the load, which can generate electric power.

The models are assessed using the Batteries and Fuel Cells Module of COMSOL Multiphysics 6.1, and physical fields, such as Secondary Current Distribution, and Transportation of Concentrated Species, are used for the computational analysis. First, the model introduces several simplification assumptions. 1) All processes operate under steady-state conditions; 2) The gas phase adheres to the ideal gas law; 3) Anisotropic electronic conductivities are used in gas diffusion layers, which have conductivities that are about an order of magnitude larger in the in-plane (x and y) directions than in the through-plane (z); 4) We assume a laminar and incompressible flows in the channels and porous layers; 5) The solid-phase electronic potential and membrane ionic potential, respectively, control the transportation of electrons and protons. The Governing equations of the model in Equations 1–16 are discussed as follows:

The convection and diffusion terms for the gaseous species in the channels, gas diffusion layers, and catalyst layers are given as: ∇ · − D i ∇ C i + u j C i = S i , (1) where S is a specific surface area ( m − 1 ), C is the molar concentration ( m o l   m − 3 ), D is diffusivity ( m 2 S − 1 , u is velocity ( m s − 1 ), the subscript i is a gaseous species, and j can be gas or liquid.

The diffusion coefficient of gaseous species in the diffusion layer is expressed as: D = D e f f · ε 1.5 , (2) where the subscript e f f represents effective value, and ε represents a gas pore volume fraction.

Momentum equation for the gaseous species in channels, gas diffusion layers, and catalyst layers: u j · ∇ u j = F − 1 ρ j ∇ P j + μ j ρ j ∇ 2 u j , (3) ∇ · u = 0 , (4) where u is velocity ( m s − 1 ), P is a pressure (Pa), ρ denotes density ( k g   m − 3 ), F is Faraday’s constant ( C   mol − 1 ), μ is a dynamic viscosity ( N   s   m − 2 ).

A continuous function that is based on Darcy’s law is used to compute the average velocity of phase j , which is given as follows: ∇ · ρ j − κ k r , j μ j ∇ P j = H j , (5) u j = − κ μ ∇ P j , (6) where H is a mass source term, κ represents gas permeability, k r is a relative permeability.

Since there is no electrochemical reaction occurring in the gas diffusion layers, the source term of the electronic current has been set to zero, and it is expressed as: − ∇ − σ s ∇ ∅ s = 0 , (7) where ∅ is a cell potential (V), and σ is electronic conductivity ( S m − 1 ).

In the same way, the source term for the ionic current in PEMFC is also zero, which is also given as follow: − ∇ − σ l ∇ ∅ l = 0 . (8)

The catalyst layer is the region in which the reaction takes place, and the ionic balance and electron balance can be defined in the following way: − ∇ − σ l , e f f ∇ ∅ l = − S · i l o c , (9) − ∇ − σ s , e f f ∇ ∅ s = S · i l o c , (10) where i l o c is the local current density. According to the Bulter-Volmer equation, the local current density by the porous electrode reaction in the catalyst layer is defined as: i l o c = i o exp α a F η R T − ⁡ exp − α c F η R T , (11) where i 0 is exchange current density ( A   c m − 2 ), α a and α c are anode and cathode transfer coefficient, respectively, and R is a gas constant m o l   c m 3   s − 1   . M H 2 , in = M H 2 + A I 2 F , (12) M O 2 , in = M O 2 + A I 4 F , (13) where A is the total active area of the MEA ( c m 2 ), M is the molar flow rate ( m o l   s − 1 ). Hydrogen and oxygen consumption are represented by the second term on the right side of Equations 12, 13, respectively. The hydrogen and oxygen molar flow rates that enter the channels are defined by their respective stoichiometric ratios ( λ ). M H 2 , in = λ H 2 + A I 2 F , (14) M O 2 , in = λ O 2 + A I 4 F . (15)

In polymer electrolytes, there is an interconnection between ion transport and water transport. Because of this intercoupling, water content affects polymer electrolyte conductivity. When modeling the water-ion-polymer system, the transport equations for both ions and water molecules must also consider the water-ion friction forces. COMSOL used Weber-Newman’s (2004) concentrated solution theory to show how the three pairs of water molecules, the charge-carrying ion, and the immobilized polymer matrix talk to each other. The model uses the chemical potential of water in the polymer μ 0   J / m o l as dependent parameter, and then the polymer water flux, N 0   m o l / m 2 · s is defined as: N 0 = − σ l ξ F ∇ ∅ l + α + σ l ξ 2 F 2 ∇ μ 0 , (16) where ξ is the electroosmotic coefficient and α   m o l 2 / J · m · s is the water transport coefficient, with the mass balance of ∇ N 0 = 0 in the domain.

The proposed methods involved four main procedures, as shown in the flow chart of Figure 2. Initially, we created a numerical model for the serpentine PEMFC using COMSOL multi-physics software. To confirm the accuracy of this model, we conducted an actual experiment using the same variables. The study’s main goal is to develop a reliability-based design optimization model for the serpentine PEMFC flow channel, considering the fuel cell’s efficiency as an objective function, cell potential (voltage), and average water activity in the membrane as a constraint function. We initially examined 17 uncertain variables to identify the most significant ones for cell potential and water activity and then employed the surrogate-based Sobol’s sensitivity method. After selecting the significant uncertainty parameter for the respective constraint functions, the study conducted reliability-based design optimization using Monte Carlo simulation. The next section discusses the theory behind the proposed methods.

A surrogate model is an estimation of a high-fidelity model that is computationally accurate. Using a surrogate model, you can estimate the response for a given input without needing to carry out any additional experiments or run a high-fidelity simulation model. In other words, it can be used to estimate the output for a given set of inputs. It is also essential that the surrogate model be computationally effective. Gaussian processes (Rasmussen, 2003; Chen et al., 2023), neural network (Tripathy and Bilionis, 2018), and support vector machines (Xiang et al., 2017) are some of the most common approaches to developing a surrogate model. In this study, a Gaussian Process Regression (GPR) was employed as a surrogate model to replace the computationally expensive numerical model of PEMFC in COMSOL. We also used GPR to determine the feasibility of estimating the propagation of uncertain responses with different sets of uncertain operating parameters.

GPR is a type of supervised machine learning technique that requires only a few parameters to yield a prediction. A Gaussian process is defined as an endless extension of the multivariate normal distribution. The relationship between the input vector and the output parameter can be expressed as: y = f x + ε (28) where f x is the function representing the independent variable for observation, ε represents the noise added to the observed variables. The aim is to use the GPR framework to model an unknown function f x , by assuming that f x at any point x is a Gaussian random variable N μ , σ 2 , where μ and σ are two constants independent of x . For any x , f x is a sample of μ + ε x , where, ε ∼ N 0 , σ n 2 . For any x , x ′ ∈ R d , ψ x , x ′ , the correlation between ε x and ε x ′ , depends on x − x ′ . More precisely ψ x , x ′ = ⁡ exp − ∑ k = 1 d θ k x k − x k ′ p k , (29) where the hyper-parameter p k is in between 1 ≤ p k ≤ 2 which is related to the smoothness of f x to x k , in this study, we took p = 2 and the hyper-parameter θ k should be above zero (0), which indicates the importance of the input variable x k on f x . To estimate the hyper-parameters μ , σ and θ 1 , θ 2 , … , θ k , let’s consider a set of N-number of input points ( x 1 , x 1 , … , x N ) ∈ R d , and their responses will be y = y 1 , y 2 , … , y N T . So, the hyper-parameters can be estimated by maximizing the likelihood function that f x = y k at x = x k , which is given as: L y μ , σ 2 , θ = 1 2 π σ 2 N / 2 Ψ 1 2 exp − y − μ 1 T Ψ − 1 y − μ 1 2 σ 2 , (30) where Ψ ∈ R N × N is an N × N correlation matrix whose elements depend on ψ x i , x j   , and 1 is an N-dimensional column vector of ones. To facilitate the hyper-parameter estimation by the method of maximum likelihood estimation, we can formulate a log-likelihood function as follows: ln ⁡ L y μ , σ 2 , θ = − N 2 ln 2 π − N 2 ln σ 2 − 1 2 ln Ψ − y − μ 1 T Ψ − 1 y − μ 1 2 σ 2 (31)

Out of the three hyper-parameters, we can use analysis to find the maximum likely values for μ and σ 2 by having the first derivatives of ln ⁡ L with respect to μ and σ 2 as follow: ∂ ⁡ ln ⁡ L ∂ μ = 0   ⇒   μ ^ = 1 T Ψ − 1 y N 1 T Ψ − 1 1 (32) ∂ ⁡ ln ⁡ L ∂ σ 2 = 0   ⇒   σ ^ 2 = y N − 1 μ ^ T Ψ − 1 y N − 1 μ ^ N (33)

By substituting Equations 32, 33 into Equation 31, the likelihood function becomes dependent only on hyper-parameter θ ,. Therefore, Equation 31 can be optimized to estimate θ . In this study, MATLAB 2022b global optimization toolbox has been used to optimize the likelihood function.

We can now generate the prediction of y x at a new point x given the maximal likelihood estimates of the hyper-parameters ( μ , σ 2 , θ ). Let r represent the correlation vector between the new point and the N data points, which is the i t h element of r = ψ x , x 1 , … , ψ x , x N T . The optimal linear unbiased estimator for y(x) can be expressed as: y ^ x = μ ^ + r T Ψ − 1 y − 1 μ ^ , (34) and its mean square error can be obtained using the following equation: s 2 = σ ^ 2 1 − r T Ψ − 1 r + 1 − 1 T Ψ − 1 r 2 1 T Ψ − 1 r . (35)

Based on this, the GPR models can predict y x at new points using the observed sample data.

To check that the approximation model is accurate/precise enough, we use numerical measures of the root mean square error (RMSE) and coefficient of determination ( R 2 ). Both the RMSE and R 2 values quantify the degree to which a regression model accurately fits a given dataset. The RMSE quantifies a regression model’s accuracy in predicting the value of a response variable, measured in absolute terms. However, R 2 measures the extent to which the predictor variables can account for the variability observed in the response variable.

The RMSE formulation is given as: R M S E = 1 n ∑ i = 1 n y ^ − y i 2 , (36) where n is the number of sample points, y is the actual value, and y ^ is the predicted one.

The formulation for the R 2 is also expressed as follows: R 2 = 1 − ∑ y ^ − y i 2 ∑ y ¯ − y i 2 , (37) where y ¯ is the mean value of the observed data.

The impact of uncertain input parameters on the cell potential and water content on the membrane is measured by Sobol’ indices. Sobol global sensitivity analysis method is a variance-based sensitivity metrics (Sobol, 1993; Bergamini et al., 2019), which can be easily derived from the surrogate model that casts the input-output relationship. The Sobol’ indices simply represents the ratio of partial variances to the total variance of y . For instance, the Sobol’ index due to the S t h -order interactive effect of the input variables x i 1 , … , x i s is estimated as S i 1 , … i s = V a r i 1 , … i s / V a r , the value of which is always between 0 and 1. With a higher Sobol’ index, the contribution of the associated parameters to the variance of y becomes more important. In particular, the first order, and total Sobol’ indices can be derived from Equations 38, 39 respectively. The first-order effect can be defined as S i = V a r x i E x \ i y | x i V a r y , (38) where x i indicates the i t h input parameter and x \ i represents all parameters except x i . While maintaining x i constant, the inner expectation value indicates that we can compute the mean of y over all potential values of x \ i . Then, the outer variance is calculated over all possible x i values. The first-order Sobol’ index S i measures the influence due to the parameter x i alone, which reflects the marginal effect of x i .

The total effect is also derived as: S T , i = 1 − V a r x \ i E x i y | x \ i V a r y , (39)

The total Sobol’s index, S T , i , summarizes the overall contribution of the input parameter, x i , by taking into account its marginal effect as well as its interactive effects. It should be noted that the more S T , i deviates from S i , the more preponderant the interaction effects among parameters. The second component in the equation can be seen as the first-order influence of all parameters, except the i t h parameter. Thus, subtracting the second term from one yields the contribution of all terms related to x i .

When there are uncertainties, the design optimization problem may involve both probabilistic and deterministic design criteria. Probabilistic criteria are frequently incorporated as constraints to limit the likelihood of failure. However, one can also employ them to establish objectives, like minimizing the cost’s function or maximizing the probability of achieving a specific value.

Given a set of optimization design variables d and a set of random (uncertain) parameters X , the fundamental RBDO formulation can be expressed as follows: m i n i m i z e   f d , μ X s . t . P g i d , X ≤ 0 ≤ P f i h j d = 0 d L ≤ d ≤ d U (40) where, μ X represent the mean value of the uncertain (random) parameters, h j is the j t h equality deterministic constraints, g i d , X represents the i t h non-equality deterministic constraint functions. If the constrain function g i d , X k is less than 0 at the design point X k , the design point X k would be a performance failure point in terms of i t h constraint of g i . Furthermore, P denotes the likelihood of failure, thus P f i means the target likelihood of failure in the i t h constraint.

The probability of failure, P f , can be expressed as follows: P f = P g X 1 , X 1 , … , X n ≤ 0 = ∫ ∫ … ∫ f X 1 , X 1 , … , X n x 1 , x 2 , … , x n d x 1   d x 2 … d x n , g X 1 , X 1 , … , X n ≤ 0 (41) where x 1 , x 2 , … , x n are the uncertain (random) parameters and f X 1 , X 1 , … , X n x 1 , x 2 , … , x n is the joint probability density function. Typically, it is not possible to obtain the analytical solution of Equation 41. Various numerical methods, including Monte Carlo simulation (MCS), importance sampling, and subset simulations, have been developed to calculate the probability of failure. However, this study has adapted the most well-known method, MCS, to determine the failure probability.

The MCS approach provides for the estimation of the probability of failure, which is provided by: p ¯ f = 1 N ∑ i = 1 N I X 1 , X 1 , … , X n , (42) where I X 1 , X 1 , … , X n is a function defined by I X 1 , X 1 , … , X n = 1  if  g X 1 , X 1 , … , X n ≤ 0 0  if  g X 1 , X 1 , … , X n > 0 (43)

According to Equation 41, N independent sets of values X 1 , X 1 , … , X n are obtained by the probability distribution for each random variable and the failure function is calculated for each sample. Using MCS, an estimated probability of the failure is obtained by N-independent sets of values, X 1 , X 1 , … , X n are derived from the probability distribution of each random variable in accordance with Equation 41, and the failure function is calculated for each sample. An estimation of the failure probability is produced by using MCS. p ¯ f = N f N , (44) where N f is the total number of cases where failure has occurred.

The numerical model of the PEMFC was developed in COMSOL 6.1 Multiphysics software, Figure 3 shows a 54 × 54 mm2 three-channel with four repeating units serpentine flow channels type. It is demonstrated that the fuel cell MEA is positioned in a sandwich configuration between two gas diffusion layers, as well as the serpentine flow channels for hydrogen and oxygen. Both the air side and the hydrogen side are placed lower than the MEA. The air side is located above the MEA. In this study, the current distribution in an LT-PEMFC is investigated for serpentine flow field patterns. This is accomplished by operating the cell in counter-flow mode, which ensures that the oxygen and hydrogen inlet flow streams are situated on opposite sides of the cell relative to the membrane’s in-plane direction, as depicted in Figures 3A, B also shows the bipolar palate used for the simulation and testing.

This model solves for the electrode and electrolyte phase potentials in the gas diffusion layers and the membrane, furthermore, it solves for the molar fractions, the gas pressure and the flow velocity on each side of the membrane. In addition, the model incorporates the permeability and electroosmotic drag mechanisms that are responsible for the transport of water through the membrane. We ran a single simulation to demonstrate the numerical model result based on the initial operating parameter value shown in Table 1.

To investigate the uncertain effect and to choose the most significant noise (random) parameters for cell potential and average water activity at the membrane, first we developed a GPR model for each response. To train the model 100 simulations were conducted by applying a Latin hypercubic sampling method (Johnson et al., 1990). For testing 30 simulations were also conducted by generating a random sample dataset. The GPR model was developed for both cell potential and water activity. We calculated the R² and RMSE values to assess the accuracy of the model. The results showed that the cell potential model had an R² of 0.9909 and an RMSE of 0.0045, while the average water activity at the membrane model had an R² of 0.9896 and an RMSE of 0.0134. Figure 7 displays the Actual vs. Prediction plot, which demonstrates a perfect match between the predicted and actual values of both cell potential and water activity. In addition, when we use the numerical simulation, the computing time for a single simulation takes approximately 4 min to yield the cell potential and water activity results. However, the GBR model will return the result in 0.001 s for a single run.

As discussed in Section 3, for sensitivity analysis, the study applied the Sobol method, which is one of the global variance-based sensitivity analysis methods. Figure 8 also displays the computed Sobol indices, demonstrating that parameters such as cell temperature (T), exchange current density at reference oxygen reduction ( i o O 2 r e f ), specific area, catalytic layers ( a C L ), and the oxygen molar fraction in the cathode relative to nitrogen ( O 2 N 2 ), significantly influence cell potential. Cell temperature (T), hydrogen flow stoichiometry ( s t o i c h H 2 ), oxygen flow stoichiometry ( s t o i c h O 2 ), inlet relative humidity at the anode side ( R H a n ), and inlet relative humidity at the cathode side ( R H c a ) also significantly influence the water content at the membrane.

Deterministic optimization has been carried out after the selection of the most significant random (uncertain) variables for the respective target responses. The mean operating parameters are applied in addition to the design variables in deterministic optimization. The PEMFC’s channel structure is then optimized using a multi-objective optimization technique. Genetic algorithms and numerical simulation are combined in the multi-objective optimization method. The height of the channel h c h , the width of the channel w c h , and the width of the ribs w r i b are chosen to be variables of the optimization process. Additionally, the objective is to maximize cell potential output and efficiency while minimizing the average water activity at the membrane. The multi-objective optimization is described by the following mathematical model. min s . t − E c e l l , − η , W a d L ≤ d ≤ d U (45) where E c e l l is cell potential, η is efficiency, W a is the average water activity at the membrane, and d is the design variables. The range of the design variables is also given in Table 3.

Based on the selected random and design variables the surrogate model has been developed for each objective function which is cell potential, average water activity at the membrane and efficiency. Firstly, we developed the GPR surrogate mode for each objective function by generating 100 training sample datasets using the hyper-Latin cube sampling method. Additionally, we generated 30 random datasets for model testing. Figure 9 shows the actual and prediction correlation plots for the testing dataset, along with the R2 value, demonstrating the accurate fitting of all three surrogate models.

As shown in Equation 45, the goal is to maximize cell potential output and efficiency while minimizing the average water activity at the membrane. In this investigation, we employed the multi-objective genetic algorithm to determine the Pareto front among the three objective functions. In this study, we utilize the Pareto fraction and the distance function to regulate the elitism of the genetic algorithm. By giving preference to candidates which are placed relatively far away from the front, the Pareto fraction option and the distance function both help to preserve variety on the front. This is accomplished by giving preference to candidates which are located at a distance. The Pareto fraction option places a limit on the number of candidates which are located on the Pareto front. In this study, we used a population size of 500 individuals, a Pareto percentage of 0.3 (which is thirty percent of the population size) and applied to 100 generations.

Several multiple-criteria decision-making methods have been suggested for different fields to find the best combination to choose from the Parato front result. Such as the analytical hierarchy process (AHP), the technique for order of preference by similarity to the ideal solution (TOPSIS), and the multi-level linguistic decision-making methodology (multi-level LDM). Various studies have also compared these methods from different perspectives (Kolios et al., 2016). This specific study employs the TOPSIS method, which is a widely used methodology in numerous fields that offers faster computational time compared to other methods (Tzeng and Huang, 2011). For the same weight the result shows that, compared to the initial (mean value), the efficiency and cell potential have increased from 0.5510 to 0.5796 and from 0.7273 to 0.7572 V, respectively. Additionally, the average water content at the membrane has decreased from 1.3416 to 1.2872. Table 4 also summarizes the TOPSIS result for different weight of weights of efficiency, cell potential, and average water activity at the membrane.

Following the completion of the sensitivity analysis and deterministic optimization, we performed the RBDO by utilizing the significant random variables for each constraint function. Here we also employed the same surrogate model from the deterministic optimization. The study’s RBDO formulation is also given as follows: m i n i m i z e − f d , μ X s . t . P G W a d , X ≤ G W a d , μ X ≤ P f W a , X = T , i o O 2 r e f , a C L , O 2 / N 2 P G E c e l l d , X ≥ G E c e l l d , μ X ≤ P f E c e l l , X = T , s t o i c h H 2 , s t o i c h O 2 , R H a n , R H c a d L ≤ d ≤ d U X ∼ N μ X , C O V 2 (46) where d denotes the design variable, X is the respective chosen random parameter, μ X is the mean value of the random parameters, f is the objective function (system efficiency), G w c the average water activity function, G c p cell potential function, P f W a and P f E c e l l target probability for average water activity and cell potential, respectively. Each constraint function’s feasibility region, or limit state, is also determined at each iteration. For instance, the cell potential’s limit state( G E c e l l d , μ X ) is obtained from the design variables d and the mean value of the random parameters μ X ; the failure state occurs when it falls less than the limit, while the average water activity’s failure value occurs when it falls above the determined limit value. The study also sets the P f W a and P f E c e l l to 0.05 to target a probability above 95% reliability for both average water activity at the membrane and cell potential. Table 5 also shows the random (uncertain) parameter values, assuming a normal distribution with a covariance of 0.05 for all random operating parameters.

To determine the reliability index for the average water activity in the membrane and the cell potential constraint, we used an MCS method, which involves generating 10,000 samples based on their distribution. For the objective function optimization, the study applied a genetic algorithm. In this investigation, we employed a maximum of 200 generations, a population size of 100, and a function tolerance of 1 e − 6 . Figure 10 depicts the evolution of the fitness value for the objective function (system efficiency) across generations until the genetic algorithm converges to the optimal value. As the result shows, the optimization ends before reaching the specified maximum number of generations. This demonstrates that the specified maximum number of generations is adequate.

Table 6 also summarizes the results for the deterministic and reliable optimal design variable, while Table 7 compares the optimal results of the deterministic optimization and RBDO with the initial values. The result demonstrates that the deterministic optimization yields higher efficiency, however, the reliability for cell potential is 60.92% and 79.31% for water activity at the membrane. The optimal system efficiency values obtained by RBDO are less than those obtained by DO; However, the reliability for cell potential is 95.01% and for water activity at the membrane is 96.85%, which is higher than the deterministic or the initial values. This suggests that the proposed RBDO provides a more reliable and robust design of the serpentine flow channel cross-section.

Figure 11 also shows the histograms of constraints (cell potentials and water activities) to show that the RBDO procedure determines a feasible, reliable design. As shown in the figure, the safe regions of the constraints are denoted by dashed black lines and an arrow, which represent the areas where the design is feasible due to the constraints being met. The histograms also demonstrate that the RBDO approach successfully preserves the system’s level of reliability despite randomness (uncertainty), which may occur at any point in the system’s operation.