Front. Energy Res. Frontiers in Energy Research Front. Energy Res. 2296-598X Frontiers Media S.A. 1206418 10.3389/fenrg.2023.1206418 Energy Research Original Research Multi-objective optimization of PEMFC performance based on grey correlation analysis and response surface method Wu and Luo 10.3389/fenrg.2023.1206418 Wu Gang Luo Na * College of Automotive and Mechanical Engineering, Changsha University of Science and Technology, Changsha, China

Edited by: Xiaohuan Zhao, Jinan University, China

Reviewed by: Qingguo Peng, Guizhou University, China

Yuandong Yang, Shandong University, China

 *Correspondence: Na Luo, LNA4316@163.com 25 05 2023 2023 11 1206418 15 04 2023 16 05 2023 Copyright © 2023 Wu and Luo. 2023 Wu and Luo

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

This paper aims to combine grey correlation analysis and response surface method to propose a fast and effective performance optimization method for PEMFC. First, based on orthogonal test data, grey correlation analysis method is used to select four variables that have significant influence on PEMFC’s comprehensive performance from eight common parameters. Secondly, based on grey correlation analysis, the multi-objective optimization problem is transformed into a single objective optimization problem about correlation degree, and applying the response surface method to build the key parameters and the correlation between the second order prediction model. Therefore, the current density, system efficiency and oxygen distribution uniformity on cathode catalyst layer of PEMFC were optimized as a whole. Finally, the optimal parameter combination was obtained by optimizing the prediction model. The simulation results show that the optimized operating conditions are significantly improved in the three performance indexes compared with the basic model, which confirms the feasibility of this method in solving the multi-objective optimization problem, and can provide some reference for the optimal design of hydrogen fuel cells.

 PEMFC grey correlation analysis response surface method performance optimization multi-objective optimization section-at-acceptance Energy Storage 1 Introduction

Due to the growing demand for energy, the consumption of fossil fuels and the total amount of CO2 emissions are rising rapidly, and “clean, low-carbon, safe and efficient” energy reform has become the general trend (Yang et al., 2023). Fuel cell is a device that directly converts chemical energy in fuel into electrical energy through electrochemical reaction without combustion, which has the advantages of high energy conversion efficiency, no pollution and no noise compared with traditional internal combustion engine, etc. Proton exchange membrane fuel cell (PEMFC) has become one of the hot spots of research in the energy field because of its advantages of low operating temperature, fast starting and high specific power (Kurnia et al., 2019). Many researchers are committed to promoting the development of hydrogen fuel cell technology from the aspects of improving performance, extending life and reducing cost (Parekh, 2022). Performance optimization has always been an important work to promote the development of fuel cells. However, the performance and stability of fuel cells are closely related to various operating parameters. Improper working conditions will seriously affect the output performance of the cell and even cause irreversible degradation inside the cell, thus greatly shortening the cell life. Therefore, it is important to study the influence of operating parameters on cell performance and determine the optimal combination of parameters for efficient and long-life operation of PEMFC.

So far, this has been studied by a number of scholars and some results have been obtained. Dehsara et al. showed that increasing the operating pressure and temperature can improve the net transport of reactants through the porous layer to the catalyst layer (CL) and improve the performance of the proton exchange membrane (Dehsara & Kermani, 2013). Lu et al. divided the pressure increase into the contribution of increasing the partial pressure of reactants and the contribution of increasing the water content of the membrane, further revealing that the improvement of cell performance is mainly due to the increase in membrane conductivity rather than enhancing the partial pressure of the reactants (Lu et al., 2019). Hasheminasab et al. pointed out that increasing the anode stoichiometric would reduce the water content of the cathode channel and the cell performance, and increase the sensitivity of generating power to changes in other operating parameters. However, increasing the cathode stoichiometric and temperature would reduce the water content in the cathode channel and increase the cell power. The sensitivity of cell power to changes in other operating parameters is reduced (Hasheminasab et al., 2020). Zhang et al. indicated that adjusting the relative humidity (RH) is an important way to influence the water management of fuel cells, and that the cells have better overall and local response current density and spatial performance uniformity when the relative humidity is 70% (Zhang et al., 2016). Xu et al. noted that the cell output performance is worst at 90°C with both cathode and anode RH at 90%. Under the condition of high current density, the lower the RH value is, the better the operation performance is, and the influence of the change of cathode RH on the cell output performance is better than that of anode RH (Xu et al., 2021). Xia et al. noted that increasing GDL porosity improves performance by increasing flow uniformity and diffusion flux, but effective electron conductivity decreases with increasing porosity. Since the cell performance depends on the combination of these three factors, the final current density will decrease with further increase in porosity (Lingchao Xia & 2021). Abraham et al. improved the cell performance by 12.5% by optimizing GDL porosity and CL porosity, as cells with optimized porosity increase the membrane water content and promote the diffusion of reactants, ultimately leading to higher current densities (Abraham B & Murugavel K, 2021).

Uneven distribution of reaction gases can lead to uneven current density and temperature distribution, which can reduce the overall cell performance and even cause local reaction gas starvation, thus enhancing carbon corrosion and loss of Pt catalyst layer. Current research on the uniformity of fuel cell reactant distribution focuses on the effects of different flow field structures and current densities. Sauermoser et al. summarized and reviewed many studies in recent years to improve the uniformity of gas distribution in the catalytic layer and the uniformity within GDL by optimizing the flow field structure (Sauermoser et al., 2020). Chen et al. designed a stepped flow field. Compared with the ordinary parallel flow field, the stepped flow field can improve the uniformity of gas concentration and current density distribution, alleviate the water flooding phenomenon, and increase the net power by 21.5% (Chen et al., 2021). Zhang et al. the degradation mechanism and mitigation measures of vehicle fuel cells under the condition of frequent start-stop are reviewed. It is concluded that the local gas starvation caused by the uneven distribution of reactive gas is one of the important reasons for the start-stop degradation of fuel cells (Zhang et al., 2018). Lim et al. summarized and reviewed the influence of flow field design on water management and reactant distribution. They pointed out that the uniform distribution of reaction gas is one of the main criteria for improving performance. In addition to the design of fuel cell pack, operating conditions such as flow rate, pressure, temperature and humidity also affect the distribution of reactants (Lim et al., 2016). Cheng et al. proposed the mass fraction difference coefficient as a uniformity evaluation index, and pointed out that hydrogen humidification has a hydrogen distribution significantly and it is appropriate to control at 25%–50%, while air humidification has less effect on oxygen distribution (Cheng et al., 2021). Zhang et al. found that the concentration distribution of anode gas tends to be more uniform than that of cathode gas, and the concentration gradually increases with the increase of channel distance with a maximum error of 15.6% (Yong et al., 2022). Therefore, it is important to improve the uniform distribution of the cathode reaction gas.

