Alexandra Pessot
email: alexandra.pessot@irt-saintexupery.com

Christophe Turpin

Amine Jaafar

Emilie Soyez

Olivier Rallières

Guillaume Gager

Julien D'arbigny

Julien D'arbigny

Contribution to the modelling of a low temperature PEM fuel cell in aeronautical conditions by design of experiments

Keywords: PEM fuel cell, aeronautical application, Design of Experiments, modelling

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

density, using response surface methodology. The authors in reference [START_REF] Giner-Sanz | Statistical Analysis of the Effect of the Temperature and Inlet Humidities on the Parameters of a PEMFC Model[END_REF] performed on a stack a full factorial plan of 2 levels for 3 factors: temperature, inlet hydrogen humidity and inlet air humidity. They established a semi-empirical model with 6 parameters using in particular an Amphlett expression for activation overpotential. They fitted this model to their experimental database and they carried out a statistical analysis in order to study the influence of each factor on cell performance. Then, they used the response surface method to link model parameters to applied operating factors, it is to say to see which operating conditions have an impact on a considered model parameter. Finally, they obtained a model which is able to estimate the experimental polarization curves. One of the limitations of the proposed model is that it is only valid for low current densities. In the thesis work of Isabelle Labach carried out at the LAPLACE laboratory, a model of H 2 /Air PEM fuel cells in variable operating conditions was proposed [8] [9]. This model was validated with single cells. Considered impacting factors were temperature, pressure, air relative humidity and air stoichiometry. The proposed model was parameterized thanks to a database created with Design of Experiments method. Experimental relations were then established to express the variation of the model parameters as a function of the operating conditions. This model allowed estimating with satisfactory results the polarization curves of a PEM single cell which operates in a given range of operating conditions. In addition, predictive behaviour of this model was investigated and results seemed promising to estimate polarization curves of cells operating outside the defined range.

In our work, a particular attention is paid to the study of low-temperature PEM fuel cells in aeronautical conditions. Indeed, nowadays, fuel cells are getting an increasing amount of interest from the actors of the aeronautical industry and their integration in airplanes is considered as a consistent solution [START_REF] Sehra | Propulsion and power for 21st century aviation[END_REF] [START_REF] Pratt | Proton exchange membrane fuel cells for electrical power generation on-board commercial airplanes[END_REF]. However, even if this technology is quite mature in automotive and stationary sectors, there are still many questions about how the operating conditions specific to aeronautical environment will influence the functioning and performance of fuel cells [START_REF] Renouard-Vallet | Improving the environmental impact of civil aircraft by fuel cell technology: concepts and technological progress[END_REF]. Pressure is an operating factor that requires a remarkable thinking in aircraft applications. On one hand, increasing pressure leads to an increase of fuel cell performance. Several works have pointed out an important decrease of performance when pressure decreases [START_REF] Horde | PEMFC application for aviation: Experimental and numerical study of sensitivity to altitude[END_REF] [START_REF] Pratt | Performance of proton exchange membrane fuel cell at high-altitude conditions[END_REF], an observation that has been largely reported by the DLR teams [START_REF] Kallo | Fuel cell System Development and Testing for Aircraft Applications[END_REF] [START_REF] Werner | A comparison of low-pressure and supercharged operation of polymer electrolyte membrane fuel cell systems for aircraft applications[END_REF]. On the other hand, operating with pressurized gases requires the use of a compressor, what means supplementary weight and volume. At an altitude of around 10 000 km, ambient pressure is close to 0.2 bar, whereas pressure is higher than 1.3 bar in a classical operating range for PEM fuel cells. A compromise could be to extract air from the aircraft cabin where ambient pressure is around 0.8 bar.

In our present project, a benchmark study has been performed in order to compare performance of several types of Membrane Electrode Assemblies (MEAs) within stacks having rigorously the same structure (bipolar plates, cooling circuit, etc.) under aeronautical conditions for cabin applications. These stacks have been supplied by ZODIAC AEROTECHNICS which develops PEM fuel cell systems dedicated to aeronautical applications. An experimental database has been constituted thanks to the Design of Experiments (DoE) method. This database is composed of several polarization curves obtained in different operating conditions for different stacks, with in particular several points at subatmospheric pressures. All these tests were carried out at the LAPLACE laboratory. The same test bench was used for all these tests.

This paper focuses on a part of the database; in other terms, this article will deal with the campaign achieved with one stack with one type of MEAs. Our objective in the present paper is indeed to describe the methodology we are following to analyze the created database, not to compare all the tested types of MEAs. Ideally, we want to define a global, robust and physics-linked model for the polarization curve able to describe its evolution when the operating conditions change. In other words, it is about ideally to obtain a physical law for each parameter of the model according to the operating conditions (T,P, gas flows, gas humidities). Having such a model permits afterwards to optimize, mainly from the efficiency point of view, the system around the fuel cell stack as the air compressor, the humidification system, the cooling system, etc. Unfortunately, this very ambitious objective has not been yet reached, but the results which are presented in this paper are a new step towards this goal.

