Data-driven sensitivity analysis of contact resistance to assembly errors for proton-exchange membrane fuel cells

The proton-exchange membrane fuel cell is a promising power source for automobile industry because of its zero pollution. However, its stack structure always faces increased contact resistance caused by assembly errors, leading to substantial energy loss during the working period. To enhance its output performance, the influence of assembly errors on contact resistance is studied for proton-exchange membrane fuel cell. The mechanical simulation model of fuel cell assembly process is established to provide contact resistance distribution with different assembly errors. An improved global sensitivity analysis method is proposed to evaluate the influence coefficient of each assembly error term on contact resistance based on a series of randomized simulation data. The case study of a single-layer fuel cell demonstrates the proposed method achieves higher efficiency than traditional sensitivity analysis methods, and finds out key assembly errors in regard to reducing contact resistance.


Introduction
The fuel cell is a promising power source for mobile transportation because of its advantages such as zero pollution, high efficiency, and portability. [1][2][3] In particular, proton-exchange membrane fuel cell (PEMFC) has been attracting extensive attentions in recent years because of its low temperature operation, high power and rapid start-up. [4][5][6] However, due to its own resistance, PEMFC has to bear serious energy loss during its working period. For instance, the contact resistance between metal bipolar plate (BPP) and gas diffusion layer (GDL) always causes about 60% energy loss of a PEMFC. 7 This research is intended to solve the analysis of contact resistance for PEMFC.
Based on theoretical analysis, assembly errors are the major factors leading to contact resistance because they make PEMFC's stack structure cannot have adequate contact area between BPP and GDL. 8 Assembly error control is an effective and practical manner to reduce contact resistance, but will spend unacceptable production cost when all the assembly errors are taken into consideration. 9,10 Since assembly errors regularly have different impacts on contact resistance, it is more efficient to implement control measures on assembly error terms having comparatively higher impacts on contact resistance. Therefore, the study of influence degrees of assembly errors on contact resistance is of great significance for improving PEMFC performance.
A great deal of research efforts has already been devoted to reveal the formation of contact resistance, especially between BPP and GDL in PEMFC. For example, Zhou et al. 11 used finite element method to establish a PEMFC contact model and studied contact resistance with different clamping forces. Lai et al. 12 established an electromechanical finite element model to predict the contact resistance between BPP and GDL. Liang et al. 13 established a three-dimensional finite element model to study the influence of manufacturing characteristics of BPP on contact resistance. They predicted the contact resistance caused by different dimensional errors of coating, weld, and metal BPP. Zhou et al. 14 designed a microscale numerical model to evaluate the contact resistance between BPP and GDL. Qiu et al. 15 investigated the contact condition between solid and porous material based on microscopic theory, and established a microscale contact model to predict the contact resistance between BPP and GDL. Li et al. 16 used finite element method to study the contact resistance between BPP and GDL under different bolt clamping forces. Wu et al. 17 studied the formation of contact resistance between BPP and GDL, and realized regression analysis of contact resistance under different assembly pressures. In general, these efforts obtain contact resistance mainly through simulation model and investigate the distribution of contact resistance when single influential parameter changes. However, these tests collect only small-sized data because simulation model construction is timeconsuming, and the analysis methods cannot give indepth analysis of contact resistance with simultaneous variation of multiple influential parameters.
Since assembly errors are controlled within a feasible interval, rather than given a fixed value, the relationship between assembly errors and contact resistance requires a quantitative analysis that can take all the assembly error terms into consideration. Of existing quantitative analysis methods including regression analysis, decision tree, Pareto diagram, and factorial design, global sensitivity analysis (GSA) method can evaluate the influences on contact resistance with the randomization of multiple assembly errors. [18][19][20][21] However, this method regularly requires a large number of data samples to achieve accurate estimation of influence degrees, which makes it impractical for some complex relationships. 22 To overcome this shortcoming, this paper attempts to collect more data samples through secondary development of the simulation model, and proposes an improved global sensitivity analysis (IGSA) method to evaluate the influence between assembly errors and contact resistance more efficiently.
The rest of the article is organized as follows. Section ''Problem description and research framework'' gives in detail the problem description and research framework. Section ''PEMFC simulation model'' presents the PEMFC simulation model. Section ''IGSA method between assembly errors and contact resistance'' proposes an IGSA method for this problem. Section ''Case study'' provides the case study. Conclusions and future research directions are discussed in section ''Conclusion.''