Sheila et al. used a weighting method to optimize the fuel cell system efficiency and size trade-off (Ang et al., 2010). Mehrdad et al. optimized the cell voltage and power density by response surface methodology (RSM) combined with non-dominated sorting genetic algorithm (NSGA ii) using operational parameters as design variables (Ghasabehi et al., 2021). Zeng et al. with flow power consumption and output power as optimization objectives, genetic algorithm (GA) was used to optimize the width of the bottom and top sides of the flow channel, and the optimal structure was trapezoidal channel (Zeng et al., 2017). Liu et al. based on the neural network proxy model and non-dominated sorting genetic algorithm (NSGA ii), the porosity of the catalytic layer (CL), the electrolyte volume fraction (CL) and the gas diffusion layer (GDL) were optimized to improve the output power density and reduce the electrolyte volume fraction (Liu et al., 2022). Liu et al. used a multi-objective genetic algorithm to optimize both the operating conditions and the channel structure of the cell, and the results showed that the optimal channel shape under optimal operating conditions was conical (Liu et al., 2017).

According to the review of previous studies, it is not found that the combination of gray correlation analysis and response surface method to achieve multi-objective optimization of PEMFC, and there is little consideration of the uniformity of reaction gas distribution in performance optimization, and almost no research on the effect of operating parameters on reaction gas distribution. Many scholars often only consider the improvement of current density or power density when performing performance optimization, without considering the influence of gas distribution. However, the uniformity of reaction gas distribution has an important impact on the improvement of cell performance and life, so it is not advisable to improve the output power by losing other performance indicators.

In view of the above deficiencies, this paper first established a three-dimensional CFD model of PEMFC with three-channel serpentine flow fields as the basic optimization model. Secondly, according to a lot of literature reading, eight common influence parameters were selected as the initial decision variables, among which six operating parameters and two diffusion layer structure parameters were included. Several parameters which have great influence on the call comprehensive performance were selected by grey correlation analysis for subsequent optimization. Then the Box-Benhnken design (BBD) is used to design a simulation scheme for the key parameters. Based on the gray correlation theory, the multi-objective optimization about the cell performance is transformed into a single-objective gray correlation degree optimization problem, and the response surface method is applied to construct a second-order regression model between the key parameters and the gray correlation degree. Finally, the optimal parameter combination is obtained by seeking the solution that satisfies the maximum gray correlation degree, and the feasibility of the optimization results is verified by simulation.

 2 PEMFC model development 2.1 Model description and assumptions

In this paper, a three-dimensional geometric model of PEMFC was established based on solidworks, as shown in Figure 1. The PEMFC consists of nine parts: anode and cathode collectors, flow channels, gas diffusion layers (GDLs), catalyst layers (CLs) and proton exchange membrane. The flow field type is three-channel serpentine, and the specific geometric parameters of this model are shown in Table 1.

Geometric model of PEMFC.

Main geometric parameters of PEMFC model.

Parameter Value
Collector thickness 1.5 mm
GDL thickness 0.23 mm
CL thickness 0.01 mm
Membrane thickness 0.025 mm
Channel width 1 mm
Channel height 0.5 mm
Rib width 1 mm
Activated area 73 × 73 mm

Due to the complexity of the actual operating process of PEMFC. On the premise of fully describing the main reaction process of PEMFC, several assumptions were made in the numerical simulations in order to simplify the model and calculations as follows (Li et al., 2017; Chen et al., 2019; Liu et al., 2022; Yao et al., 2022).

1) PEMFC operation state is steady state;

2) The initial temperature of the fuel cell is constant, and the thermal effect caused by the temperature gradient is not considered;

3) The reactant gases are ideal and incompressible, and the gas flow in the flow channel is laminar;

4) The gas diffusion layers, catalyst layers and membrane are homogeneous models with isotropic porous media;

5) Only protons can permeate the proton exchange membrane, but not the reaction gases.

2.2 Governing equations 2.2.1 Mass conservation equation

∂ ε ρ ∂ t + ∇ ⋅ ε ρ u → = S m ρ , ε , u → , S m are density, porosity, velocity vector, the mass source term in that order. For both cathode and anode flow channels and GDLs, S m is 0. The mass source term for the anode and cathode CLs can be expressed as follows: S m = S H 2 = − M H 2 2 F i a S m = S H 2 O + S O 2 = M H 2 O 2 F i c − M O 2 4 F i c Where M and F represent molar mass and Faraday’s constant, respectively, i a and i c represent the exchange current density of the anode and cathode, respectively, which can be calculated by Eqs 10, 11.

 2.2.2 Momentum conservation equation

∂ ε ρ u → ∂ t + ∇ ⋅ ε ρ u → u → = − ε ∇ p + ∇ ⋅ ε μ ∇ u → + S u p , μ , S u are pressure, viscosity, and the momentum source term, respectively. Assuming that the source term in the flow channels and membrane is zero, the momentum equation in the porous electrode can be expressed as follows: S u = − ε μ K u →

Where K is the absolute permeability.

 2.2.3 Energy conservation equation

∂ ε ρ c p T ∂ t + ∇ ⋅ ε ρ c p u → T = ∇ ⋅ k e f f ∇ T + S Q c p , k e f f , T are the specific heat capacity at constant pressure, effective thermal conductivity and temperature, respectively. S Q is the energy source term, including ohmic heat, chemical reaction heat, gas-liquid phase transition heat and heat generated by over-potential. It can be expressed as follows: S Q = I 2 R o h m + β S H 2 O h r e a c t i o n + r W h l + i a , c η