Firstly in this paper, the carried out DoE is presented, followed by the formulation of the modelling we want to obtain. These two points being presented, the difficulty we have to overcome is that the physical laws describing the dependence of the model parameters to the operating conditions are not really formulated in a general point of view. This is linked to the complexity of the fuel cell operation and the strong couplings in terms of impact between all the operating conditions. That is the reason why the approach proposed in this paper is very progressive, in order to stay general and applicable to other characterization campaigns. We have adopted a step-by-step approach by introducing progressively hypotheses to analyze the experimental results. Thus, we start by analyzing one by one each measured polarization curve to identify trends concerning the dependence of the model parameters to the operating conditions. From these trends, we simplify for each model parameter the number of links with the operating conditions. Afterwards, we re-analyze all the polarization curves simultaneously with the model simplified (even if it remains complex). These last results being obtained, we start to identify laws describing the dependence of the model parameters to the operating conditions. Each time, an exploitation typical of the DoE methodology is proposed, that means the determination of an empirical law assuming a linear dependence to the considered operating conditions. If possible, another exploitation is proposed in parallel by introducing a physical law. We will see that this was only possible for a single parameter related to chemical kinetics. At the end of the paper, the obtained model is tested concerning its predictive ability.

Presentation of the experimental approach and of the model used

Experimental approach

Tests were organized following the Design of Experiments approach on the same test bench at the LAPLACE laboratory with a 1kW PEMFC stack. The impact on fuel cell performance was studied for four operating conditions: fuel cell temperature T, gas pressure P, air relative humidity RH air and hydrogen relative humidity RH H2 . Tests were conducted at constant stoichiometries whatever the conditions: λ H2 = 1.5 for hydrogen and λ air = 2 for air.

Except for pressure, two levels were chosen. Four levels were selected for pressure which impact is particularly not well known amongst all the aeronautical conditions. Table 1 shows the selected values for each factor. The lowest pressure (0.8 bara) is typical of the cabin pressure. The highest pressure is in the range of classical values used in terrestrial applications. The values of temperature are in a classical range for PEMFC operating at low temperatures. Gas relative humidity can seem rather low compared with classical values, but the idea was here to investigate performance at low relative humidity. In practice, it was rather difficult, taking into account our humidification systems, to apply this range of relative humidity. We estimate that, depending on temperature and pressure, we were between 13% and 18% when 16% was the objective and between 22% and 30% when 26% was the target. In other words, the precision was from ±10 to ±20% for our gas relative humidity. In the following, the supposed reached experimentally value of HR is given with the initial target "16%" or "26%".

In the DoE approach, the normalized operating conditions are classically considered. For an operating factor X i (X i = T, P, RH air or RH H2 ), the corresponding normalized value x i (x i = t, p, rh air or rh H2 ) is given by:

x X X X 2 X X 2 ( 1 )
Where x i is a value comprised between 0 and 1, X max and X min are respectively the maximal and the minimal value of X i .

Concerning the DoE experiments achieved, only 31 polarization curves have been measured whereas 32 tests should have been carried out theoretically (multiplication of all the levels for all the considered factors: 2x4x2x2=32). Indeed, curve corresponding to the operating point (t, p, rh air , rh H2 ) = (1, -1, -1, -1) could not be tested because of too low performances (cell voltage security level has been reached). In other words, the operating conditions were too constraining: too high temperature, too low pressure and too low relative humidities.

Let us underline that Electrochemical Impedance Spectroscopies (EIS) were also carried out at the same time that the measurement of a polarization curve for each current, step by step. These EIS will be exploited here to obtain the global ohmic resistance (named Rohm in the following) including ionic and electronic phenomena.

Only the significant and representative results of the carried out DoE are presented in the following. Figure 1 presents all the curves obtained for the highest gas relative humidities supposed to be around 26%. At given temperature, the ranking of measured performance according to pressure is very clear: the more the pressure increases, the better the performance. Unfortunately, performance decreases very strongly when pressure decreases under subatmospheric pressures. These results are very consistent with the state-of-art [START_REF] Horde | PEMFC application for aviation: Experimental and numerical study of sensitivity to altitude[END_REF] [15] [START_REF] Yan | Steady state and dynamic performance of proton exchange membrane fuel cells (PEMFCs) under various operating conditions and load changes[END_REF]: the O 2 concentration is too low at low pressures to obtain satisfying operation.

At given pressure, the ranking according to temperature is not really the expected one: each time, the best performance is obtained for 65°C the lowest temperature tested. Generally, an increase of temperature has indeed a benefic effect on the performance of a fuel cell because a lot of phenomena are favored: kinetics of the involved chemical reactions, protonic conduction and the involved species transport because water in liquid form decreases. That means that the operation of the tested MEAs is more complex than that of the classical ones. At low current densities where dynamics is dominated by the kinetics of the involved chemical reactions (activation phenomena), it appears that the operation is clearly sensitive to pressure, and in a more complex way to temperature. These trends were not expected illustrating all the operation complexity of a fuel cell.

For some selected points of the carried out DoE, Figure 2 illustrates the influence of gas humidities.