Problem description and research framework
As shown in Figure 1, the stack structure of a singlelayer PEMFC mainly consists of BPP that provides the flow channels, a membrane exchange assembly (MEA), some sealing originals and two end plates. [23][24][25] The MEA is composed of a proton-exchange membrane, catalyst layers and GDL.
In most practical applications, multiple layers of PEMFC are connected together to provide enough power and voltage. 26 However, the PEMFC cannot afford enough output power when the energy loss caused by contact resistance is too large in some cases. Since assembly errors prevent sufficient contact between BPPs and MEA and are the major cause of contact resistance, it is important to find out effective assembly error control measures for contact resistance reduction. For this purpose, this research deals with the influence analysis problem between assembly errors and contact resistance for single-layer PEMFC. Table 1 lists some symbols used in this problem.
For a single-layer PEMFC, the assembly errors regularly include dimensional errors of BPP, sealing originals and GDL, as well as assembly deviation between BPPs. Dimensional errors of BPP and assembly deviation between BPPs are major causes of contact resistance and therefore considered in this research. 27,28 Dimensional errors of BPP Dimensional errors of BPP refer to the inconsistence of channel heights, which is caused by the elastic recovery of metal material and uneven pressing force during the stamping process. To describe dimensional errors, two reference lines are defined. The lower datum line describes the bottom of BPP, and the upper datum line describes the standard height of flow channels. The difference between a flow channel's height and the standard height is defined as the dimensional error Dt where H r BPP represents the height of a flow channel and H n BPP represents the standard height of all flow channels, as shown in Figure 2.

Assembly deviations between BPPs
Due to inaccurate positioning size, the disturbance from external vibration, and the compressive deformation of stack structure, horizontal displacement happens between the flow channels of two adjacent BPPs, which causes assembly deviation of single-layer PEMFC. To describe assembly deviation, the central axis of each flow channel of BPP gives the reference line, as shown in Figure 3. The distance d between reference lines of upper and lower BPPs is defined as assembly deviation.
Taking above assembly errors into consideration, this research includes two parts: 1. Based on qualitative analysis of PEMFC structure, the parametric mechanical simulation  Assembly deviations between BPPs in flow channels 1,2, ..., K P 1 , P 2 , . . . , P M Contact stress of nodes 1,2, ..., M in finite element model Contact resistance of nodes 1,2, ..., M in finite element model X 1 , X 2 , . . . , X D Input parameters in the influence analysis problem X i1 , X i2 , . . . , X iS Sub-set of input parameters, fi 1 , i 2 , :::, i S g 2 f1, 2, :::, Dg u Contact resistance of PEMFC V Total variance of contact resistance caused by X 1 , X 2 , :::, Local variance of contact resistance caused by X i1 , X i2 , . . . , X iS e 1 , e 2 , . . . , e N Randomized sample points in Monte Carlo estimation r 1 , r 2 , . . . , r D Impact degrees of assembly errors X 1 , X 2 , . . . , X D on contact resistance u MinE Constant number, the minimum sample points in each subdomain   model between BPP and GDL is constructed in finite element analysis software ABAQUS. The variation of assembly errors is enabled by secondary development of this simulation model using Python language. According to assembly errors, contact stresses are simulated at GDL, and contact resistance is calculated using an empirical formula. 2. Based on mathematical description of assembly errors, the calculation formula of influence degrees of assembly errors on contact resistance are defined. Considering the complex relationship between assembly errors and contact resistance, Monte Carlo algorithm (MCA) is used to calculate the influence degrees through contact resistance variance estimation, of which sample data are collected from the mechanical simulation model. In the collection of sample data, an adaptive Sobol sequence is specially designed for assembly error randomization to improve the efficiency of MCA.

PEMFC simulation model
This section establishes the quarter-stack simulation model of a single-layer PEMFC to provide contact resistance with different assembly errors.