Where I is current, R o h m is resistance, β is the ratio of chemical energy to heat energy, S H 2 O is liquid water generation rate, h r e a c t i o n is the enthalpy of the reaction, r W is the phase transition rate of water, h l is the enthalpy of the phase transition of water, η is the overpotential.

 2.2.4 Species conservation equation

∂ ε c k ∂ t + ∇ ⋅ ε u → c k = ∇ ⋅ D k e f f ∇ c k + S k c k , D k e f f , S k are concentration, the effective diffusion coefficient and the species source term, respectively, and the subscript k represents each species. The species source term is 0 in the flow channels and GDLs, and the source term of hydrogen, oxygen and water in the CLs can be expressed as: S H 2 = − M H 2 / 2 F i a , S O 2 = − M O 2 / 4 F i c , S H 2 O = M H 2 O / 2 F i c .

 2.2.5 Current conservation equation

∇ ⋅ σ e ∇ ∅ e + S e = 0 ∇ ⋅ σ m ∇ ∅ m + S m = 0 ∅ e and ∅ m are the solid phase potential and the membrane phase potential, respectively, S e and S m are the electron current source term and the proton current source term, respectively, and S e + S m = 0, σ is conductivity.

The Butler-Volmer equations of anode and cathode in the CLs can be respectively expressed as follows: S a = j a , r e f C H 2 C H 2 , r e f γ a e α a F R T η a − e − α c F R T η c S c = j c , r e f C O 2 C O 2 , r e f γ c e α a F R T η c − e − α c F R T η a

Where η is overpotential, j r e f is the reference exchange current density, C i is the molar concentration of species, C i , r e f is the reference molar concentration of species, γ is the concentration index, α is the exchange coefficient.

 2.3 Model verification 2.3.1 Grid independence test

Grid independence test can check whether the simulation results will be affected by the number of grids, so as to determine the appropriate number of grids. This paper uses the grid division software ICEM to divide the PEMFC model with four different numbers of grids, and then the four models are simulated and calculated under the same boundary conditions. The error between simulation results and experimental results is shown in Figure 2. When the number of grids is 3,470,688 and 5,958,942, the error is 2.4% and below. Considering the time cost and simulation accuracy, 3,470,688 mesh quantities are selected, and the whole calculation domain is shown in Figure 3.

Simulation accuracy of different mesh numbers.

(A) Diagram of the overall computing domain (B) Partial diagram of the entrance.

 2.3.2 Experimental verification

The mesh file of the PEMFC model is imported into ANYSY FLUENT 2020R2 and the add-on PEMFC module is called to solve the model. The relevant parameters in the model were set as shown in Table 2. The type of inlet boundary was mass-flow inlet, and the type of outlet boundary was pressure-out. The inlet mass flow rate was calculated according to the anode and cathode stoichiometric ratio of 2, relative humidity of 100%, temperature of 80 °C and operating pressure of 1atm.The output current density under different working voltages was calculated under constant voltage. The structure and numerical simulation conditions of the PEMFC used in this paper are basically consistent with the experimental conditions of Li et al. (Li et al., 2017), so the model in this paper is verified by the experimental data of Li et al. The comparison between simulation results and experimental results is shown in Figure 4. At low to medium current densities, the numerical calculation results agree well with the experimental data. When the current density increases to 0.9A/cm2, the simulation results begin to deviate from the experimental results. And the simulation results are always slightly larger than the experimental results. This is because this paper uses a single-phase model for numerical simulation, and the generated water exists in gaseous form by default. In fact, when the current density is larger, the liquid water tends to accumulate in the flow channel, and the liquid water tends to block the pores in the gas diffusion layer and catalyst layer, thus hindering oxygen transport and eventually leading to performance degradation at high current densities. Overall, the numerical simulations are consistent with the trend of experimental results, and the PEMFC model used can be considered effective.

Parameter Settings of PEMFC numerical model.

Item Value Item Value
Anode References exchange Current density/(A.m-2) 10,000 Cathode References exchange Current density/(A.m-2) 17
H2 References concentration/(kmol.m3) 0.005 O2 References concentration/(kmol.m3) 0.01
Anode concentration exponent 0.5 Cathode concentration exponent 1
Anode exchange coefficient 0.5 Cathode exchange coefficient 1
H2 diffusion coefficient/(m2.s-1) 1.1*10–4 O2 diffusion coefficient/(m2.s-1) 3.2*10–5
H2O diffusion coefficient/(m2.s-1) 7.35*10–5 Other diffusion coefficient/(m2.s-1) 1.1*10–5
GDL porosity 0.78 CL porosity 0.5
CL surface/volume ratio/m-1 200,000 Membrane equivalent weight/(kg.kmol-1) 1,100
Membrane electrical conductivity/ohm.m 1*10–16 Collector electrical conductivity/ohm.m 1000000
Protonic conduction exponent 1 Open-circuit voltage/V 0.92

Polarization curves of model simulation and experimental data.

 3 Research method 3.1 Analysis of significance

In previous parametric optimization studies, people usually select three to six optimization parameters based on experience, without stating the reasons for the selection, so it is somewhat subjective. In this study, we use gray correlation analysis to analyze the significance of the variables before optimization, and select the variables that have a greater impact on the optimization target, which not only can improve the optimization efficiency, but also can get more accurate results.

Gray correlation analysis, a systematic scientific theory first proposed by Jurong Deng in 1982, is feasible for analyzing the degree of influence of multiple variable factors on a specific quality indicator by analyzing the degree of similarity between the geometry of the reference sequence and the comparative sequence to determine the closeness of the association (Yin et al., 2015). Since the units and ranges of each series are different, the original data need to be normalized by dimensionless preprocessing before performing gray correlation analysis. The normalization process can be performed with the following two equations depending on the type of quality characteristics (Lian et al., 2022). If the actual data series has the quality characteristic of “bigger is better”, then equation (13) is used for preprocessing. If the actual data series has the quality characteristic of “smaller is better”, then equation (14) is used for preprocessing (Xu et al., 2020). x i k ′ = x i k − min x i k max   x i k − min x i k x i k ′ = max x i k − x i k max x i k − min x i k