Firstly, let us underline that, in this figure, are drawn the best and the worst polarization curves obtained during this campaign. The best one corresponds to 65°C, 1.5 bara and gas humidities of 26%: except temperature, that was obtained for the highest values for the operating conditions. The worst one was obtained for 80°C, 0.8bara, and gas humidities of 16%/26% (air/H 2 ): except temperature, this was for the lowest values tested in terms of operating conditions (let us recall that the case 80°C, 0.8bara and RH of 16%/16% (air/H 2 ), the worst case possible in our DoE, was not possible to measure).

Figure 2 -Selected experimental polarization curves of the DoE: illustration of the influence of gas humidities

Let us now observe the influence of gas humidities. In Figure 2, three series of curves can be clearly separated corresponding respectively to the three couples of operating conditions (65°C, 1.5bara), (65°C, 1bara) and (80°C, 0.8bara). For each couple, the ranking for the humidity influence is always the same: the best performance is when air relative humidity has its highest value (26%). The highest value for hydrogen relative humidity (26%) improves also performance, but the influence air relative humidity is priority. The impact of gas humidities is not chaotic, but ordered or even logical.

Figure 3 proposes the evolution of R ohm according to the current for selected operating conditions. Let us recall that R ohm is measured by EIS. When no value was measured at low current densities, the value considered for R ohm is that of the lowest current density measured; this value is arbitrarily blocked to that value if necessary.

The evolution according to the current is globally classical: at low current densities, less water is produced and moreover it is difficult to maintain constant the gas overstoichiometries, leading to a progressive drying of the membranes and thus an increase of R ohm when the current tends towards zero. On the other hand, the behavior of R ohm is a little bit amazing concerning temperature: indeed temperature is usually favorable to decrease Rohm, contrarily to the trends observed here. Pressure has a nonlinear and very important impact on the value of R ohm as it can be seen clearly in Figure 3.a: the more the pressure increases, the more R ohm increases. The influence of gas humidities (Figure 3.b) is again here coherent and logical: the best results are obtained when gas humidities have their highest value (26%) and the worst ones for their lowest value (16%).

The complex modelling of Rohm is not an objective in this paper. In the following, these measures of Rohm will be directly injected into our modelling works.

Presentation of the model used for the polarization curve

The global model equation is:

( 2 )

Where: -U cell : fuel cell voltage (V) -E rev : theorical thermodynamic reversible voltage (V) -η act , η diff , η ohm : activation, diffusion and ohmic losses (V):

The thermodynamic reversible voltage E rev can be defined thanks to the following relation:

° ! " #$ % & '$ ( & ) * +, ( 3 )
The activation losses, represented by η act , are linked to the chemical reactions' kinetics. Thanks to a simplified Butler-Volmer equation, η act can be written: Where the parameters α, J n and J 0 are defined in Table 2.

The diffusion losses η diff refer to the reactant gas transport in the Gas Diffusion Layers (GDLs) and in the Active Layers (ALs). These losses are mainly occuring at the cathode side. Combining Fick and Faraday laws, and assuming that diffusion losses are mainly occuring at the cathode side, the following relation can be proposed for η diff :

(.) 3 4 ! " '1 . . +3 ( 5 )
Where the parameters β and J lim are defined in Table 2.

The ohmic losses η ohm include voltage drops linked to the membrane resistance, to the electronic resistance of the other elements of the cell (electrodes, bipolar plates, etc.) and to contact resistances existing between the different layers. The ohmic losses η ohm can be expressed by: (.) .

(

) 6
Where the parameter R ohm is defined in Table 2.

Table 2 -Signification of terms presented in the previous equations

The fuel cell voltage can be thus expressed with the following relation:

(.) 0 ! " ' . . 1 . 2 + 3 4 ! " '1 . . +3 . ( 7 )
The parameters to be identified are a priori α, J 0 , J n , β, J lim and R ohm .

In accordance with many bibliographic studies [18] [19], in this paper, the charge transfer coefficient α was assumed constant with current, invariant with the operating conditions and equal to 0.5.

Moreover, as already explained in part 2.1, R ohm is not identified but measured thanks to Electrochemical Impedance Spectroscopy (EIS); it corresponds to the intersection of the spectrum at high frequencies with the abscissa axis. R ohm is measured at every current density point of the polarization curve. Finally, the model parameters to be identified are a priori J 0 , J n , β and J lim . The objective in the following of this paper is to determine their dependence to the operating conditions from the analysis of the carried out DoE.

Modelling of the results for the DoE experimental database

General comments

The diffusion losses are nonlinear from a general point view as in part 2.2 (equation ( 5)). In the database (see Figure 1 and Figure 2), as there is no inflection/bending of the polarization curves at high current density, β and J lim are impossible to determine because an infinity of solution exists. That is why a linear relation is used for diffusion losses instead of the equation ( 5):

(.) . ( 8
)
Where the diffusion resistance R diff (Ω.cm²) is expressed by 67 189: ;<= . So the following relation will be finally used to model a polarization curve for our DoE:

(.) 0 ! " ' . . 1 . 2 + . . ( 9 )
In this case, parameters which have to be identified are J 0 , J n , and R diff .