Finite element model
The contact model between BPP and GDL is constructed in the finite element analysis software. In this model, the end plates are considered as rigid bodies because their thickness and stiffness are very large when compared with BPPs. 29,30 A rigid line is built on the top of BPP to simulate the implement of assembly pressures. As shown in Figure 4, quadrilateral unit CPE4R of two-dimensional plane strain is used in the gridding of BPP and GDL. The stress concentration parts use subdivided grids to ensure simulation accuracy, while other parts use larger grids to save computational time.
To restrict the movement of this model, zero displacement is applied to the left and right boundaries in the x direction, and fixed constraint is applied to the bottom of GDL in the y direction. The displacement load of 0.1 mm was applied to the rigid line to simulate assembly pressure.
Moreover, the secondary development in Python library package included in ABAQUS realizes the read of multiple sets of BPP and GDL parameters in Excel files, the generation of finite element model with specific geometric parameters and material property, and the simulation of contact stresses between BPP and GDL through finite element analysis.

Contact resistance calculation
For contact nodes between BPP and GDL, the contact resistance of each node is calculated based on following regression equation 31 where P m represents the mth node's contact stress obtained by finite element analysis, and a and b are fitting parameters, as shown in Figure 5.
By regarding the contact interface between BPP and GDL as an equipotential surface, the parallel resistance of these contact nodes gives the contact resistance of PEMFC 32 where R m represents the mth node's contact resistance, and M represents the number of contact nodes.

Assembly error description
The IGSA method describes two kinds of assembly errors, including dimensional errors of BPP and assembly deviations between BPPs.
Dimensional errors of BPP. Based on equation (1), dimensional errors of BPP are represented by Dt 1 , Dt 2 , . . . , Dt K , where K represents the number of flow channels in BPP.
where D = 23K represents the number of input parameters, U D represents the normalized space of input parameters, u represents contact resistance, X 1 , . . . , X D represent dimensional errors and assembly deviations, f 0 is a constant, and f i (X i ) represents the local variance of contact consistence owing to assembly error X i . Local variance of contact resistance resulted by some assembly errors and total variance of contact resistance resulted by all assembly errors are calculated by equations (5) and (6), respectively 34 Based on contact resistance variance, the first-order sensitivity coefficient and the total sensitivity coefficient are calculated in equations (7) and (8), respectively Monte Carlo estimation. Due to nonlinear and complicated relationship between assembly errors and contact resistance, the mathematical model u = f(X) can hardly be expressed by analytic equations. For this reason, regression analysis of accurate formula of this model is almost infeasible in the calculation of contact resistance variances. On the contrary, MCA can evaluate the distribution of function values based on a series of randomized sample points located in the feasible domain of input parameters, without depending on analytic equations. 35 And, the simulation model can provide the oneto-one correspondence of contact consistence based on randomized assembly errors. Thereupon, MCA is used in the estimation of integrals in equations (5) and (6) 1 N 1 N where U D represents the integral domain, N represents the number of random points, e n = (e 1 n , e 2 n , . . . , e d n , . . . , e D n ) represents a random point in the integral domain, e d n represents a randomized value of the dth assembly error X d , u n = (u i 1 n , u i 2 n , . . . , u i s n ) represents randomized values of assembly errors X i 1 , . . . , X i s and h n represents randomized values of remaining assembly errors of X = (X 1 , ..., X D ).
Random points distribution. In MCA, the contact resistance of each random point is calculated using PEMFC simulation model. Massive random points can realize accurate estimation of contact resistance variance, but will spend impractical simulation time. To save simulation time, the Sobol sequence generates uniformly distributed pseudo-random points in the integration domain and achieves more accurate estimation of integrals than other quasi-random methods with the same number of random points. 36 In PEMFC simulation model, specific assembly error distribution can realize the minimum contact resistance, which corresponds to the valley point of u = f(X) in the integration domain. 37 For integration domains near valley points, that is, subdomain C in Figure 6, the Sobol sequence needs to generate more random points to reduce estimation errors.
However, this will cause the increase in random points in other subdomains because the Sobol sequence needs to keep uniform distribution of quasi-random points in the entire integral domain. These increased sample points will lead to longer simulation time, but realize very limited improvement of estimation accuracy in subdomains without valley points, that is, subsection A. To improve analysis efficiency, the IGSA method generates random points adaptively in Monte Carlo estimation. In IGSA procedure shown in Figure 7, the MCA finds out valley regions based on the variance of norm vectors between each pair of random points in subdomain, and generates new randomized sample points only in these regions to reduce estimation errors efficiently.
After IGSA procedure, the influence degree matrix ½r 1 , . . . , r d , . . . , r D of assembly errors (X 1 , . . . , X d , . . . , X D ) on contact resistance u is obtained for single-layer PEMFC, where r d represents influence degree of assembly error X d on contact resistance u.