Where x i k is the actual data, min x i k is the minimum value in the ith series, max   x i k is the maximum value in the ith series, and x i k ′ is the pre-processed comparison sequence, where k = 1, 2, 3 … … 27; i = 1, 2, 3.

After normalizing the test results, the gray correlation coefficient ξ corresponding to each comparison sequence and the reference sequence (all response reference sequence are 1) can be calculated by (Eq. 15). ξ x i k ′ , x 0 k = ∆ min + ς ∆ max ∆ i k + ς ∆ max

Where ς is the resolution coefficient, ς ϵ 0,1 , which is generally taken as 0.5, x 0 k is the reference sequence, ∆ i k is the deviation sequence, ∆ i k = x 0 k − x i k ′ , ∆ min is the minimum deviation, and ∆ max is the maximum deviation.

The gray correlation degree is a weighted sum of gray correlation coefficients, which can be calculated by (Eq. 16), and here the weight of all three indicators is considered to be 0.3333. λ = ∑ k = 1 n β k ξ

Where β k is the weight, ∑ k = 1 n β k = 1 , and n is the number of response targets, n = 3.

Eight common operating parameters and gas diffusion layer structure parameters, namely, temperature, pressure, anode and cathode stoichiometry, anode and cathode relative humidity, anode and cathode gas diffusion layer porosity, were selected for significance analysis to investigate the extent of their effects on three performance indexes, namely, current density, oxygen distribution uniformity and system efficiency of the PEMFC. In order to reduce the calculation time and improve the calculation efficiency, eight factor and three level tests were designed based on the orthogonal design method. The setting of factor levels was shown in Table 3. The three levels of these eight parameters were selected based on the general range and experience of PEMFC parameters (Öztürk et al., 2017; Abraham B & Murugavel K, 2021), and the specific test arrangement is shown in Table 4.

Parameter settings of orthogonal test.

Factors Abbreviation Level 1 Level 2 Level 3
Temperature (°C) T 70 80 90
Pressure (atm) p 1 2 3
Stoichiometric Ratio of Anode ζ a 1 3 5
Stoichiometric Ratio of Cathode ζ c 1 3 5
Relative Humidity of Anode RHa 30% 60% 100%
Relative Humidity of Cathode RHc 30% 60% 100%
Porosity of anodic GDL ε a 0.3 0.5 0.8
Porosity of cathode GDL ε c 0.3 0.5 0.8

L27 (38) Orthogonal test table.

Number T P ζ a ζ c RHa (%) RHc (%) ε a ε c
1 70 3 1 1 100 60 0.8 0.5
2 70 2 1 1 60 100 0.5 0.8
3 70 3 5 3 30 60 0.5 0.3
4 70 2 5 3 100 100 0.3 0.5
5 70 1 3 5 100 30 0.5 0.5
6 70 2 3 5 30 100 0.8 0.3
7 70 1 1 1 30 30 0.3 0.3
8 70 1 5 3 60 30 0.8 0.8
9 70 3 3 5 60 60 0.3 0.8
10 80 2 5 5 30 60 0.5 0.8
11 80 1 5 5 100 100 0.3 0.3
12 80 1 3 1 30 100 0.8 0.8
13 80 2 3 1 60 60 0.3 0.5
14 80 3 5 5 60 30 0.8 0.5
15 80 2 1 3 100 60 0.8 0.3
16 80 3 1 3 30 30 0.3 0.8
17 80 3 3 1 100 30 0.5 0.3
18 80 1 1 3 60 100 0.5 0.5
19 90 2 3 3 100 30 0.5 0.8
20 90 3 3 3 30 100 0.8 0.5
21 90 2 1 5 30 30 0.3 0.5
22 90 1 5 1 30 60 0.5 0.5
23 90 1 3 3 60 60 0.3 0.3
24 90 3 5 1 100 100 0.3 0.8
25 90 2 5 1 60 30 0.8 0.3
26 90 1 1 5 100 60 0.8 0.8
27 90 3 1 5 60 100 0.5 0.3
3.2 Multi-objective optimization method

Response surface method is a mathematical statistical method commonly used to solve multivariate problems, which can construct a functional relationship between the target and the variables based on the existing measured data, and further solving can obtain the best combination of parameters to satisfy the conditions. Before response surface analysis, the central combination design or BBD method is usually used to design the simulation cases and simulate to get certain data. In this paper, we use BBD method to design simulation cases as shown in Table 7. This method not only can avoid extreme conditions at the test point, but also can perform continuous analysis of the test at all levels. According to the results of significance analysis, the four key influencing parameters of cathode and anode stoichiometry, pressure and cathode relative humidity are used as independent variables, and the optimization objectives are current density, system efficiency and oxygen distribution uniformity. Current density could be directly obtained through simulation calculation, the system efficiency can be calculated by the following formula: η s y s = W − W p r s W f u e l

Where W is the output power of PEMFC, W p r s represents the parasitic power, and W f u e l is the inherent power of the fuel. They are expressed as follows: W = w A W p r s = W c o m p + W o t h e r s W c o m p = c p T e η c η m t P P i n 0.286 − 1 m a i r m a i r = 3.57 × 10 − 7 ζ c i A W o t h e r s = 0.05 W W f u e l = ζ a i A 2 F L H V

Where W c o m p is the power consumption of the compressor, and W o t h e r s is the other power losses. Some parameters in the above equations are assumed to be constant, but in practice they may vary with compressor size and full load ratio. However, this does not affect the optimization process and a similar approach has been used in other literature (Na & Gou, 2007; Li et al., 2021). The values of these constants can be found in the literature (Ang et al., 2010).

The oxygen distribution uniformity is expressed using the standard deviation of the molar concentration of oxygen at the interface between the cathodic gas diffusion layer and the catalyst layer, and its expression can be defined as (Chen et al., 2019): s d C = n ∑ i = 1 n C i 2 − ∑ i = 1 n C i 2 n 2 Where C i represents the oxygen Mole concentration value at each point on the catalyst layer.