At this step of our study, a general relation taking into account the parameter dependence to the operating conditions can be formulated:

(., , ?, @ , @ %* ) ( , ?) 0 ! " A . . 1 ( , ?, @ , @ %* ) . 2 ( , ?, @ , @ %* ) B ( , ?, @ , @ %* ) . ( , ?, @ , @ %* ) .

In order to simplify these dependences to all the operating conditions for J 0 , J n , and R diff , a first modelling is going to be performed firstly curve by curve. It consists in identifying, with an optimization algorithm, the three values J 0 , J n , and R diff curve by curve for the 31 measured curves of the DoE. The approach is afterwards to observe if the impact of an operating condition on a parameter can be neglected.

The optimization algorithm used in our modelling approach is the Covariance Matrix Adaptation Evolution Strategy, called "cmaes". An objective function gives the error to be minimized during the optimization process. Several starting points are randomly drawn in the fixed variation domain (lower and upper bounds are defined for every parameter) and at the end, only the best solution is considered. In order to converge more quickly, parameters used in the optimization process are expressed in the logarithmic scale and are normalized. In the following sections, parameters are always presented with their normalized values for confidentiality reasons.

Curve-by-curve identification

In this section, a curve-by-curve identification is carried out: a parameter set (J 0 , J n , R diff ) has to be identified for each polarization curve, meaning that 31 J 0 , 31 J n and 31 R diff have to be identified.

The error ε which is calculated in the objective function and which has to be minimized in the optimization process is based on a least square criterion:

C D D A _ F (.) _ GH (.) _ GH (.) B * : 1 IJKLMN ( 11 )
Where n curves is the number of curves, U cell_exp is the experimental cell voltage (V) and U cell_est is the estimated cell voltage calculated using the parameters obtained in the identification process (V). Estimated cell voltage is the calculated cell voltage using the model equation ( 9) presented in section 3.1.

In the following work, identified parameters will be presented in their normalized value x norm compared to their lower and upper bounds fixed in the identification process:

O 1 P O Q O RQ O Q ( 12 )
With X the parameter (X = J 0 , J n or R diff ), x lb the lower bound of the parameter, x ub the upper bound of the parameter and x norm the normalized value of the parameter.

Results of the curve-by-curve identification are presented in Table 3. Even if not directly visible in Table 3, we can affirm that the orders of magnitude of identified parameter values are coherent with the typical values found in the literature for each parameter: identified J 0 are around 5.10 -7 A/cm², J n are around 1.10 -4 A/cm² and R diff are around 0.2 Ω.cm². Maximal relative error between experimental cell voltage and estimated cell voltage of the curve 1 is 0.42%. Maximal relative error of the curve 31 is 0.92%. Similar results are obtained for every polarization curve of the DoE. Mean relative error between experimental and estimated polarization curves calculated for the 31 DoE curves is 0.24 %. In all, maximal relative error at a given current density is 1.26 % (obtained for curve 30).

Curve

As a conclusion concerning this first identification, we can say that the model, with parameters identified curve by curve, enables to estimate with a satisfactory precision the polarization curve of the studied fuel cell at a given operating point. However, the simulation time is important: it takes around 3 hours to make the identification for all the 31 curves (93 parameters have to be identified in all in this case).

Let us try now to put in evidence simplifications in terms of model parameters' dependence to operating conditions.

The results in table 3 show that for experiments performed at given temperature and pressure (red rectangles), the identified J 0 values are really similar regardless the value of gas humidities. That shows that it is possible to make the assumption that J 0 depends only on temperature T and on pressure P. In the same way, the assumption considering that J n depends on temperature T, on pressure P and on air relative humidity RH air can be proposed. Looking at the identified values of R diff , it appears that it is quite difficult to detect a trend concerning the dependency of this parameter to particular operating conditions. That is why R diff is considered to be dependent on all operating conditions (T, P, RH air and RH H2 ). Theses hypotheses are consistent with first analyzes of the experimental polarization curves presented in section 3.1.

Consequently, the previous equation [START_REF] Sehra | Propulsion and power for 21st century aviation[END_REF] proposed in section 3.2 can be modified in order to taking into account these new hypotheses, globally consistent with the literature [START_REF] Wahdame | Study of gas pressure and flow rate influences on a 500W PEM fuel cell, thanks to the experimental design methodology[END_REF]

(., , ?, @ , @ %* ) ( , ?) 67 S18 " T :U: V (7,W,6% X<K ) : Y (7,W) Z ( , ?, @ , @ %* ) .

( , ?, @ , @ %* ).

( 13 )

Identification by applying the new hypotheses

In this section, new model identification is going to be carried out by applying the new hypotheses that emerged at the end of section 3.2. This time, the identification is not carried out curve by curve since some parameters are supposed to be common for several curves. A total of 8 J 0 , 16 J n and 31 R diff have to be identified.

Table 4 -Results of identification with the new hypotheses for every curve of the DoE

The results of this identification are presented in Table 4. Table 4 has to be compared with Table 3. For example, the new values of J 0 are very close to the average value of the previouslygrouped values in Table 3. This observation is the same for J n . The values R diff are very close to the ones of Table 3.