PEMFC parameters
Based on practical data collected from an automobile enterprise, the contact model between a standard nominal-size metal BPP of stainless steel of SS316L and GDL of TGP-H-060 is constructed in ABAQUS, as shown in Figure 8.

Experiments and discussion
Comparative experiments. Table 3 illustrates estimation errors of contact resistance variance using traditional GSA method and IGSA method, separately. The initial number of sample points changes from 600 to 1200 in IGSA method, and GSA method generates the same number of sample points with IGSA method to compare the estimation accuracy. With the increase in the sample points, the estimation of either method decreases, and IGSA method realizes smaller estimation errors than GSA method. As shown in Figure 9, when more than 1100 sample points are distributed in the integration domain, IGSA method realizes the estimation error less than 0.001, which ensures accurate evaluation of influence degrees.
Analysis results and discussion. Based on comparative experiments, IGSA method gives influence degrees of assembly errors on contact resistance by generating 1100 sample points, as shown in Table 4.
The results demonstrate that the seventh, ninth, 12th, and 16th dimensional error parameters and the third, ninth, and 19th assembly deviation parameters have comparatively greater influences on contact resistance. To improve PEMFC performance, the height deviation of the seventh, ninth, 12th, and 16th flow channels of BPP should be strictly controlled in the stamping and forming process. Besides, the positioning accuracy of ridge centers of the third, ninth, and 19th flow channels should be ensured in the assembly of BPPs.
In Table 5, the control of assembly errors realizes contact resistance reduction effectively. According to contact resistances of different trials, dimensional error control of Dt 7 , Dt 9 , and Dt 12 is more effective than assembly deviation control, which validates the influence degrees in Table 5. This is because the decrease in the negative dimensional error parameters realizes more adequate contact between BPP and GDL, and thus reduces contact resistance. On the contrast, Table 6    provides the contact resistance with Dx j j= 0:02mm decrease of two assembly error parameters. Table 6 demonstrates that the control of two assembly error parameters is more effective than that of only one assembly error parameter. Moreover, the combination of assembly error controls in the 12th-dimensional error parameter and the ninth assembly deviation parameter realizes minimum contact resistance. The reason might be the assembly error control in the middle part of BPP and GDL realizes more complete contact not only in the ninth to 12th flow channels but also in other flow channels. Moreover, the averaged contact resistance of No. 2 This comparison result demonstrates that dimensional error control is more effective than assembly deviation control, which validates the difference in impact degrees between Dt 12 and d 9 in Table 4.

Conclusion
The influence degrees of assembly errors on contact resistance is analyzed for PEMFC, which informs assembly error control measure that can reduce contact resistance effectively and realize enough output power under given voltage. The key contribution is an IGSA   Table 6. Contact resistance with two decreased assembly error parameters.
No. method that designs an adaptive Sobol sequence to collect randomized data more effectively from PEMFC simulation model and evaluates the influence degrees via Monte Carlo estimation of contact resistance variance based on randomized data. The case study validates the effectiveness of this IGSA method over traditional method, and demonstrates the application of influence degrees of assembly errors on contact resistance in certain scenarios. Our further research includes influence analysis for PEMFCs with 5-10 layers. In addition, we are also trying the substitution of simulation model with some machine-learning methods, for example, deep neural network, since the simulation in ABAQUS becomes much more time-consuming for multiple-layer PEMFC.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.