The sub-objectives of multi-objective optimization problems are often contradictory to each other. It is impossible to make multiple sub-objectives reach the optimal value together at the same time, only coordination and compromise between them, and the optimization solution is not unique. If a final solution is to be sought. The relative importance degree of each objective, namely, the weight, can be decided in advance. Based on this, the algorithm converts the multi-objective problem into a single-objective problem for solving. For different practical applications, each objective has different weights. In this study, the three indicators are considered equally important. The multiple responses of current density, system efficiency and oxygen distribution uniformity are first transformed into a single objective gray correlation based on gray correlation analysis. Then a second-order prediction model between the Key influencing parameters and the gray correlation is constructed, and the final optimized solution is obtained by solving for the solution that satisfies the maximum gray correlation.

 4 Results and discussion 4.1 Key influence parameter selection

The results of the three performance indexes are listed in Table 5 by simulating 27 different sets of working conditions in the orthogonal design. Then according to Eqs 13–16, the grey correlation analysis is carried out on the statistical results, and the grey correlation coefficient and correlation degree are calculated, as shown in Table 5.

Results of orthogonal test and grey correlation analysis.

Test results Standardized Grey correlation coefficient Gray correlation degree
i σ o 2 η / % i σ o 2 η / % ξ 1 ξ 2 ξ 3
1 0.86551 0.0043189 47.70 0.38142 0.00000 0.90740 0.4470 0.3333 0.8437 0.541,303
2 0.80581 0.0027232 49.42 0.32206 0.44526 0.94497 0.4245 0.4740 0.9008 0.599,728
3 1.29078 0.0033527 7.84 0.80429 0.26960 0.03691 0.7187 0.4064 0.3417 0.488,888
4 1.15217 0.0018796 8.88 0.66647 0.68066 0.05962 0.5999 0.6103 0.3471 0.519,028
5 1.00649 0.0019165 17.31 0.52161 0.67034 0.24372 0.5110 0.6027 0.3980 0.50385
6 1.21748 0.0021428 13.12 0.73141 0.60721 0.15222 0.6505 0.5600 0.3710 0.527,135
7 0.57948 0.0014097 51.94 0.09701 0.81176 1.00000 0.3564 0.7265 1.0000 0.69422
8 0.74750 0.0013269 10.39 0.26408 0.83487 0.09260 0.4046 0.7517 0.3553 0.503,802
9 1.48760 0.0026747 10.25 1.00000 0.45879 0.08954 1.0000 0.4802 0.3545 0.611,507
10 1.09434 0.0018864 7.87 0.60896 0.67875 0.03756 0.5611 0.6088 0.3419 0.503,904
11 0.79032 0.0017373 10.39 0.30666 0.72037 0.09260 0.4190 0.6413 0.3553 0.471,815
12 0.55498 0.0008491 17.31 0.07265 0.96821 0.24372 0.3503 0.9402 0.3980 0.562,782
13 0.77043 0.0025771 16.47 0.28688 0.48604 0.22538 0.4122 0.4931 0.3923 0.432,473
14 1.33319 0.0029804 6.15 0.84647 0.37349 0.00000 0.7651 0.4439 0.3333 0.514,035
15 0.99930 0.0016224 44.38 0.51446 0.75243 0.83490 0.5073 0.6688 0.7518 0.64258
16 0.99318 0.0025716 39.22 0.50837 0.48757 0.72221 0.5042 0.4939 0.6428 0.546,921
17 0.84358 0.0039184 15.90 0.35961 0.11175 0.21293 0.4384 0.3602 0.3885 0.395,658
18 0.74312 0.0013123 51.94 0.25973 0.83896 1.00000 0.4031 0.7564 1.0000 0.719,768
19 1.10602 0.0019825 14.79 0.62058 0.65193 0.18869 0.5686 0.5896 0.3813 0.513,092
20 1.23882 0.0024318 13.07 0.75262 0.52656 0.15112 0.6690 0.5136 0.3707 0.517,725
21 0.99754 0.0018864 39.35 0.51271 0.67875 0.72505 0.5064 0.6088 0.6452 0.586,764
22 0.48192 0.0007351 10.39 0.00000 1.00000 0.09260 0.3333 1.0000 0.3553 0.56281
23 0.73065 0.0010655 17.31 0.24732 0.90783 0.24372 0.3991 0.8443 0.3980 0.54711
24 0.83747 0.0034450 9.54 0.35354 0.24385 0.07403 0.4361 0.3980 0.3506 0.394,896
25 0.71375 0.0025214 9.88 0.23052 0.50158 0.08146 0.3939 0.5008 0.3525 0.415,668
26 0.83387 0.0017981 51.94 0.34996 0.70340 1.00000 0.4348 0.6277 1.0000 0.687,409
27 1.38757 0.0024134 30.74 0.90053 0.53171 0.53702 0.8341 0.5164 0.5192 0.62316

According to the grey correlation analysis theory, the greater the grey correlation degree, the closer the corresponding target sequence is to the optimal (Guo et al., 2019). Figure 5 shows the average grey correlation degree of each parameter level. It can be seen from Figure 5 that the parameter combination with the maximum grey correlation degree, that is, the parameter combination close to the better cell performance, is:T = 70,p = 1, ζ a = 1, ζ c = 5,RHa = 30%,RHc = 60%, ε a = 0.5, ε c = 0.8. In addition, the range of the average grey correlation degree of the influencing parameters reflects the influence of the parameters on the multi-objective response. The larger the difference, the more significant the influence of the independent variable on the dependent variable (Wang & Li, 2021). It can be seen from Figure 5 that the comprehensive influence degree of eight influencing parameters on current density, oxygen distribution uniformity and system efficiency is ranked as: ζ a (0.140,776)> P (0.06883) > ζ c (0.04778)> RHc(0.03822)> RHa (0.03572)> ε a (0.0268)> T (0.02216)> ε c (0.01309). Therefore, in order to improve the optimization efficiency, four key parameters, ζ a , ζ c ,P and RHc, which have a large impact on the comprehensive performance of the cell, are selected for the next optimization study. While the remaining four variables, which have a smaller impact, are fixed at a better level and kept constant:T = 70 °C, RHa = 30%, ε a = 0.5, ε c = 0.8.