Figures 6 and7 show the comparison between experimental polarization curve and estimated polarization curve obtained with the identification for curves 1 and 31 of the DoE.

Figure 6 -Experimental and estimated polarization curves (identification applying new assumptions on parameter dependences on operating conditions)

Curve 1: T = 80°C, P = 1.5 bara, RH air = 30 % (/"26%"), RH H2 = 30% (/"26%")

Figure 7 -Experimental and estimated polarization curves (identification applying assumptions on parameter dependences on operating conditions)

Curve 31: T = 65°C, P = 0.8 bara, RH air = 13 % (/"16%"), RH H2 = 13 % (/"16%")

Relative error between experimental and estimated cell voltage for every current density point is similar to relative error obtained previously for curve-by-curve identification in the case of curves 1 and 31. Similar results are obtained for other curves of the DoE. The maximal error obtained between experimental and estimated polarization curves for one current density point is 1.7 % (obtained for the curve 31).

Table 5 -Comparison between error values obtained with results of identification by applying the new assumptions on parameters' dependences in comparison to previous results of curve-by-curve identification

The results presented in Table 5 point out that the results of identification by applying the new assumptions on parameters' dependences are really close to the results of curve-by-curve identification for presented curves. Moreover, as shown in table 4, identified parameter values appear to be consistent with the previous identified parameters (Table 3) using a curve-by-curve identification. It confirms that hypotheses which have been done in this identification work are relevant.

One of the advantages is that the identification time has been considerably reduced (around 30 minutes for 10 different starting points with this identification compared to 3 hours with curve-bycurve identification). Moreover, the identified values are less dependent on the quality of the experimental curves. That means that if there is an aberrant point or a lack of data in some parts of an experimental curve, this default will be smoothed.

Construction of parameter variation laws as a function of operating conditions

In this section, the objective is to determine variation laws as a function of operating conditions for each model parameter: J 0 (T,P), J n (T,P,HR air ) and R diff (T,P,HR air, HR H2 ) respectively starting from 8 J 0 , 16 J n and 31 R diff identified in the section 3.3.

For each parameter, an exploitation typical of the DoE methodology will be proposed, that means the determination of an empirical law assuming a linear dependence to the considered operating conditions. If possible, another exploitation is proposed in parallel by introducing a physical law. In practice, according to our knowledge, this is only possible for J 0 .

In this following, 50 starting points have been randomly drawn in the fixed variation domain of each parameter. For every variation law identification, calculations have lasted less than 10 minutes.

Variation laws for J 0 (T, P)

8 values for J 0 are available to carry out the variation law identification. Two laws are successively evaluated for J 0 : one purely empirically, the second one based on physical considerations.

Empirical variation law for J 0 (T, P):

Considering that J 0 depends on temperature and on pressure, an empirical variation law can be proposed:

. 2 ( , ?)

. 2_

1

[ :2 × ] ^:2 × $ _ :2 × ] × $ ( 14
)
Where J 0_mean , a J0 , b J0 and c J0 are parameters to be identified

The error ε J0 which is calculated in the objective function and which has to be minimized in the optimization process is based on a least square criterion:

C :2 D A . 2_
where J 0_ref is the reference exchange current density, meaning the previous identified J 0 which are presented in the part 3.3, and J 0_est is the estimated exchange current density calculated using the variation law parameters obtained in the identification process.

Table 6 gives a comparison of reference J 0 and estimated J 0 (empirical law case). The maximal relative error using the proposed empirical variation law is 4.9%.

T p J 0_ref (norm) J 0_est (norm)
Relative error (%) between J 0_ref and J 0_est Table 6 -Comparison of reference J 0 (previously identified) and estimated J 0 using the defined empirical variation law

The orders of magnitude (rough values) of variation law parameters are:

J 0_mean ≈ 5.10 -7 A/cm²; a J0 ≈ 2.10 -7 ; b J0 ≈ 1.10 -7 ; c J0 ≈ 4.10 -8
These coefficient values show that temperature and pressure have a similar impact on the parameter J 0 .

Physics-based variation law for J 0 (T, P):

A physics-based relation constructed with an Arrhenius law is given in references [START_REF] Labach | Caractérisation et modélisation de Piles à combustible et d'Electrolyseurs PEM à conditions opératoires variables, en vue de leur association[END_REF] and [START_REF] Labach | Steady-state Semi-empirical Model of a Single Proton Exchange Membrane Fuel Cell (PEMFC) at Varying Operating Conditions[END_REF]:

. 2 ( , ?) . 2_ A ? ? B `aO$ b A1 Bc

Where the reference exchange current density J 0_arr (A/cm²), the reaction order with respect to oxygen γ and the activation energy E act (J/mol) are parameters to be identified.

Table 8 gives a comparison of reference J 0 and estimated J 0 (physics-based law case). The maximal relative error using the proposed physics-based variation law is 7.3 %. It appears that the error is a slightly more important in this case than for the estimation using the empirical variation law.

t p J 0_ref (norm) J 0_est (norm)
Relative error (%) between J 0_ref and J 0_est The orders of magnitude of variation law parameters are: J 0_arr ≈ 2.10 -7 A/cm²; γ ≈ 0.8; E act ≈ 69.10 3 J/mol

Variation law for J n (T, P, RH air )

Considering that J n depends on temperature, pressure and air relative humidity, an empirical variation law can be proposed for the parasitic equivalent current density J n :

. 1 . 1_ 1 [ :1 × ] ^:1 × $ _ :1 × ℎe f :1 × ] × $ a :1 × ] × ℎe g :1 × $ × ℎe ( 17 )
Where J n_mean , a Jn , b Jn , c Jn , d Jn , e Jn and f Jn are parameters to be identified Table 9 gives a comparison of reference J n and estimated J n . The maximal relative error using the proposed empirical variation law is 6.7%. The orders of magnitude of variation law parameters are: J n_mean ≈ 1.10 -4 A/cm²; a Jn ≈ 3.10 -5 ; b Jn ≈ 1.10 -5 ; c Jn ≈ 3.10 -5 ; d Jn ≈ 6.10 -19 ; e Jn ≈ 8.10 -6 ; f Jn ≈ 3.10 -18 It appears that temperature, pressure and air relative humidity have a comparable impact on J n . It also shows that certain crossed terms can be neglected.

Variation law for R diff (T, P, RH air , RH H2 )

Considering that R diff depends on temperature, pressure and gas relative humidities, an empirical variation law can be proposed for diffusion resistance R diff : The orders of magnitude of variation law parameters are: R diff_mean ≈ 0.2 Ω.cm²; a Rdiff ≈ 2.10 -2 ; b Rdiff ≈ -9.10 -2 ; c Rdiff ≈ -4.10 -2 ; d Rdiff ≈ -3.10 -2 ; e Rdiff ≈ -3.10 -3 f Rdiff ≈ -5.10 -3 ; g Rdiff ≈ -2.10 -3 ; h Rdiff ≈ 3.10 -2 ; i Rdiff ≈ 2.10 -2 ; j Rdiff ≈ 2.10 -3

=MXV [ 6 × ] ^6 × $
It can be observed that temperature, pressure and gas humidities present a similar impact on R diff value. Moreover, influence of interaction terms between operating conditions appears also to be important and cannot be neglected. It highlights that the assumption to consider that R diff is dependent on every operating condition is coherent.

Reconstitution of the DoE curves using the identified variation laws in function of operating conditions

Using the previously-defined variation laws, it is possible to reconstitute all the DoE polarization curves. Two cases are possible: in the first case, the empirical variation law for J 0 is used whereas in the second case, the physics-based variation law for J 0 is used. Unchanged variation laws are used in both cases for J n and for R diff . Our objective is to quantify the errors by using these laws to simulate the DoE polarization curves. Figure 8 and Figure 9 present the experimental polarization curve and the reconstructed polarization curve using the defined variation laws respectively for curve 1 and 31. The results show that it is possible to reconstitute the DoE curves using the defined variation laws with a satisfactory precision. Results in the two cases are quite similar. The maximal error obtained between the experimental cell voltage and the estimated cell voltage at a given current density is close to 5.5% (obtained for curve 12 in both cases) and the mean error is lower than 1%.

Testing of the predictive behaviour of the model using the identified variation laws in function of operating conditions

In this part, first elements of predictive behaviour testing of the model are presented. The objective is to estimate the error when predicting a polarization curve that is taken outside the operating domain of the previous DoE.

As R ohm is measured with Electrochemical Impedance Spectroscopy in our approach, R ohm is an input data of the model, even this evaluation phase. In order to be fully predictive, it would be necessary to determine a variation law of R ohm as a function of operating conditions.

Curve 32 corresponding to T = 80°C, P = 1.7 bara, RH air = 19 %, RH H2 = 31 % is firstly proposed for this testing.

Figure 11 presents the identification results for curve 32. Two cases are considered: in the first case, J 0 empirical variation law is used whereas in the second case, J 0 physics-based variation law is used. In both cases, the defined empirical variation laws for J n and R diff are used. It appears that the prediction of the polarization curve in the case of curve 32 is quite coherent with the experimental polarization curve obtained. The maximal relative error between experimental and predictive curves is close to 8 % and the mean relative error is close to 2 %. The predictive polarization curve is close to the experimental curve at low current densities. However, the error between the two curves increases as the current density increases and is maximal at the higher current densities. This can be explained by the difficulty of modelling the diffusion phenomena by an empirical law as observed in Table 9.

Comparable results are obtained for estimation using J 0 empirical variation law and J 0 physicsbased variation law.

Curve 33 corresponding to T = 80°C, P = 0.65 bara, RH air = 12 %, RH H2 = 20 % is secondly studied in the same way. As for the previous curve, R ohm is considered known (R ohm is not estimated using a variation law).

Figure 11 presents the identification results for curve 33. As explained before, two cases are considered for J 0 variation law.

For this curve 33, the conclusions are similar to those for curve 32: the results show again that the model is able to predict with satisfactory precision cell voltage at low current densities. At higher current density, an "unhooking" appears between the experimental and estimated curves. Similar results are obtained for estimation using J 0 empirical variation law and J 0 physics-based variation law. However, evaluate the predictive behaviour of the model is limited with this curve 33 because only a restricted domain of current density were explored. In order to evaluate the predictive behaviour of the model, a curve-by-curve identification is carried out for curves 32 and 33. It is to say that a set of parameters (J 0 , J n , R diff ) is identified directly for each of these two curves. The results of these curve-by-curve identifications are compared to prediction results presented before and obtained using applying variation laws. Tables 10 andTable 11 present the results respectively for the curve 32 and the curve 33. These tables show that estimation of J 0 and J n gives quite acceptable errors. The error is more important for prediction of R diff . These results point out that it is necessary to work again on the definition and identification of variation laws particularly for diffusion phenomena, maybe on the parameter dependences on operating conditions.

Conclusion

In this paper, the complex behavior a PEMFC stack in function of various operating conditions was experimentally investigated by the Design of Experiments methodology. The various operating conditions explored in this DoE are representative of aeronautical conditions within an aircraft cabin. The experimented subatmospheric pressures are one of the highlights of our works.

After being discussed, these experimental results have been studied via a modelling approach of the polarization curve. Firstly, a curve-by-curve study was carried out. These first modelling results were exploited to evaluate if the dependence of each model parameter was strong or not to each involved operating condition. Thus, assumptions have been introduced concerning the parameter dependence to the operating conditions. The consequences were a simplification of the model (even if it remains complex) and a reduction of the identification time. Under these simplification hypotheses, the polarization curves were estimated with no significant loss of precision compared to the case of curve-by-curve identification.

Using the results of this first modelling step, variation laws for model parameters J 0 , J n and R diff have been constructed. The aim of defining these laws is to be able to have a predictive model in function of operating conditions. Our initial goal was to use and/or define physics-based laws describing this parameter dependence to the operating conditions. We were only able to use a physicsbased law for the activation current J o issued from previous works at the LAPLACE laboratory. A result is that this law appears also valid for subatmospheric pressures. The modelling with physicbased laws for J n and R diff , even for R ohm , remains an issue. Only empirical laws were proposed for J n and R diff .

The predictive behaviour of the model was tested on two curves taken outside the operating field of the DoE. The estimation of these two polarization curves seems promising. However, at high current densities, the error between the predictive curve and the real experimental curve is more important, illustrating that the diffusion phenomena modelling is not easy. For their modelling, a direct exploitation with a typical empirical law generally used with a DoE (assumption of a linear dependence to the operation conditions), as done in this paper, gives good results, but the limits in terms of precision have been quickly reached regarding its predictive capability. Future works will be devoted to correct this error and to improve variation law definition. It could be interesting to complete the carried out DoE with higher current densities in order to investigate more deeply the nonlinear behaviour of the diffusion phenomena, and perhaps in this way to improve the predictive character of our model.

More generally, the behaviour of a fuel cell with subatmospheric pressures remains also an issue: we have indeed observed in our works the big impact of pressure on the fuel cell performance and new exploratory tests would be necessary to go further.

Figure 1 -

 1 Figure 1 -Selected experimental polarization curves of the DoE: illustration of the influence of temperature and pressure, the gas humidities being the highest possible (around 26%)

Figure 3 -

 3 Figure 3 -Evolution of R ohm (measured by EIS) according to the current for selected operating conditions of the DoE -(a) illustration of the influence of temperature and pressure, the gas humidities being the highest possible (around 26%) -(b) illustration of the influence of gas humidities

Figure 4 and

 4 Figure4and Figure5give the experimental polarization curve (points in red), the estimated polarization curve (in blue) and the relative error at each current density between the experimental polarization curve and the estimated polarization curve, respectively for curves 1 and 31 of the DoE.

Figure 4 -Figure 5 -

 45 Figure 4 -Experimental polarization curve and estimated polarization curve (curve-by-curve identification)Curve 1: T = 80°C, P = 1.5 bara, RH air = 30 % (/"26%"), RH H2 = 30 % (/"26%")

Figure 8 -Figure 9 -

 89 Figure 8 -Experimental polarization curve and estimated polarization curve using the variation lawsCurve 1: T = 80°C, P = 1.5 bara, RH air = 30 % (/"26%"), RH H2 = 30% (/"26%")

Figure 10 -

 10 Figure 10 -Experimental polarization curve and estimated polarization curves using the different variation laws in function of operating conditions (2 cases for J 0 variation law) Curve 32: T = 80°C, P = 0.65 bara, RH air = 19 %, RH H2 = 31 %

Figure 11 -

 11 Figure 11 -Experimental polarization curve and estimated polarization curves using the different variation laws in function of operating conditions (2 cases for J 0 variation law) Curve 33: T = 80°C, P = 0.65 bara, RH air = 12 %, RH H2 = 20 %

Table 1 -

 1 Operating conditions applied for our Design of Experiments

	Operating factor	Levels	Number of level
	T	65 -80 (°C)	2
	P	0.8 -1 -1.3 -1.5 (bara)	4
	RH air	16 -26 (%)	2
	RH H2	16 -26 (%)	2

Table 3 -

 3 Results of curve-by-curve identification for every curve of the DoE

Table 7 -

 7 Comparison of reference J

0 (previously identified) and estimated J 0 using the physics-based variation law

Table 8 -

 8 Comparison of reference J n (previously identified) and estimated J n using the defined empirical variation law

  Where R diff_mean , a Rdiff , b Rdiff , c Rdiff , d Rdiff , e Rdiff , f Rdiff , g Rdiff , h Rdiff , i Rdiff and j Rdiff are parameters to be identified.Table9gives a comparison of reference R diff and estimated R diff . The maximal relative error using the proposed empirical variation law is 17.1%.

		t	p	hr air	hr H2	R diff_ref (norm)	R diff_est (norm)	Relative error (%) between R diff_ref and R diff_est
		1.00	1.00	1.00	1.00	0.6183	0.617	1.42
		1.00	1.00	-0.44	1.00	0.6317	0.628	4.68
		1.00	1.00	1.00	-0.44	0.6248	0.622	3.65
		1.00	1.00	-0.44 -0.44	0.6442	0.634	11.22
		-1.00	1.00	1.00	1.00	0.5977	0.597	1.06
		-1.00	1.00	-0.44	1.00	0.6038	0.598	6.98
		-1.00	1.00	1.00	-0.44	0.6017	0.596	6.51
		-1.00	1.00	-0.44 -0.44	0.6144	0.600	15.34
		-1.00	0.43	0.79	0.79	0.6109	0.616	5.51
		-1.00	0.43	-0.56	0.79	0.6212	0.633	14.92
		-1.00	0.43	0.79	-0.56	0.6175	0.629	13.84
		-1.00	0.43	-0.56 -0.56	0.6321	0.646	17.13
		1.00	0.43	0.79	0.79	0.6316	0.636	5.70
		1.00	0.43	-0.56	0.79	0.6493	0.657	9.38
		1.00	0.43	0.79	-0.56	0.6442	0.650	7.45
		1.00	0.43	-0.56 -0.56	0.6629	0.670	8.14
		1.00	-0.43	0.40	0.40	0.6772	0.673	4.48
		1.00	-0.43 -0.80	0.40	0.6947	0.696	1.95
		1.00	-0.43	0.40	-0.80	0.6882	0.690	2.62
		1.00	-0.43 -0.80 -0.80	0.6989	0.711	14.30
		-1.00 -0.43	0.40	0.40	0.6667	0.655	12.84
		-1.00 -0.43 -0.80	0.40	0.6779	0.678	0.60
		-1.00 -0.43	0.40	-0.80	0.6733	0.673	0.06
		-1.00 -0.43 -0.80 -0.80	0.6902	0.694	4.08
		1.00	-1.00	0.07	0.07	0.7077	0.700	8.04
		1.00	-1.00 -1.00	0.07	0.7336	0.721	13.71
		1.00	-1.00	0.07	-1.00	0.7275	0.716	12.54
		-1.00 -1.00	0.07	0.07	0.6918	0.684	8.08
		-1.00 -1.00 -1.00	0.07	0.7140	0.706	9.18
		-1.00 -1.00	0.07	-1.00	0.7078	0.701	7.41
		-1.00 -1.00 -1.00 -1.00	0.7227	0.719	3.71
	a 6 j 6	× ] × $ g 6 × $ × ℎe %* k 6	× ] × ℎe × ℎe × ℎe %* h 6	_ 6 × t × ℎe %* ℎ 6 × ℎe	f 6 × $ × ℎe × ℎe %*	( 18 )

Table 9 -

 9 Comparison of reference R diff (previously identified) and estimated R diff using the defined empirical variation law

Reference value (norm) Curve-by-curve identification Estimated value (norm) Prediction with the model Error between reference and estimated values (%)

	J 0	0.487	Empirical law: 0.499 Physical law: 0.503	Empirical law: 12 Physical law: 17
	J n	0.739	0.746	4.8
	R diff	0.621	0.593	27

Table 10 -

 10 Comparison of curve-by-curve identification and prediction results for curve 32

Reference value (norm) Curve-by-curve identification Estimated value (norm) Prediction with the model Error between reference and estimated values (%)

	J 0	0.427	Empirical law: 0.428 Physical law: 0.416	Empirical law: 0.25 Physical law: 9.7
	J n	0.686	0.674	8
	R diff	0.759	0.739	21

Table 11 -

 11 Comparison of curve-by-curve identification and prediction results for curve 33

26%, 26% 26%, 16% 16%, 26% 16%, 16% RHair, RHH2 26%, 26% 26%, 16% 16%, 26% 16%, 16% 26%, 26% 26%, 16% 16%, 26%

Acknowledgments

The authors would like to express their sincere thanks to all the members of the FUCHYA project, it is to say the IRT Saint Exupéry, ZODIAC AEROTECHNICS, SAFRAN, AIRBUS, the LAPLACE laboratory and to the ANR (Agence Nationale de la Recherche).