Average grey correlation degree of each parameter level.

 4.2 Optimization model building and verification

In this section, the optimization model between the key parameters and the response target is mainly combined with the gray correlation analysis method and the response surface method, and the model is validated. In order to reasonably select the test site, BBD was used to design the four-factor and three-level test arrangement, as shown in Table 6. According to the test arrangement, the three performance indicators were numerically simulated and calculated, and then the grey correlation degree was calculated by using the above grey correlation analysis method. The results were collected in Table 6. In order to obtain the optimal parameter combination, it is necessary to accurately construct the mapping relationship between the grey correlation degree and the key influencing parameters. Therefore, Design-expert software is applied to construct the second-order regression model between the response value of the grey correlation degree and the input variables by using the response surface method. The response function is shown in Eq. 25. λ = 1.098 − 0.162898 ζ a + 0.007019 ζ c − 0.330919 P + 0.001188   R H c + 0.002971 ζ a × ζ c + 0.015976 ζ a × P + 0.037622 ζ c × P + 0.000328 P × R H c + 0.01407 ζ a 2 − 0.011274 ζ c 2 + 0.02861 P 2 − 9.57546 E − 6 R H c 2

BBD test arrangement and grey correlation analysis results.

Number ζ a ζ c P RHc i /A.cm -2 σ o 2 η / % λ
1 1 3 2 30 0.996,732 0.002081441 44.38 0.629,784
2 1 1 2 65 0.802,702 0.002871944 49.42 0.601,121
3 3 5 2 30 1.018406 0.0021159314 13.12 0.503,459
4 1 5 2 65 0.991,851 0.001567094 39.35 0.646,323
5 3 3 2 65 1.120,158 0.00202,838 14.79 0.535,993
6 3 3 1 30 0.662,398 0.001272287 17.31 0.565,158
7 3 5 2 100 1.234,363 0.001996927 13.12 0.567,831
8 5 3 3 65 1.319,101 0.0030638101 7.84 0.523,428
9 3 1 1 65 0.603,492 0.001193815 17.31 0.574,062
10 3 1 3 65 0.866,141 0.0047093979 15.9 0.377,269
11 1 3 1 65 0.889,993 0.001354796 51.94 0.781,628
12 5 3 2 30 0.921,569 0.001980908 8.88 0.490,753
13 1 3 3 65 0.99359 0.0019819768 39.22 0.601,039
14 5 3 1 65 0.813,987 0.001223425 10.39 0.576,206
15 3 3 2 65 1.112,436 0.002112043 14.79 0.527,077
16 3 3 3 100 1.345,636 0.0029256631 13.07 0.55163
17 5 3 2 100 1.111,258 0.001865607 8.88 0.537,736
18 3 3 1 100 0.828,013 0.001341491 17.31 0.571,762
19 3 3 3 30 1.265,085 0.0033588643 13.07 0.499,071
20 3 5 1 65 0.951,094 0.001862389 17.31 0.523,439
21 1 3 2 100 0.99468 0.001526076 44.38 0.688,264
22 5 1 2 65 0.78824 0.002924677 9.88 0.412,334
23 3 1 2 30 0.771,814 0.0032176832 16.47 0.410,546
24 3 3 2 65 1.130,049 0.0019544 14.79 0.545,057
25 3 3 2 65 1.130,207 0.00195,102 14.79 0.545,408
26 3 5 3 65 1.45988 0.0027087119 10.25 0.627,619
27 5 5 2 65 1.220,229 0.002686816 7.87 0.505,073
28 3 1 2 100 0.798,285 0.002786262 16.47 0.433,602
29 3 3 2 65 1.112,479 0.00212,392 14.79 0.526,119

Applying the resulting regression model for data prediction, the results of comparing the experimental and predicted values of gray correlation degree are shown in Figure 6. As shown in the figure, the model predicted values and calculated values are in high agreement, and the average relative error is only 2.23%. Figure 7 is the residual diagram of the prediction model. All data points in the figure are randomly distributed around the 0 value, indicating that the predicted value has a good fit with the sample value. Meanwhile, the ANOVA results of the regression equation in Table 7 show that p < 0.0001, indicating that the prediction model is very significant. R-Sq indicates the fit between the established model and the experimental data, R-Sq of the established model is 96.7%, indicating that the prediction model fits the experimental data very well. R-Sq (adj) is the adjusted R-Sq. R-Sq (adj) = 94.23% is very close to R-Sq, indicating that the prediction model is very reliable. In summary, the regression model established based on gray correlation analysis has a good fitting goodness and predictive ability, and can be solved by the next step of parameter optimization.

Comparison of model predicted value and test value of grey correlation degree.

Residual diagram of grey relational degree model.

Results of variance analysis of regression model.

Source Degrees of freedom Sum of squares Mean square F Value p-Value
Model 12 0.1903 0.0159 39.09 <0.0001
Residual 16 0.0065 0.0004
Cor total 28 0.1968
R-Sq = 96.7% R-Sq (adj) = 94.23%
4.3 Optimization results and verification

The influence of key parameters on different targets is reflected by the change of grey correlation degree, so the overall optimization of different targets can be achieved by optimizing the grey correlation degree. The response surface optimizer in Design-expert software is used to optimize and solve the established second-order prediction model within the range of each parameter, and the optimal solution is λ = 0.773. The corresponding parameter combination is as follows: ζ a = 1, ζ c = 2.113, p = 1, RHc = 79%. The optimized operating conditions were numerically simulated and compared with the performance of the basic model. The comparison results are shown in Table 8. Compared with before the optimization, the current density increased by 3.8%, the system efficiency increased by 25.97%, and the standard deviation of oxygen concentration was almost unchanged.

Performance comparison between the optimization results and the base model.

ζ a ζ c P RHc (%) i σ o 2 η
Base model 2 2 1 100 0.762,076 0.001150605 25.97%
Optimization results 1 2.113 1 79 0.792,395 0.001156882 51.94%

The polarization curve comparison before and after optimization is shown in Figure 8 It can be clearly seen from the figure that the output current density of the optimized model is higher than that of the basic model in the whole voltage range. Performance improvement is more obvious under low pressure conditions, with better polarization performance. Under the condition of high current density, obvious concentration polarization appears in the basic model, but this phenomenon is not obvious in the optimized model, mainly because the optimized model has a higher concentration of reactants in the catalytic layer, resulting in less concentration polarization loss.

Comparison of polarization curves between the basic model and the optimized model.

Figure 9 shows the current flux nephogram of the cathode membrane interface before and after optimization. It can be seen from the figure that the optimized model has a higher current density. This in turn confirms the results of Figure 8. Figure 10 represents the cloud plot of oxygen distribution at the cathode GDL-CL interface of the two models. From the figure, it can be seen that the oxygen distribution uniformity of the optimized model is almost the same as that of the basic model, but the molar concentration of oxygen in the optimized model is higher than that of the basic model. This is mainly because the cathode inlet air stoichiometric increases, and more reactants are added to the cell. Moreover, the mass fraction of oxygen increases with the decrease of relative humidity, resulting in more oxygen being transported to the catalyst layer. Thus enabling the optimized model to generate more current, which is also consistent with Figure 9. On the basis of ensuring the uniformity of oxygen distribution on CL, increasing the molar concentration of oxygen can not only improve the output performance of the cell, but also alleviate the aging of PEMFC. The system efficiency increases because the anode stoichiometric decreases while the output power increases, which can be deduced from equations 17–23.

Current density distribution nephogram at the m-cl interface of cathode: (A) Basic model; (B) The optimization model.

Oxygen distribution nephogram at the gdl-cl interface of cathode (A) Basic model; (B) The optimization model.

 5 Conclusion

In this paper, a multi-objective performance optimization design method of PEMFC based on gray correlation analysis and response surface method is proposed with current density, oxygen distribution uniformity and system efficiency as optimization objectives. Firstly, the parameters that have significant influence on the comprehensive objectives are selected as decision variables by using gray correlation analysis. Then the influence of key parameters on different objectives is reflected as a whole through the gray correlation degree by assigning weights. The second-order model between the key influence parameters and the gray correlation degree is constructed by applying the response surface method, so as to seek the optimal combination of parameters. The following conclusions are drawn.

1) The influence of the eight influencing parameters considered in this paper on the comprehensive performance of PEMFC is ranked as: ζ a (0.140,776)> P (0.06883) > ζ c (0.04778)> RHc(0.03822)> RHa (0.03572)> ε a (0.0268)> T (0.02216)> ε c (0.01309). Among them, anode and cathode stoichiometric, pressure and cathode relative humidity are the main parameters that affect the optimization objective. Selecting these four key parameters as control variables can reduce the space of optimization design.

2) The established gray correlation prediction model shows by ANOVA that its prediction accuracy is high and can better describe the response of the gray correlation degree on the design variables.

3) The optimal combination of operating parameters was obtained as: ζ a = 1, ζ c = 2.113,P = 1,RHc = 79%. The overall performance of the PEMFC was verified by numerical simulation to be much improved compared with that before optimization, in which the current density was improved by 3.8%, the system efficiency was improved by 25.97%, and the oxygen distribution uniformity remained basically the same.

4) The optimization method combining grey correlation analysis and response surface method can transform multi-objective optimization into single-objective optimization, and quickly find a suitable unique optimal solution. It can provide a better reference for the optimization of PEMFC parameters because the multi-objective optimization solution is unique and the overall performance is improved compared with previous studies.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

 Author contributions

GW provided research ideas and participated in the revision and finalization of the paper. NL participated in the entire research process and the writing and revision of the first draft. All authors contributed to the article and approved the submitted version.

This work is supported by National Natural Science Foundation of China (Grant No. 51976016) and supported by the Natural Science Foundation of Hunan Province, China (Grant No. 2023JJ10043).

 Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

 Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

 References Abraham B. P. Murugavel K. K. (2021). Influence of catalyst layer and gas diffusion layer porosity in proton exchange membrane fuel cell performance. Electrochimica Acta 389, 138793–138806. 10.1016/j.electacta.2021.138793 Ang S. M. C. Brett D. J. L. Fraga E. S. (2010). A multi-objective optimisation model for a general polymer electrolyte membrane fuel cell system. J. Power Sources, 195 (9), 2754–2763. 10.1016/j.jpowsour.2009.10.095 Chen H. Liu B. Zhang T. Pei P. (2019). Influencing sensitivities of critical operating parameters on PEMFC output performance and gas distribution quality under different electrical load conditions. Appl. Energy 255, 113849–113862. 10.1016/j.apenergy.2019.113849 Chen X. Chen Y. Liu Q. Xu J. Liu Q. Li W. (2021). Performance study on a stepped flow field design for bipolar plate in PEMFC. Energy Rep., 7, 336–347. 10.1016/j.egyr.2021.01.003 Cheng Z. Luo L. Huang B. Jian Q. (2021). Effect of humidification on distribution and uniformity of reactants and water content in PEMFC. Int. J. Hydrogen Energy 46 (52), 26560–26574. 10.1016/j.ijhydene.2021.05.129 Dehsara M. Kermani M. J. (2013). The effects of operating parameters on the performance of proton exchange membrane fuel cells. Mechanics 19 (6), 649–656. 10.5755/j01.mech.19.6.5989 Ghasabehi M. Shams M. Kanani H. (2021). Multi-objective optimization of operating conditions of an enhanced parallel flow field proton exchange membrane fuel cell. Energy Convers. Manag. 230, 113798–113819. 10.1016/j.enconman.2020.113798 Guo J. Shi Y. Chen Z. Yu T. Zhao P. Shirinzadeh B. (2019). Optimal parameter selection in robotic belt polishing for aeroengine blade based on GRA-RSM method. Symmetry 11 (12), 1526–1538. 10.3390/sym11121526 Hasheminasab M. Kermani M. J. Nourazar S. S. Khodsiani M. H. (2020). A novel experimental based statistical study for water management in proton exchange membrane fuel cells. Appl. Energy 264, 114713–114729. 10.1016/j.apenergy.2020.114713 Kurnia J. C. Sasmito A. P. Shamim T. (2019). Advances in proton exchange membrane fuel cell with dead-end anode operation: A review. Appl. Energy 252, 113416–113434. 10.1016/j.apenergy.2019.113416 Li H. Xu B. Lu G. Du C. Huang N. (2021). Multi-objective optimization of PEM fuel cell by coupled significant variables recognition, surrogate models and a multi-objective genetic algorithm. Energy Convers. Manag. 236, 114063–114075. 10.1016/j.enconman.2021.114063 Li W. Zhang Q. Wang C. Yan X. Shen S. Xia G. (2017). Experimental and numerical analysis of a three-dimensional flow field for PEMFCs. Appl. Energy 195, 278–288. 10.1016/j.apenergy.2017.03.008 Lian G. Liu Z. Zhang Y. Que L. Chen C. Yue K. (2022). Grey relational analysis–based composite coating property optimization fabricated by laser cladding. Int. J. Adv. Manuf. Technol. 120, 7985–7996. 10.1007/s00170-022-09305-x Lim B. H. Majlan E. H. Daud W. R. W. Husaini T. Rosli M. I. (2016). Effects of flow field design on water management and reactant distribution in PEMFC: A review. Ionics 22 (3), 301–316. 10.1007/s11581-016-1644-y Lingchao Xia M. N. Xu Q. Cheng C. Cheng C. (2021). Optimization of gas diffusion layer in high temperature PEMFC with the focuses on thickness and porosity. Appl. Energy 300, 117357–117369. 10.1016/j.apenergy.2021.117357 Liu S. Chen C. Tan J. Hu H. Xuan D. (2022). Multi‐objective optimization of porous layers for proton exchange membrane fuel cells based on neural network surrogate model. Int. J. Energy Res. 46 (14), 19796–19813. 10.1002/er.8503 Liu Z. Zeng X. Ge Y. Shen J. Liu W. (2017). Multi-objective optimization of operating conditions and channel structure for a proton exchange membrane fuel cell. Int. J. Heat Mass Transf. 111, 289–298. 10.1016/j.ijheatmasstransfer.2017.03.120 Lu J. B. Wei G. H. Zhu F. J. Yan X. H. Zhang J. L. (2019). Pressure effect on the PEMFC performance. Fuel Cells 19 (3), 211–220. 10.1002/fuce.201800135 Na W. Gou B. (2007). The efficient and economic design of PEM fuel cell systems by multi-objective optimization. J. Power Sources 166 (2), 411–418. 10.1016/j.jpowsour.2007.01.022 Öztürk A. Fıçıcılar B. Eroğlu İ. Bayrakçeken Yurtcan A. (2017). Facilitation of water management in low Pt loaded PEM fuel cell by creating hydrophobic microporous layer with PTFE, FEP and PDMS polymers: Effect of polymer and carbon amounts. Int. J. Hydrogen Energy 42 (33), 21226–21249. 10.1016/j.ijhydene.2017.06.202 Parekh A. (2022). Recent developments of proton exchange membranes for PEMFC: A review. Front. Energy Res. 10. 10.3389/fenrg.2022.956132 Sauermoser M. Kizilova N. Pollet B. G. Kjelstrup S. (2020). Flow field patterns for proton exchange membrane fuel cells. Front. Energy Res. 8. 10.3389/fenrg.2020.00013 Wang Z. Li L. (2021). Optimization of process parameters for surface roughness and tool wear in milling TC17 alloy using Taguchi with grey relational analysis. Adv. Mech. Eng. 13 (2), 168781402199653–168781402199658. 10.1177/1687814021996530 Xu M. Zhou C. Huang X. Zhang Z. Wang T. (2020). Multiobjective optimization of 316L laser cladding powder using gray relational analysis. J. Mater. Eng. Perform. 29 (12), 7793–7806. 10.1007/s11665-020-05278-8 Xu Y. Chang G. Zhang J. Li Y. Xu S. (2021). Investigation of inlet gas relative humidity on performance characteristics of PEMFC operating at elevated temperature. World Electr. Veh. J. 12 (3), 110–121. 10.3390/wevj12030110 Yang Y. Chen Z. Sun X. Yao S. Zhang X. Liu W. (2023). Li4SiO4 adsorbent derived from industrial biomass fly ash for high-temperature CO2 capture. Fuel 337, 126853–126862. 10.1016/j.fuel.2022.126853 Yao J. Wu Z. Wang H. Yang F. Xuan J. Xing L. (2022). Design and multi-objective optimization of low-temperature proton exchange membrane fuel cells with efficient water recovery and high electrochemical performance. Appl. Energy 324, 119667–119685. 10.1016/j.apenergy.2022.119667 Yin Y.-m. Cui H.-y. Hong M. Zhao D.-y. (2015). Prediction of the vertical vibration of ship hull based on grey relational analysis and SVM method. J. Mar. Sci. Technol. 20 (3), 467–474. 10.1007/s00773-014-0299-5 Yong Z. Shirong H. Xiaohui J. Yuntao Y. Mu X. Xi Y. (2022). Performance study on a large-scale proton exchange membrane fuel cell with cooling. Int. J. Hydrogen Energy 47 (18), 10381–10394. 10.1016/j.ijhydene.2022.01.122 Zeng X. Ge Y. Shen J. Zeng L. Liu Z. Liu W. (2017). The optimization of channels for a proton exchange membrane fuel cell applying genetic algorithm. Int. J. Heat Mass Transf. 105, 81–89. 10.1016/j.ijheatmasstransfer.2016.09.068 Zhang Q. Lin R. Técher L. Cui X. (2016). Experimental study of variable operating parameters effects on overall PEMFC performance and spatial performance distribution. Energy 115, 550–560. 10.1016/j.energy.2016.08.086 Zhang T. Wang P. Chen H. Pei P. (2018). A review of automotive proton exchange membrane fuel cell degradation under start-stop operating condition. Appl. Energy 223, 249–262. 10.1016/j.apenergy.2018.04.049