André De Palma
email: andre.depalma@ens-cachan.fr

Karim Kilani
email: karim.kilani@cnam.fr

Moshe Ben-Akiva

Gilbert Laond

Jen Lainé

Claude Lefèvre

Ordered choice probabilities in random utility models *

Keywords: Generalized Roy's identity, Logit, Ordered utilities, Order statistics, Probit, Random utility models. JEL classication: D11

We prove a general identity which states that any element of a tuple (ordered set) can be obtained as an alternating binomial weighted sum of rst elements of some sub-tuples. The identity is then applied within the random utility models framework where any alternative's ordered choice probability (the probability that it has a given rank) is expressed with respect to standard best choice probabilities. The logit and the logsum formulas are extended to their ordered choice counterparts. In a symmetric case, we compare for the probit and the logit, the surplus loss due to the withdrawal of a product with the damage due to the loss of a rank.

Introduction

Structural demand models have been very much used to estimate disaggregated data collected at the individual level. Probabilistic choice models, and in particular behavioral choice models (such as the logit, nested logit, mixed logit, probit, etc.) play here an important role. They are known as discrete choice models (DCM), which are very much used in many elds of economics and beyond (trade, industrial organization, labor economics, demography, urban and transport economics, marketing, etc. (see [START_REF] Mcfadden | Economic choices[END_REF]).

These tools describe and evaluate the behavior of individuals selecting their best alternative from a discrete set (usually nite and known) of alternatives (see Anderson, de Palma, and Thisse (1992)).

Most researchers in DCM and in particular in random utility models (RUM) focused their attention on the choice of the best alternative. The selected option, or best alternative, corresponds to the one yielding the highest realized utility. Theoretical formula for the probability of the best choice are crucial for the modeler mainly for building demand models, to study the theoretical properties of the demand system; simulation of choice probabilities are needed for estimating the parameters entering into the utility (relying on standard econometric package). The ranking of preferences has attracted much less attention. However, ranks matter. [START_REF] Anderson | Discrete Choice Theory of Product Dierentiation[END_REF] Quite often in surveys, elicitation of preference is achieved via stated preference questions, i.e. by asking the respondent to rank alternatives, such as by Netix. (2) Alternatively, Amazon, tour operators, doctors, health satisfaction surveys, investment advisers do not suggest the best alternative, but an array of ordered alternatives (see [START_REF] Greene | Modeling Ordered Choices: A Primer[END_REF]). ( 3) For a variety of reasons, the agents may not choose their best alternative since it may not be available (or known) when the decision was made (moreover, how much is lost in this case?). [START_REF] Ben-Akiva | Structure of Passenger Travel Demand Models[END_REF] Lack of full rationality or lack of full information provide other justications for studying lower-order choices.

We believe that rankings are more likely to matter with increasingly complex technologies, larger variety of goods and services, which make individual choices for the best alternative more dicult. A recent body of the literature has made signicant progresses in the direction of integrating lower-order or sub-optimal choices (see e.g. [START_REF] Berry | Dierentiated products demand systems from a combination of micro and macro Data: The new car market[END_REF]). These authors use rst and second-best choices revealed by respondents to survey questionnaires. They show that it is possible to obtain much better estimates in the context of the choice of dierentiated products and to reduce the burden associated to surveys. In the marketing literature, Marley and Louviere (2005) consider the extreme situation where respondents are asked to reveal both their best and their worse preferred alternative. They provide new insights on model parameter estimations obtained from this type of data.

Researchers who wish to incorporate more than just the optimal alternative in their study are interested to know how to derive expressions for the probabilities of full rank-ordered choice proles, where the alternatives are ordered from the best choice to the worst one. Unfortunately, this problem is dicult and explicit expressions exist in very few cases, such as for the rank-ordered logit (see [START_REF] Beggs | Assessing the Potential Demand for Electric Cars[END_REF]). However, even when such expressions exist (as in the logit case) they do not allow to recover simple expressions for single lower-order choices. This is mainly due to intricate combinatorics inherent to these problems. Some researchers have added new assumptions to the main axioms of standard discrete choice approach, to get for example, expressions for the worst choice probabilities (see [START_REF] Marley | Some probabilistic models of best, worst, and bestworst choices[END_REF]). Other approaches, such as the ordered probit model, which use threshold value, has also been proposed. Ordered choice modeling, mixed and hierarchical models include dierent ways to introduce thresholds, which represent one natural manner to model ranking probabilities. Finally, an alternative solution is the use of EM algorithms, Monte-Carlo simulation techniques exploiting Gibbs sampler, al-lowed by faster computer's processors in order to deal with lower-order choice models (see [START_REF] Train | Discrete Choice Methods with Simulation[END_REF]). This solution presents an interest for estimation, but is bounded to be slow for large choice sets, and it remains less attractive from a theoretical point of view.

In this paper, we adopt a very conservative view. We wish to keep the standard discrete choice random utility approach and study the ranking probabilities and welfare of lower ranked alternatives. The simple analysis of the logit model suggests some research agenda. It shows that the probabilities of ranks can be expressed as a sum of products of standard logit models.

However, these expressions seem very much related to the strong hypothesis behind the logit (i.i.d. and Gumbel error terms). Extension to the probit, even with i.i.d. normal distribution appears to be an impossible task.

Moreover, the standard logit analysis remains totally silent about the welfare properties of lower ranks, while the welfare of the best choice is well known since 40 years ago (see the derivation of the nested logit and logsum by Ben-Akiva (1973)).

The literature remains silent a fortiori for any other RUM. These questions will be addressed in this paper. Moreover, the relation between choice probabilities and welfare given by Roy's identity, is shown to remain true for lower ranks. Clearly, these results have to make abstraction of the error terms underlying the construction of the RUM. We proceed more generally with a system of choice probabilities. A rather long detour to order statistics will be required. The main mathematical background beyond this paper is drawn from order statistics. In our setting, we focus on the r-th order statistics, i.e. on r-th highest value (standard order-statistics consider the k-th lowest values). We focus on the relation with the maximum (the reasoning for the minimum would be the same).

The maximum plays a key role in economics, since the maximum of random variables corresponds to the alternative chosen by the individual in a RUM, and since the expected maximum utility is the consumer surplus (when there are no income eects). We start with a deterministic approach, with ordered sets. We wish to compare elements in dierent order and in dierent ordered sets. To xed idea consider an example. Assume that a student is ranked 5-th in a class of 20 students. S/he should be encouraged if 6 months later, s/he turns out to be 4-th out of 20. However, s/he should also be encouraged if the size of the class increases from 20 to 25 students, while maintaining the 5-th rank. However, it is harder to decide whether or not s/he should be more encouraged in the rst situation than in the second.

To address this question, we start with a deterministic version of the triangular inequalities, used in order statistics with random variables (see e.g., [START_REF] Arnold | A First Course in Order Statistics[END_REF]). Doing that, we are able to recover a deterministic version of the triangular inequality removing the restrictions on the correlation patterns. This approach will play a key role in this paper to obtain results valid for random utility model with any correlation patterns. Our rst result (Lemma 1) relates the gain from higher rank with the gain from the same rank in larger set. This identity relates ranks and size. Such result is standard for random variables in order statistics (triangular identity). Here it is shown to be valid for: (a) any element (not only for random variables but also for real numbers); and (b) set of elements not necessarily ordered. In this case, elements in the choice set are denoted by labels (which play the role of the rank). This deterministic triangular identity (Lemma 1) corresponds to a recurrence which can be solved. The main result of this paper (Theorem 2) provides the element of order r in a set C as a function the highest order element in dierent subsets of C appropriately chosen. Then, we show (Theorem 3) that the same identity is valid for any operator T (•). When T (•) is the expectation operator, and when elements are random utilities, Corollary 6 shows that any expectation of the r-th ordered utility can be expressed as expectations of maximum utility within appropriately chosen subsets. Similarly, Theorem 6 allows to compute the probability of rank r as a function of the standard (best) choice probabilities. We then provide a new version of Roy's identity, which is valid for any rank, and not, as in the standard case, for the best choice probabilities. To sum-up, our setting allows to generate the following economic results, valid for any random utility models (with the nested logit and the GEV as special case): (1) Relation between ranking probabilities and maximum probabilities. (2) Welfare of the ordered choice as a function of best-choice expected utility. (3) Relation between ranking probabilities and welfare of ordered choices, extending Roy's identity to ordered random utility models.

The paper is organized as follows. Section 2 introduces the setting for tuples of ordered objects and derived our ground results, Theorem 2, where an identity is obtained. Section 3 applies these results in the RUM framework where standard results appearing in the order statistics literature are derived as corollaries of Theorem 3. We also derive a relation between expectation of order utilities and expectation of maximum utilities. Choice probabilities in the RUM context are introduced in Section 4 where a new identity relating any order choice probability to best choice probabilities is obtained. It allows to derive a generalization of Roy's identity, valid for any welfare of expected ranked utility (Theorem 8). Application to the standard multinomial logit, to GEV models, and an illustrative example in the symmetric case are discussed in Section 5.

Preliminaries results

2.1

The main theorem Consider an n-tuple (ordered set, list) l n ≡ (l n (1) , . . . , l n (n)) of n elements of a vector space (real or complex numbers, real matrices, random variables, etc.), n ≥ 2. Let l p ≡ (l p (1) , . . . , l p (p)) denotes a p-tuple of p dierent elements drawn from l n , 1 ≤ p ≤ n, assuming that the order of appearance of the elements in l p remains the same as in l n . We denote by L p the set of such p-tuples, whose cardinality is n p .

Note that any l n-1 of L n -1 can be obtained by dropping a single element

l n (k) from l n , 1 ≤ k ≤ n. For convenience, such (n -1)-tuple is instead denoted by l [k] n , with r-th order element l [k] n (r) , 1 ≤ r ≤ n -1. Our rst result is a technical lemma.
Lemma 1. The following recurrence relation holds:

n k=1 l [k] n (r) = rl n (r + 1) + (n -r) l n (r) , 1 ≤ r ≤ n -1. (2.1) Proof. When building l [k] n , if the dropped element from l n has a rank k in l n such that k ≤ r (resp. k > r), 1 ≤ r ≤ n -1, the r-th order element of l [k]
n matches the (r + 1)-th (resp. r-th) order element of l n . This observation yields the following equation

l [k] n (r) = 1 {k≤r} l n (r + 1) + 1 {k>r} l n (r) ,
where 1 {•} is the indicator operator. Summing up over k both sides of the above equation, we obtain

n k=1 l [k] n (r) = n k=1 1 {k≤r} l n (r + 1) + n k=1 1 {k>r} l n (r) .
Using the fact that:

n k=1 1 {k≤r} = n -n k=1 1 {k>r} = r, Eq. (2.1) is obtained.
In what follows, a m q ≡ (-1) m-q m q will denote alternating binomial coecients, where m and q are two integers such that 0 ≤ q ≤ m. The cornerstone result of this paper is now stated. Theorem 2. The following identity, relating the r-th order element of the tuple to rst elements of some of its sub-tuples, holds:

l n (r) = n p=n-r+1 a p-1 n-r lp∈Lp l p (1) , 1 ≤ r ≤ n. (2.2)
Proof. Identity (2.2) states that l n (r) = R n (r), where R n (r) is the the RHS of (2.2): R n (r) ≡ n p=n-r+1 a p-1 n-r lp∈Lp l p (1). For any n ≥ 2 and r = 1, the sum contains only one term. Using the fact that a n-1 n-1 = 1 and that L n contains the single tuple l n , we obtain l n (1) = R n (1). For r = 2, the rst sum runs from p = n -

1 to p = n. Since a n-1 n-2 = -(n -1), we get R n (2) = n k=1 l [k] n (1) -(n -1) l n (1) . Eq. (2.1) with r = 1 implies n k=1 l [k] n (1) = l n (2) + (n -1) l n (1)
, so that the above equation yields l n (2) = R n (2). Therefore, (2.2) is true for any n ≥ 2 and for r = 1, 2. For n = 2, note that (2.2) is true for all the possible ranks (r = 1, 2). For n ≥ 3, by induction, assume that (2.2) is valid for any (n -1)-tuple and all possible rank r verifying 1 ≤ r ≤ (n -1). This implies that, for any k such that 1 ≤ k ≤ n, we have

l [k] n (r) = n-1 p=n-r a p-1 n-1-r lp∈L [k] p l p (1) , where L [k] p is the set of all p-tuples, 1 ≤ p ≤ n -1, drawn from l [k]
n in the same manner as the tuples drawn from l n , i.e. by keeping the same order for the selected elements. We prove that (2.2) is true for n and any possible rank r veryng 1 ≤ r ≤ n. As seen above, (2.2) is true for r = 1, 2. By induction on r, assume that R n (r -1) = l n (r -1) for a given rank r,

1 ≤ r ≤ n -1. The tuples of L [k]
p consists of tuples of L p but those containing the element l n (k), their set being denoted hereafter by L ←k p , must be subtracted. This observation allows us to rewrite the last above equation as

l [k] n (r) = n-1 p=n-r a p-1 n-1-r   lp∈Lp l p (1) - lp∈L ←k p l p (1)   .
By summation over all k, we get

n k=1 l [k] n (r) = n-1 p=n-r a p-1 n-1-r   n lp∈Lp l p (1) - n k=1 lp∈L ←k p l p (1)   .
Within the double summation appearing into the brackets, any tuple l p of L p , n -r ≤ p ≤ n -1, is accounted for as many times as its number of elements, i.e. p times. Consequently, the following expression is obtained

n k=1 l [k] n (r) = n-1 p=n-r (n -p) a p-1 n-1-r lp∈Lp l p (1)
.

Note that the summation of the RHS can be extended to n, obtaining

n k=1 l [k] n (r) = n p=n-r (n -p) a p-1 n-1-r lp∈Lp l p (1)
.

It is worth to rewrite it as

n k=1 l [k] n (r) = rR n (r + 1) + n p=n-r (n -p -r) a p-1 n-1-r lp∈Lp l p (1)
.

Note that for p = n -r, the corresponding term of the sum of the RHS is null, allowing to rewrite the above expression as

n k=1 l [k] n (r) = rR n (r + 1) + n p=n-r+1 (n -p -r) a p-1 n-1-r lp∈Lp l p (1)
.

Using an elementary property of binomial coecients, the last equation be-

comes n k=1 l [k] n (r) = rR n (r + 1) + (n -r) n p=n-r+1 a p-1 n-r lp∈Lp l p (1) , yielding n k=1 l [k] n (r) = rR n (r + 1) + (n -r) R n (r) .
By induction, we have assumed that R n (r) = l n (r), so that

n k=1 l [k] n (r) = rR n (r + 1) + (n -r) l n (r) .
Comparing the above equation with Eq. (2.1), we get l n (r + 1) = R n (r + 1).

The above theorem, which provides an identity relating tuple's elements to rst order elements of sub-tuples, is now applied in a RUM framework.

We derive new identities among ordered utilities, their distribution, their expectation, etc., and then derive an identity for ordered choice probabilities.

Application to random utility models

There are n≥ 2 objects of choice referred by their labels. The set of n labels of the objects of choice (alternatives) is denoted by C. In a RUM framework, the utility U i derived by an individual selecting the object with label i (alternative i), i ∈ C, is modeled as a real-valued random variable dened on a common probability space (Ω, F, P ). For one particular realization {U i (ω)} i∈C , ω ∈ Ω, of the nite set of random utilities {U i } i∈C , the alternatives referred by their label, can be rankordered according to a decreasing order of magnitude of their associated utilities. The top rank alternative has the highest utility while the bottom rank alternative has the lowest utility. When two or more alternatives have the same utility, the order among them can be chosen arbitrary, e.g. using their alphabetic order.

Denote by (U 1:n (ω) , ..., U n:n (ω)) the n-tuple of ordered utilities, where U 1:n (ω) ≡ max i∈C U i (ω) , ..., U n:n (ω) ≡ min i∈C U i (ω). We refer to U r:n , the random variable with realizations U r:n (ω), ω ∈ Ω, as the r-th ordered utility,

1 ≤ r ≤ n.
The results which follows are the most general and use operators denoted by T (•) which are mappings from the vector space of real random variables dened on (Ω, F, P ) into another vector space. This general setting allows us to apply Theorem 2 to ordered utilities, to their CDFs, their expectations, etc. Theorem 3. For any operator T (•) mapping the space of real random variables into another vector space, the following identity holds among ordered utilities:

T (U r:n ) = n p=n-r+1 a p-1 n-r Ap⊆C T max k∈Ap U k , 1 ≤ r ≤ n, (3.1)
where A p denotes subsets of C with p cardinality.

Proof. Dene the n-tuple l n ≡ (T (U

1:n ) , • • • , T (U n:n ))
, where its r-th order element is T (U r:n ), 1 ≤ r ≤ n. Thanks to Theorem 2, we have

T (U r:n ) = n p=n-r+1 a p-1 n-r lp∈Lp l p (1)
,

where l p (1) is the rst element of a p-tuple l p drawn from l n . Let L p (1) ≡ {l p (1) : l p ∈ L p } be the set of elements of l n appearing as rst elements of some l p ∈ L p . Note that L p (1) has n p elements. Recall that the ordered utilities used in l n verify:

U 1:n ≥• • • ≥ U n:n . Therefore, L p (1) = T max k∈Ap U k : A p ⊆ C .
Note that there are n p subsets A p . Therefore, the sum lp∈Lp l p coincides with Ap⊆C T max k∈Ap U k . Putting together the above observations, the required equation is obtained.

The above theorem can be applied to derive an identity among ordered utilities.

Corollary 4. The ordered utilities verify the following identity:

U r:n = n p=n-r+1 a p-1 n-r Ap⊆C max k∈Ap U k , 1 ≤ r ≤ n. (3.2)
Proof. Apply Theorem 3 using as T (•) operator the identity of the space of real random variables dened on (Ω, F, P ), to obtain the required identity.

Eq. (2.1) implies:

k∈C U r:n -U [k] r:n r = U r:n -U r+1:n =, 1 ≤ r ≤ n -1, (3.3)
where U

[k] r:n denotes the r-th order utility when alternative k is dropped from the choice set C. The LHS of (3.3) is an average utility gain due to a larger choice. It coincides with the RHS which measures the utility gain derived from a better ranking.

Recurrence relations and identities appearing in the order statistics literature apply to the distributions of the order statistics rather than the order statistics themselves (or their realizations). In this vein, let F i (x) ≡ P (U i ≤ x), x ∈ R, be the marginal CDF of U i , i ∈ C. Under assumptions kept the most general for now, we show that a standard result of the order statistics literature can be seen as a corollary of Theorem 3.

Corollary 5. The CDFs of the ordered utilities verify, for any real x, the following identity:

P (U r:n ≤ x) = n p=n-r+1 a p-1 n-r Ap⊆C P max k∈Ap U k ≤ x , 1 ≤ r ≤ n. (3.4)
Proof. Let T (•) be the following operator T (X) = P (X ≤ x) , where X is any random variables dened on (Ω, F, P ). It maps the space of real random variables into [0, 1]. Application of Theorem (3) with such operator yields Eq. (3.4).

We now apply Theorem 3 to derive an identity among the expectations of the ordered utilities.

Corollary 6. The expectations of the ordered utilities, µ r:n ≡E (U r:n ), when they exist, verify the following identity:

µ r:n = n p=n-r+1 a p-1 n-r Ap⊆C E max k∈Ap U k , 1 ≤ r ≤ n. (3.5)
Proof. We apply Theorem 3 with T (•) chosen as being the expectation operator: T (X) = E (X) , which maps the space of real random variables into the real line.

Note that no linearity is required for the T (•) operator. Similar corollaries can also be stated for any m-th order moment (taking T (X) = E (X m )), or even for variances of the order statistics (taking T (X) = Var (X) , where Var (•) is the variance operator). [START_REF] Ben-Akiva | Structure of Passenger Travel Demand Models[END_REF] Ordered choice probabilities 4.1

Choice probabilities

Assumptions about the probability distribution of the random utilities {U i } i∈C are slightly strengthened. Its CDF given by P (U i ≤ x i , i ∈ C) is assumed to be absolutely continuous with respect to Lebesgue measure over R n . One consequence of absolute continuity is that ties among utilities occur with zero probability, so that almost surely, there is a strict ranking among the alternative utilities.

Consider the event (U r:n = U i ) where it occurs that i ∈ C corresponds to the r-th order choice in the choice set C, 1 ≤ r ≤ n. Order choice probabilities P r:n (i) are dened, for i ∈ C , as being the probability of this event P r:n (i) ≡ P (U r:n = U i ) . The sequence (indexed by the alternatives) of events {(U r:n = U i )} i∈C form a partition of Ω (up to a null-measure set), implying that:

i∈C P r:n (i) = 1, 1 ≤ r ≤ n.
The same property prevails for the sequence (indexed by the ranks) of events {(U r:n = U i )} 1≤r≤n , so that:

n r=1 P r:n (i) = 1, i ∈ C.
An important topic in RUMs is the derivation of an expression for the best choice probabilities P Ap (i)

≡ P max k∈Ap U k = U i , i ∈ A p , A p ⊆ C.
Recall they can be obtained by performing the following integration (see e.g. Anderson, de Palma, and Thisse (1992))

P Ap (i) = ˆ+∞ -∞ ∂P (U i ≤ x i , i ∈ A p ) ∂x i x i =x dx, (4.1)
where the symbol x i =x means that the partial derivative has to be taken at

x i = x, i ∈ A p .
In the logit case where alternative utilities are independent Gumbels, a closed form recalled in Section 5.1 is obtained. This is not usually the case for any RUM, as for the probit where utilities are normally distributed, requiring simulation techniques exposed in [START_REF] Train | Discrete Choice Methods with Simulation[END_REF], in order to compute the integral (4.1).

Identity for choice probabilities

Theorem 3 is now adapted to derive a new identity for order choice probabilities.

Theorem 7. For any RUM, for any alternative i ∈ C, the ordered choice probabilities are related to best choice probabilities by:

P r:n (i) = n p=n-r+1 a p-1 n-r {i}⊆Ap⊆C P Ap (i) , i ∈ C, 1 ≤ r ≤ n. (4.2)
Proof. Consider the operator given by T (X) = P (X = U i ) , mapping the set of random variables X dened on (Ω, F, P ) into [0, 1]. Application of Theorem 3 with such operator yields

P (U r:n = U i ) = n p=n-r+1 a p-1 n-r Ap⊆C P max k∈Ap U k = U i .
For any subset A p verifying i / ∈ A p , thanks to the absolute continuity of the distribution of the utilities {U i } i∈C , we have P max k∈Ap U k = U i = 0, implying the required identity.

Williams-Daly-Zachary theorem extended

Additive random utility models (ARUMs) assume that utilities have an additive form:

U i = v i + ε i , i ∈ C
, where v i is the systematic part of the utility and ε i is the random error term. Utilities are assumed here to have nite expectation.

The Williams-Daly-Zachary theorem (see e.g. [START_REF] Mcfadden | Econometric Models of Probabilistic Choice[END_REF]) states that the derivative with respect to v i of the expected maximum utility (rst-order moment) within a choice subset A p , allows to recover the best choice probabilities

∂E max k∈Ap U k ∂v i = P Ap (i) , i ∈ A p . (4.3)
More generally, we prove that the derivation of the expected order utilities µ r:n with respect to v i allows to recover the order choice probabilities.

Theorem 8. For any ARUM where utilities have nite expectation, we have:

∂µ r:n ∂v i = P r:n (i) , i ∈ C, 1 ≤ r ≤ n. (4.4)
Proof. Deriving Eq. (3.5), with respect to v i , we obtain

∂µ r:n ∂v i = ∂ ∂v i n p=n-r+1 a p-1 n-r Ap⊆C E max k∈Ap U k .
Inverting the derivation and the sum signs of the RHS, we get Thanks to Identity (4.2), the RHS of the above equation is P r:n (i), obtaining the required identity.

5 Applications

The multinomial logit model

The multinomial logit model (MNL) is an ARUM with independent Gumbel distributed utilities with marginal CDFs given, for any real x, by

P (U i ≤ x i ) = exp -e v i -x i -γ , i ∈ C, (5.1)
where γ is Euler's constant. Note that v i is the expectation of the utility: (5.3)

E (U i )= v i .
Recall that the best choice probabilities also have a closed form given, for i ∈ A p , by the logit formula

P Ap (i) = e v i k∈Ap e v k , i ∈ A p .
(5.4) Application of Eq. (4.2) of Theorem 7 combined with the above logit formula yields a generalized logit formula for r-th order choice probabilities

P r:n (i) = n p=n-r+1 a p-1 n-r {i}⊆Ap⊆C e v i k∈Ap e v k
, 1 ≤ r ≤ n.

(5.5)

For example, the second-best logsum is given by

µ 2:n = i∈C ln k∈C-{i} e v k -(n -1) ln k∈C e v k , (5.6)
while the second-best logit choice probabilities are given by (5.8)

In accordance with formulas of the rank-ordered logit (see [START_REF] Beggs | Assessing the Potential Demand for Electric Cars[END_REF]), every term of the sum of the RHS corresponds to the joint probability that j is the best choice and i is the second best in C.

On the other hand, the expected utility corresponding to the worst choice is given by

µ n:n = n p=1 (-1) p-1 Ap⊆C ln k∈Ap e v k .
(5.9)

It provides a useful benchmark for welfare measures. Moreover, using (5.5), the worst choice probabilities where minimum utility is achieved can be written as an alternating sum of logit expressions

P n:n (i) = n p=1 (-1) p-1 Ap⊆C e v i k∈Ap e v k , i ∈ C.
(5.10)

This formula has been derived by de Palma, Kilani, and Laond (2013) making use of the Inclusion-Exclusion principle in probability theory.

Correlation among utilities

The MNL can be extended to allow for correlation among the random utilities which are still assumed to have Gumbel margins with CDFs given by Eq.

(5.1). Correlation is introduced via absolutely continuous copulas (for an introduction to copulas, see e.g. Nelsen (1999)) Θ (•) which are CDF functions dened over the unit n-cube [0, 1] n . The multivariate CDF of the random utilities is becoming

P (U i ≤ x i , i ∈ C) = Θ exp -e v 1 -x 1 -γ , .
.., exp -e vn-xn-γ . (5.11) We focus our attention to the class of extreme value copulas (EVC, also referred to as max-stable copulas), which verify the following property (see [START_REF] Salvadori | Extremes in Nature: An Approach Using Copulas[END_REF]

, p. 192) Θ φ λ 1 , ..., φ λ n = C λ (φ 1 , ..., φ n ) , ∀λ > 0.
(5.12)

Under the above assumptions, the maximum utility within a choice subset A p ⊆ C has a CDF which can be written as

P max k∈Ap U i ≤ x = exp -G Ap e -x , (5.13)
where G C (referred to as tail dependence functions) is dened as G C ≡ -ln Θ (exp (-e v 1 ) , ..., exp (-e vn )) ,

while any other G Ap is obtained by setting v k = -∞ for all k / ∈ A p in the RHS of the above expression.

A consequence of Eq. (5.13) is that the maximum utility is also Gumbel distributed with expected value given by E max k∈Ap U k = ln G Ap .

(5.15)

The above formula provides a generalization of the logsum formula (5.2) to the GEV framework. Thanks to the Williams-Daly-Zachary theorem, best choice probabilities can be derived using the following (5.17) Moreover, thanks to Theorem 7 and using (5.16), the r-th order choice probabilities are given by

P Ap (i) = ∂ ln G Ap ∂v i , i ∈ A p . ( 5
P r:n (i) = n p=n-r+1 a p-1 n-r {i}⊆Ap⊆C ∂ ln G Ap ∂v i , i ∈ C, 1 ≤ r ≤ n (5.18)
Example 1 The Gumbel-Hougaard (or logistic) family copula, the only Archi-median family of max-stable copulas, has the following form (5.20)

Θ (φ 1 , ..., φ c ) = exp    - i∈C (-ln φ i ) θ 1 θ    , θ ≥ 1.
As θ = 1, they coincide with the logit order choice probabilities given by (5.5).

Example 2 Let {C 1 , ..., C m } be a partition of the choice set C into m subsets or groups of alternatives, and construct a mixture of a Gumbel-Hougaard copula given by (5.19) and the independent copula in the following way . Therefore, order choice probabilities of any alternative belonging to group g are given, for 1 ≤ r ≤ n, and i ∈ C g , The two gures look very similar. Several comments are in order. The expected maximum utility is increasing and concave while the expected minimum utility is decreasing and convex. For a given number of alternatives, the r-th order utility decreases with r and the utility loss when losing one rank is rst decreasing and then increasing with the rank. This means that the loss from the penultimate to the last rank is signicant. The expected utility of the r-th rank (computed from the top) increases with c while the expected utility of the r -th rank (computed from the bottom) decreases with n, consistently with the triangular condition (5.23). This means in particular that the penultimate for n alternatives is better than the penultimate for n + 1 alternatives. Finally, about the Julius Caesar's quote I had rather be

P

  ∂E max k∈Ap U k ∂v i .For any A p such that i / ∈ A p , since ∂E max k∈Ap U k /∂v i = 0, the above equation can be written as ∂µ rAp (i) .

Ae v k , 1 ≤

 1 closure property of the Gumbel distribution ensures that the CDF of max k∈Ap U k , the maximum utility among the alternatives of A p , remains a Gumbel with an expectation given by E max k∈Ap U k = ln k∈Ap e v k , A p ⊆ C, (5.2) the celebrated logsum formula often used as a welfare measure in empirical works.Applying Identity (3.5) of Corollary 6 and using the logsum formula (5.2), a generalized logsum formula is obtained for the expected r-th order utility within the MNL framework µ r:n = r ≤ n.

P 2

 2 the above expression yields an alternative formP 2:n (i) = j∈C-{i} e v j k∈C e v k e v i k∈C-{j} e v k .

  .16) They correspond to the GEV probabilities obtained by McFadden (1978). Corollary 6 allows a generalization of (5.15) to any order choice µ r:n = ln G Ap , 1 ≤ r ≤ n.

( 5 .= k∈Ap e θv k 1 θ

 51 19)Note that the boundary case θ = 1 coincides with the product (independent) copula: Θ (φ 1 , ..., φ c ) = i∈C φ i . The associated tail dependence functions are G Ap . Order choice probabilities for such model are given by P r:n (i) = n p=n-r+1 a p-1 n-r {i}⊆Ap⊆C e θv i k∈Ap e θv k , i ∈ C, 1 ≤ r ≤ n.

Θ (φ 1

 1 , ..., φ c ) = exp

( 5 .

 5 21)Tail dependence functions are given by G Ap =

Figure 5 . 1 :Figure 5 . 2 :

 5152 Figure 5.1: Expected utility vs. rank order choice: logit

the rst in this village than second in Rome, we see that the answer depends on the relative sizes of Rome and the village. For example, the rst among six alternatives is better than the second among fteen.

by Ben Akiva (1973).

5.3

The symmetric case

In the symmetric case, arising for example with i.i.d. utilities, using (2.1), the following recursion rule is obtained

(5.23)

Only the computation of the expected best utilities (µ 1:2 , . . .,µ 1:n ) are needed, the remaining expectations are then computed iteratively from (5.23), which can also be written as

(5.24)

Hence, the expected damage due to a loss in one rank is higher than the one due to the loss of an alternative.

We consider a symmetric MNL where (utilities are centered) µ 1:n = ln n.

We also consider the case of i.i.d. normal utilities with zero mean and variance π 2 /6, in order to get the same variance as for the Gumbel case. Note that for the normal case, explicit values of µ 1:n can be obtained for n ≤ 5. 1 For larger values of n, numerical integration is needed.

Expected order utilities are computed using a spreadsheet for the Gumbel case and the R statistical software for performing numerical integration for the normal case. The results are displayed in the two gures below. The upper curve represents the expected maximum utility, which is used to deduce the lower gures representing expected lower-order utilities.

1 For the standard normal distribution, µ 1:2 = π -1/2 ; µ 1:3 = 1.5π -1/2 ; µ 1:4 = 6π -3/2 tan -1 √ 2; µ 1:5 = 15π -3/2 tan -1 √ 2 -2.5π -1/2 (see [START_REF] Arnold | A First Course in Order Statistics[END_REF]).


Mikaël Dumortier
email: mikael.dumortier@gmail.com

José Sanchez

Michel Keddam

Olivier Lacroix

Energy transport inside a three-phase electrode and application to a proton conducting solid oxide electrolysis cell

This work focuses on the modelling of thermal processes inside a planar high temperature steam electrolyzer that use cermets as electrodes. While the continuity equation for mass and charge have been demonstrated in a previous publication, energy balance for thermal transfers inside the electrode assembly is established via a control volume method. A non-dimensional number is built from different criterion used in the literature in order to validate the local thermal equilibrium assumption (LTE) inside the porous electrodes. A parametric analysis is carried out on a proton-conducting solid oxide electrolysis cell in galvanostatic mode. The results show that the heat sources are mainly ohmic and that their locations are not dependent on inlet current and inlet velocity of gases. This observation allows us to build an original thermal resistance network in order to analytically evaluate the temperature inside each component of the cell. This modelling strategy reduces computation time, allows reverse physical analysis and gives a precise estimation on the maximum temperatures attained in the components of the cells.

Nomenclature c p

Specific heat capacity at constant pressure (J K -1 kg -1 )

d p Pore diameter (m)
e Width of the electrode (m)

F Faraday constant (C mol -1 )
h g Specific enthalpy (J K -1 kg -1 )

h th(i-j) Heat transfer coefficient between phase i and j (W K -1 m -2 ) High temperature electrolysis of water by using proton conducting ceramic membranes is an innovative and attractive process for producing pure hydrogen. This process, which can be carried out without noble catalysts, consumes less electricity than classical low temperature electrolysis in aqueous medium. The electrolytic membrane reactor is composed basically of an anode and a cathode compartment, separated by a dense solid oxide membrane made of perovskite. This electrode assembly is commonly referred as the Positive electrode/Electrolyte/Negative electrode system (PEN).Water vapour is fed into the anode compartment and hydrogen is produced at the cathode with the recombination of protons conducted through the electrolyte (Fig. 1). Thus it is not necessary to separate water and hydrogen in contrast to oxygen conducting high temperature electrolysis cells. In both compartments, the electrochemical reaction takes place inside cermets at high temperature (500 to 600 °C) and pressure (10 to 20 bar) [START_REF] Mingyi | Thermodynamic analysis of the efficiency of hightemperature steam electrolysis system for hydrogen production[END_REF]. Cermets are porous composite materials made of a mixture of solid oxides and metal: the metal ensures electron transport; protons are transferred through solid oxides and gases flow through porosity (voids). The chemical reaction occurs at the junction of the three phases: metal, oxide and gas, called the triple phase boundary.

H

Fig. 1 -Principle of a proton conducting solid oxide electrolysis cell (SOEC).

The study of transport of transfer phenomena inside those electrolyzers is closely related to the study of solid oxide fuel cells (SOFC) and anion conducting solid oxide electrolysis cells. They all operate under hard operating conditions (1-10 bars, 500-1000°C), uses gases as reactants and porous media as electrodes. As a consequence, detailed in situ measurements during operation are very difficult to perform in those devices. Therefore, numerous numerical models for fuel cells and electrolyzers have been developed in the past two decades in order to acquire by calculation the operating conditions inside such membrane reactors [START_REF] Pramuanjaroenkij | Mathematical analysis of planar solid oxide fuel cells[END_REF][START_REF] Ferguson | Three-dimensional numerical simulation for various geometries of solid oxide fuel cells[END_REF][START_REF] Kakaç | A review of numerical modeling of solid oxide fuel cells[END_REF][START_REF] Andersson | Review on modeling development for multiscale chemical reactions coupled transport phenomena in solid oxide fuel cells[END_REF][START_REF] Soukup | Comparison of Wicke-Kallenbach and Graham's diffusion cells for obtaining transport characteristics of porous solids[END_REF][START_REF] Hussain | Mathematical modeling of planar solid oxide fuel cells[END_REF][START_REF] Demin | High temperature electrolyzer based on solid oxide coionic electrolyte: A theoretical model[END_REF][START_REF] Deseure | Theoretical optimisation of a SOFC composite cathode[END_REF]. Among the topics studied by modelling, temperature management inside high temperature electrolyzers and fuel cells is a key issue. Indeed, accurate temperature prediction within solid oxide cells is essential for calculating and optimizing the overall

Code de champ modifié

cell performance since the material properties, chemical kinetics, and transport properties depend on temperature [START_REF] Andersson | Review on modeling development for multiscale chemical reactions coupled transport phenomena in solid oxide fuel cells[END_REF].

Moreover, the operating temperature determines the nature of charge carriers in the electrolyte. For example, proton conducting solid oxide cells operate at below 600°C, while anion conducting oxide cells operate at above 600°C. A change in temperature can modify the type of charge carriers inside the electrolyte, hence switching the technology of the electrolyzer from one to another [START_REF] Slodczyk | Structural modifications induced by free protons in proton conducting perovskite zirconate membrane[END_REF].

Furthermore, strong temperature gradients can also cause thermo-mechanical degradation, such as delamination in the PEN, whose harmful incidence on the performance of the system was investigated by Xinfang and Xingjiang in 2010 [START_REF] Jin | Computational fluid dynamics analysis of solid oxide electrolysis cells with delaminations[END_REF]. Electrode material can also be damaged under high temperature gradients or thermal fatigue, allowing the hydrogen to mix with oxygen. This can be very dangerous as the auto-ignition temperature of hydrogen is close to the operating temperature of these electrolyzers [START_REF] Thiessen | The autoignition in air of some binary fuel mixtures containing hydrogen[END_REF]. Therefore, temperature monitoring is an important issue not only for the electrolyzer performance, but also for safety considerations.

Among the numerous models available to simulate the working operation of such membrane reactors, some consider the process as isothermal [START_REF] Izzo | Modeling of gas transport through a tubular solid oxide fuel cell and the porous anode layer[END_REF][START_REF] Lehnert | Modelling of gas transport phenomena in SOFC anodes[END_REF][START_REF] Xue | Dynamic modeling of single tubular SOFC combining heat/mass transfer and electrochemical reaction effects[END_REF][START_REF] Suwanwarangkul | Performance comparison of Fick's, dusty-gas and Stefan-Maxwell models to predict the concentration overpotential of a SOFC anode[END_REF][START_REF] Hussain | A general electrolyte-electrode-assembly model for the performance characteristics of planar anode-supported solid oxide fuel cells[END_REF][START_REF] Bhattacharyya | Dynamic modeling and validation studies of a tubular solid oxide fuel cell[END_REF] and mainly investigate on the effects of mass transfer on the performance of the cell. Other models take heat generation and transfer into account [START_REF] Kakaç | A review of numerical modeling of solid oxide fuel cells[END_REF][START_REF] Jin | Computational fluid dynamics analysis of solid oxide electrolysis cells with delaminations[END_REF][START_REF] Iora | Comparison of two IT DIR-SOFC models: Impact of variable thermodynamic, physical, and flow properties. Steady-state and dynamic analysis[END_REF][START_REF] Ni | Computational fluid dynamics modeling of a solid oxide electrolyzer cell for hydrogen production[END_REF][START_REF] Hwang | Heat/mass transfer in porous electrodes of fuel cells[END_REF][START_REF] Jeon | A comprehensive CFD model of anode-supported solid oxide fuel cells[END_REF][START_REF] Ho | Effects of heat sources on the performance of a planar solid oxide fuel cell[END_REF][START_REF] Yu | Permeability of fractal porous media by Monte Carlo simulations[END_REF][START_REF] Colclasure | Modeling and control of tubular solidoxide fuel cell systems. I: Physical models and linear model reduction[END_REF][START_REF] Murthy | Radiation heat transfer analysis of the monolith type solid oxide fuel cell[END_REF][START_REF] Wang | 3-D model of thermo-fluid and electrochemical for planar SOFC[END_REF][START_REF] Bove | Modeling solid oxide fuel cell operation: Approaches, techniques and results[END_REF][START_REF] Ni | Modeling of a solid oxide electrolysis cell for carbon dioxide electrolysis[END_REF][START_REF] Dutta | Numerical analysis of laminar flow and heat transfer in a high temperature electrolyzer[END_REF][START_REF] Hajimolana | Mathematical modeling of solid oxide fuel cells: A review[END_REF][START_REF] Authayanun | Maximizing the efficiency of a HT-PEMFC system integrated with glycerol reformer[END_REF] in order to evaluate the temperature in the cell and the influences of thermal effects. The planar design is the most studied , but other design are studied such as the tubular design [START_REF] Colclasure | Modeling and control of tubular solidoxide fuel cell systems. I: Physical models and linear model reduction[END_REF][START_REF] Dutta | Numerical analysis of laminar flow and heat transfer in a high temperature electrolyzer[END_REF] and the monolithic design [START_REF] Murthy | Radiation heat transfer analysis of the monolith type solid oxide fuel cell[END_REF]. The cell is generally divided cell into several subsystems in order to treat them with a specific set of equations: the gas channels, the porous electrodes, the interconnectors and the electrolyte.

A large majority of studies include heat convection and heat conduction in the PEN and in the channel, as well as heat generation caused by ohmic heat loss and by the electrochemical reaction. However, including radiation transfer inside heat transfer models is a subject of debate that has led to the production of numerous studies [START_REF] Pramuanjaroenkij | Mathematical analysis of planar solid oxide fuel cells[END_REF]. Indeed, radiation transfer inside solid oxide cells entailed a number of challenges in terms of modelling. First, the governing equations are integrodifferential and depend on time and many geometric variables, which lead to increased computational calculation time.

In addition, the electrolyte and the porous electrode are semi-transparent materials, whose optical properties are generally difficult to determine accurately [START_REF] Ho | Effects of heat sources on the performance of a planar solid oxide fuel cell[END_REF]. In addition, the effects of radiation appear to be very design-dependent: in 2007, Calise et al [START_REF] Calise | Simulation of a tubular solid oxide fuel cell through finite volume analysis: Effects of the radiative heat transfer and exergy analysis[END_REF] showed that radiation transfer contributes to a large majority of heat transfer inside tubular designs. For planar systems, Damm and Fedorov have conducted an analysis which concluded that radiation transfer could be safely neglected [START_REF] Damm | Radiation heat transfer in SOFC materials and components[END_REF] for thick and opaque electrodes. Other studies also concluded that heat radiation effects could be neglected in planar geometries [START_REF] Hajimolana | Mathematical modeling of solid oxide fuel cells: A review[END_REF]. On the contrary, other studies think radiation heat transfer must be considered in SOFC modelling [START_REF] Kakaç | A review of numerical modeling of solid oxide fuel cells[END_REF]. Finally, transient analysis has been performed by some groups [START_REF] Xie | Transient modeling of anode-supported solid oxide fuel cells[END_REF]. In 2005, Damm and Fedorov published a transient thermal model [START_REF] Damm | Reduced-order transient thermal modeling for SOFC heating and cooling[END_REF] showing the influence of temperature rise during startup/shutdown of the system. Their analysis also provided a characteristic time for startup, and a criterion to evaluate the maximum spatial and time thermal gradient in the cell.

Finally, Sunil Murthy and Andrei G. Fedorov [START_REF] Murthy | Radiation heat transfer analysis of the monolith type solid oxide fuel cell[END_REF] developed in 2003 a simplified radiation model that was tested on a monolith cell design in order reduce the computational cost by 10 times.

However, the majority of studies focus their work on the variations of temperature, concentration and current density on the directions of the flow, which is often parallel to the electrode. Although the PEN is considered as a volumetric electrode, only a few studies [START_REF] Xie | Transient modeling of anode-supported solid oxide fuel cells[END_REF][START_REF] Fardadi | Feedback control of solid oxide fuel cell spatial temperature variation[END_REF][START_REF] Martinez | Modeling and comparison to literature data of composite solid oxide fuel cell electrode-electrolyte interface conductivity[END_REF][START_REF] Autissier | CFD simulation tool for solid oxide fuel cells[END_REF] have given interest on the temperature variation on the direction perpendicular to the electrode, where strong temperature gradients can occur. Moreover, the reaction is generally supposed to occur in an active layer close to the interface with the electrolyte [START_REF] Andersson | Review on modeling development for multiscale chemical reactions coupled transport phenomena in solid oxide fuel cells[END_REF][START_REF] Ho | Effects of heat sources on the performance of a planar solid oxide fuel cell[END_REF]. However, our team [START_REF] Dumortier | Theoretical considerations on the modelling of transport in a three-phase electrode and application to a proton conducting solid oxide electrolysis cell[END_REF] previously showed that this assumption is valid under certain conditions : low ionic conductivity of ceramics, no limitation in the diffusion of species inside the electrode, metal grains radius not too large compared to ceramic grain radius. Our previous model also allowed predicting reaction profiles inside the electrodes. A reasonable question would be to know whether the reaction profile have an influence on heat release.

While models consider the porous medium as an equivalent medium with specific properties, other models deal separately with the energy transport inside the fluid and the solid phase. While those models are more descriptive and general, the equivalent medium models are supported by simplifying assumptions, such as the local thermal equilibrium (LTE) inside the PEN, thus disregarding the heat transfer by radiation. It has become relevant to investigate whether these phenomena should be considered inside the porous electrodes of a high temperature electrolyzer.

Another general issue of modelling is time calculation. The use of complex, computational time consuming models is not suitable for the performance prediction of larger systems such as stacks [START_REF] Iora | Comparison of two IT DIR-SOFC models: Impact of variable thermodynamic, physical, and flow properties. Steady-state and dynamic analysis[END_REF].

Therefore, it becomes necessary to find assumptions that simplify those models and reduce the resources required to perform calculations. It is however crucial to clearly define the scope in which those assumptions are valid. Another challenge for future solid oxide cells modelling development is the concept of inverse physics modelling, where the parameters of the process and microstructure of the electrode are calculated from the system requirements [START_REF] Ji | Monte-Carlo simulation and performance optimization for the cathode microstructure in a solid oxide fuel cell[END_REF]. This requires the building of analytical relationship between variables that do not require the help of numerical simulation to be solved.

The work presented in this paper focuses on planar solid oxide cells with electrodes made of cermets. The study is divided into three parts. In the first part, a continuity equation for thermal transfers inside the PEN in steady state is demonstrated. The assumptions for the establishment of this equation are checked with regards to the results of the literature. In the second part, the model is completed with a continuity equation set developed previously by our group for the description of transport phenomena inside a proton conducting high temperature electrolyzer [START_REF] Dumortier | Theoretical considerations on the modelling of transport in a three-phase electrode and application to a proton conducting solid oxide electrolysis cell[END_REF]. Few models have investigated on the thermal management of solid oxide electrolyzer [START_REF] Jin | Computational fluid dynamics analysis of solid oxide electrolysis cells with delaminations[END_REF][START_REF] Ni | Computational fluid dynamics modeling of a solid oxide electrolyzer cell for hydrogen production[END_REF][START_REF] Ni | Modeling of a solid oxide electrolysis cell for carbon dioxide electrolysis[END_REF], therefore the proposed model is tested throughout a parametric analysis inside the PEN of a proton-conducting cell in galvanostatic mode. The parameters tested are the inlet velocity of gases and inlet current density.

The results of this analysis allow modelling the thermal behaviour of the PEN as a thermal resistance network, which is presented in the third part. This analytical model reduces computing costs and offers the opportunity to do inverse physics modelling. One result of this inverse modelling is to state whether the process should be considered isothermal as regards to the operating conditions.

Modelling

The electrodes studied in this paper are based on cermets, which are considered to be among the most efficient candidates as electrode material for high-temperature steam electrolysis [START_REF] Deslouis | Characterization of porous texture of cermet electrode for steam electrolysis at intermediate temperature[END_REF]. A cermet is a mixed media made of three phases: an electron-conducting phase, which is usually a metal, an ionconducting phase, which is usually a ceramic, and a gas phase filling the voids between solids allowing the transport of reactants and reaction products (Fig. 2). The metallic and ceramic phases are motionless and the charge inside these solids is transported by migration. On the contrary, the gas phase does not carry any charge but transports chemical species by diffusion and convection. These transport phenomena inside the PEN are modelled mathematically by continuity equations obtained by applying the so-called volume averaging method [START_REF] Worster | Interfacial conditions between a pure fluid and a porous medium: implications for binary alloy solidification[END_REF] to a representative volume of electrode.

Continuity equations for charge and mass transfers have been derived by our group in a previous publication [START_REF] Dumortier | Theoretical considerations on the modelling of transport in a three-phase electrode and application to a proton conducting solid oxide electrolysis cell[END_REF]; as explained earlier, this paper focuses on the conservation of energy inside a PEN made of cermet.

Fig. 2 -The different phases of the cermet medium.

Heat transfers occur inside the fluid phase and the solid phase of the cermet. Even if the description of heat transfer is quite simple inside the solid phase, it is much more complex inside the fluid phase.

Hence, several assumptions have been made:

 Steady state.
 Because of low values of flow velocities and associated pressure gradient in the cermet and low dynamic viscosities of gas, compressibility effects and viscous dissipation are neglected.

 Thermal effects due to the mixing of gaseous species are neglected.

 Gravity effects are neglected.

Fluxes and sources

In the phase i (metal, ceramic or gas), heat is exchanged or transported by several mechanisms, expressed by heat fluxes (in W m -2 ):

 Convection flux inside the gas:

(1)

 Conduction flux inside the solid and the gas:

(2)

 Interfacial heat exchange with the solid j:

(3)

Inside charge conductors, heat is released by Joule's effect (in W m -3 ):

 Joule effect inside the solid i: (4)

In the fluid phase, the main heat sources come from the chemical reaction and thermal effects associated to the electronic transfers caused by the electrochemical reaction at the triple phase. We can assume that this source is produced inside the reactive volume and released inside the fluid phase.

Using the same demonstration process as featured by Coeuret [START_REF] Coeuret | Elements de génie électrochimique[END_REF], we can express the volumetric heat source p chim at the triple points as:

 Electrochemical heat source: [START_REF] Andersson | Review on modeling development for multiscale chemical reactions coupled transport phenomena in solid oxide fuel cells[END_REF] This volumetric source can be divided into two contributions:

 i r η accounts for the irreversibility caused by the charge transfer from one conductor to another [START_REF] Hwang | Heat/mass transfer in porous electrodes of fuel cells[END_REF]. It is in the form of a "Joule effect of the reaction" and is always positive. η is the overpotential defined by the generalized Butler-Volmer law.

 represents the electrical energy that is not utilized by the reaction. It accounts for thermodynamic losses that cannot be prevented during the process where ∆ r S is the molar entropy of reaction.

Code de champ modifié

Code de champ modifié

Finally, the heat radiation flux inside the PEN can be modelled as following: [START_REF] Soukup | Comparison of Wicke-Kallenbach and Graham's diffusion cells for obtaining transport characteristics of porous solids[END_REF] σ is the Stefan Boltzmann constant (5.67×10 -8 W m -2 K 4 ), n i is the refractive index of the medium.

Also, Damm and Fedorov [START_REF] Damm | Radiation heat transfer in SOFC materials and components[END_REF] have shown that the thermal radiation inside what? for thick PEN could be modelled by a heat conduction equation. They further demonstrated that radiation transfer between the PEN and other components of the cell could be neglected in planar systems because the thermal conductivity of the electrodes representing radiative heat transfer is much smaller than their overall thermal conductivity.

Local thermal equilibrium

The interfacial heat exchange flux is difficult to estimate as it involves the temperature of the phases. Hence, the objective of this model is to find simplifying assumptions that are valid in the case of solid oxide cells in order to consider the cermet as a single medium with respect to the energy conservation. One major issue that is currently discussed nowadays is the validity of the local thermal equilibrium assumption (LTE), which states that all the phases inside the PEN are locally at the same temperature. In our case, this would mean that:

T g = T io = T e = T ( 7
)
The applicability of this assumption for fuel cells has been quantitatively studied by Hwang and

Chen [START_REF] Hwang | Heat/mass transfer in porous electrodes of fuel cells[END_REF]. They assumed that the LTE is valid when the local difference in temperature between the fluid and the solid is small compared to the difference in temperature of the fluid between the inlet and the outlet of the electrode. Their analysis showed that the LTE can be assumed when with:

(8)

The mean velocity u mean of the fluid inside the electrode is the ratio of the volumetric flow rate produced to the surface of the electrode. This is equal to:

Code de champ modifié

Code de champ modifié [START_REF] Deseure | Theoretical optimisation of a SOFC composite cathode[END_REF] For the calculation of this criterion, gases can be considered as ideal:

( 1 0 ) (11)
The possibility to calculate the heat transfer coefficient h th(g-i) for flows in a porous media was investigated in 2001 by Kuwahara et al [START_REF] Kuwahara | A numerical study of interfacial convective heat transfer coefficient in two-energy equation model for convection in porous media[END_REF] and updated in 2010 by Pallares and Grau [START_REF] Pallares | A modification of a Nusselt number correlation for forced convection in porous media[END_REF].

Furthermore, Damn and Fedorov proposed in 2006 [START_REF] Damm | Local thermal non-equilibrium effects in porous electrodes of the hydrogen-fueled SOFC[END_REF] a conservative estimation of the heat transfer coefficient:

(12)
This is the most conservative estimation of h th(g-i) since it considers that heat transfer between the solid and the fluid is only made by conduction. Combining equations ( 8) to [START_REF] Jin | Computational fluid dynamics analysis of solid oxide electrolysis cells with delaminations[END_REF] gives a new expression of the LTE criterion:

(

) 13
For a sample calculation of H, we take values for solid oxide cells that could maximise this quantity:

J=10 4 A m -2 , | νox + νred |=1/2, λ g =0.05 W m -1 K -1 , ε g =0.5, d p =2×10 -5
m , e ce = 5×10 -5 m. This gives a value of H = 6×10 -5 , which confirm that the LTE can reasonably be considered valid inside the PEN of solid oxide cells. As a consequence, all the phases of the cermet are locally at the same temperature T.

Thus, the interfacial heat exchange between the solid and the fluid phase in Eq. ( 3) can safely be neglected as well as radiative heat transfer in Eq. ( 6).

We suppose that the mass that is produced inside the fluid phase is at the same temperature as the fluid phase. Thus, mass production only appears indirectly inside the energy balance. The conservation of energy in the fluid phase can be then written by the following integral equation:

(

) 14
Applying the same averaging method presented in previous publication [START_REF] Dumortier | Theoretical considerations on the modelling of transport in a three-phase electrode and application to a proton conducting solid oxide electrolysis cell[END_REF], the local equation of energy for the fluid phase can be written as follows:

(

For the solid phase i:

(

Note that without additional approximation, this equation can also be written as follows :

For numerical calculations we can use a more suitable expression by stating the last equation with equivalent variables for the porous medium:

(

This equation shows that the cermet can be modelled as a single phase equivalent medium for the thermal phenomena, with the equivalent properties:

( 2 0 )

Code de champ modifié

( 2 1 ) (22)
Equation [START_REF] Jeon | A comprehensive CFD model of anode-supported solid oxide fuel cells[END_REF] shows that the equivalent thermal resistance of the medium is comparable to the parallel connection of the resistance of each phase. Thermal energy conservation in the porous electrodes is now established and can be modelled into one continuity equation that features equivalent properties. This model will be now tested numerically on a proton conducting solid oxide electrolysis cell.

Numerical model

The properties of gases (hydrogen, oxygen and water vapour) are taken from the literature [START_REF] Kestin | Thermophysical Properties of Fluid H2O. JPhysChem[END_REF][START_REF] Kunz | The GERG-2004 Wide-Range Reference Equation of State for Natural Gases and Other Mixtures[END_REF][START_REF] Laesecke | Transport Properties of Fluid Oxygen[END_REF][START_REF] Leachman | Fundamental Equations of State for Parahydrogen, Normal Hydrogen, and Orthohydrogen[END_REF][START_REF] Schmidt | A New Form of the Equation of State for Pure Substances and its Application to Oxygen[END_REF]. For a given pressure, each quantity can be written as a polynomial function of temperature, whose coefficients are given in Table 1.

(

) 23
Table 1 -Coefficient for polynomial functions to calculate the properties of gases at 10 bars between 773 K and 973 K.

Quantity Coefficients Steam Hydrogen Oxygen a h - - - b h 2.22×10 4 1.48×10 4 1.07×10 3 Specific enthalpy (J kg -1 ) c h 1.76×10 7 -6.09×10 5 -8.84×10 4 a k - - - b k 1.31×10 -4 4.96×10 -4 6.97×10 -5 Thermal conductivity (W m -1 K -1 ) c k -3.41×10 -2 2.86×10 -2 1.02×10 -2 a s - - - b s 2.56×10 -3 1.70×10 -2 1.23×10 -3 Specific entropy (J kg -1 K -1 ) c s 5.73 44.6 5.79 a c 3.74×10 -4 1.51×10 -3 -1.71×10 -4 b c -5.49 ×10 -2 -1.31 0.482 Specific heat capacity (J kg -1 K -1 ) c c 1.99 ×10 3 1.48×10 4 7.80×10 -2 a μ - - - Dynamic viscosity (Pa s) b μ 4.00×10 -8 1.44×10 -8 3.45×10 -8

Code de champ modifié

Supprimé :

à c μ -2.34×10 -6 5.89×10 -6 1.48×10 -5
The thermal model is then coupled to the mass transfer model published earlier by our group [START_REF] Dumortier | Theoretical considerations on the modelling of transport in a three-phase electrode and application to a proton conducting solid oxide electrolysis cell[END_REF].

Calculations are run on COMSOL Multiphysics 4.2. considering the same reactor dimensions (Table 2) and boundary conditions, parameters and reference values, except for those presented in Table 3.

The conductivity of protons inside the ceramic was measured by the work team by impedance spectroscopy. Measurements were carried out between 748K and 873K -10 bars.

The values considered for the parametric study encloses realistic operating data for solid oxide cells.

Some values are out of the operating range for those cells : for example a value for inlet current density of J = 1.5×10 4 A m -2 is considered whereas solid oxide cells usually operated under 1×10 4 A m -2 . Also, an operating temperature of T = 973 K is studied, which is too high for the conduction of protons inside the electrolyte [START_REF] Slodczyk | Structural modifications induced by free protons in proton conducting perovskite zirconate membrane[END_REF]. These values are tested for theoretical purpose to evaluate the influence of the parameter or to predict the behaviour of the system in the case of an accidental malfunctioning.

Unless another value of those parameters is specified, the default value that is taken for any calculation is the one presented in Table 3. The indexes of the boundary conditions are given in Fig. 3 and the additional boundary conditions for the thermal problem are given in Table 4.

Parametric analysis

An analysis is carried out on the parameters that are suspected to have a strong influence on power and temperature, i.e. the inlet current density, which generates ohmic heat losses, and the inlet velocity.

The influence of the operating temperature is also tested in order to predict the behaviour of the cell at different working temperatures that can be reached. Inlet current is obviously important as it determines the mass consumed or produced in the process. The conduction of charge also generates ohmic heat losses that are proportional to the square of the value of inlet current density. Inlet velocity is a very influential parameter on the global energetic and economic performance of solid oxide cells.

Low velocities gas flows allow larger reactant utilization, but can cause local reactant depletion and cold spots that exacerbate temperature non-uniformities. On the contrary, higher gas flows tend to increase the uniformity of the reaction rates across the active area but decrease reactant utilization.

Therefore, management of gas flows is critical to the stable operation of the cell [START_REF] Wang | 3-D model of thermo-fluid and electrochemical for planar SOFC[END_REF]. The Reynolds number in this study ranges from 0.9 to 44, therefore the hydraulic conditions remains globally laminar.

It should be mentioned that the numerical results presented in this study are dependent on the design and dimensions of the device. Those dimensions are not clearly defined as there is currently no definitive design for those cells [START_REF] Andersson | Review on modeling development for multiscale chemical reactions coupled transport phenomena in solid oxide fuel cells[END_REF]. Therefore the results obtained with the parametric study are not meant to be realistic; rather they are used to show tendencies and relationships between the different phenomena of the system. From these observations that are presented in this section, a series of assumptions are built for the modelling of the PEN into a network of thermal resistances.

Code de champ modifié

Code de champ modifié

Influence of operating temperature

The influence of the operating temperature on the distribution of the transfer current density is tested in order to determine whether this parameter should be monitored. The temperatures range from 773 K to 873 K, which corresponds to a reasonable range for which the perovskite can conduct protons. An additional value of 973 K is tested for the reasons presented in section 3.

From Figure 4 we notice that in this range of temperature variation, the oxygen molar fraction profiles at the anode are not affected (Fig. 4). Indeed, the ratio between the diffusion coefficients of oxygen and steam, which are transported in opposite direction mainly by Knudsen diffusion, remains the same with temperature:

( 2 4 )

The influence of the operating temperature on the transfer current density across the electrode is shown in Fig. 5. This influence is minor compared to other influential parameters that have been tested in a previous study, such as the grain radius of the porous media and the ratio between the conductivity of the metal and the ceramic σ io /σ e [START_REF] Dumortier | Theoretical considerations on the modelling of transport in a three-phase electrode and application to a proton conducting solid oxide electrolysis cell[END_REF]. Operating temperature appears in the expressions of many quantities of the system. The little variations shown in Fig. 5 are suspected to be caused by the effect of the temperature in the Butler Volmer equation and by the slight change in σ io /σ e , which ranges from 0.9×10 -5 at 773K to 2×10 -5 at 973K.

Effect of inlet current on heat release and temperature profile

An important question about heat sources is to know where and how they take place. In order to underline this problem the effect of the inlet current on heat release and temperature profile is studied.

The results of the calculations are presented in Tables 5 and6. Firstly, we observe that a very large majority of this energy release is caused by Joule's effect inside the ceramics (98%). This is caused by the low ionic conductivity of this phase compared to the electronic conductivity of metal.

Secondly, the energy is also mainly released in the electrolyte (89%) despite its low thickness compared to the electrode Il y aquelque chose que je ne comprends pas (effet de copier-coller?) difference entre ceramics et electrolyte(Table 1). Indeed, the distance covered by protons, which is directly linked with heat power generation, is thinner in the electrodes than in the electrolyte because the electrochemical reaction occurs close to the interface with the electrolyte. Actually, protons tend to minimize their path through the ceramic when the ionic conductivity of the ceramic is low: charge carriers choose the easiest conduction path il me semble qu'il y a un principe de physique général derrière cette affirmation, " principle of minimum rate of energy dissipation in dynamics equilibrium"

Voir papiers attachés . However, the most important observation about these results is that the location and the nature (ohmic or electrochemical) of the heat sources remain roughly independent of the value of the current density (±4% disparity).

This means that for a given problem, the values of the heat sources in the cell can be predicted with only few parameters: the dimensions of the electrolyte on voit mal ce que c'est!, the value of the ionic conductivity, the inlet current density and the proportion of the total power released in the electrodes.

These observations should however be confirmed in special cases, for example if the electronic conductivity of the metal is close to the ionic conductivity of the ceramics and in the case where diffusion phenomena become limiting.

Note also that the electrochemical heat source is negligible when inlet current density is larger than 1000 A m -2 . This result is valid for the electrolysis of steam, but should be checked for technologies using methane or carbon dioxide, for which enthalpies and entropies of reaction are different [START_REF] Ni | Computational fluid dynamics modeling of a solid oxide electrolyzer cell for hydrogen production[END_REF].

Effect of velocity on heat release and temperature profile

Different inlet velocities inside the channels are tested from 0.1 m s -1 to 10 m s -1 . The flow remains laminar for the velocities studied. Like current density, the value of velocity has little influence on the location and the nature of the heat sources inside the cells shown in Table 5, mainly because the phenomena inside the electrode are driven by the electrode structure and the value of the inlet current density [START_REF] Dumortier | Theoretical considerations on the modelling of transport in a three-phase electrode and application to a proton conducting solid oxide electrolysis cell[END_REF]. However, this has a huge influence on the maximum temperature reached inside the electrolyzer. Indeed, the fluid extracts the heat from the electrodes and conveys it to the outlet. As gases do not convey heat as efficiently as liquids, the management of the mass rate, the pressure and the composition of the gas flows are important parameters for monitoring the maximum temperature inside the cell. Fig 6 shows that the maximum temperature reached inside the electrolyte is dependent on the inlet velocity. The value of the maximum temperature is not meant to be representative for real electrolyzer, rather it will be used for building the resistance model that is presented. The maximum temperature reached in the cell is 913 K, i.e. an elevation of 40K, for the maximum inlet current density J=1.5×10 4 A m -2 and the minimum inlet velocity U=0.1 m s -1 .

It can be noticed that temperature remains relatively constant in the PEN and varies strongly in the channels. Indeed, the thermal resistance of the PEN (approximately 5×10 -3 K W -1 ) is low compared to the high thermal resistance of the gas flow (approximately 2 K W -1 ), therefore the temperature profile is more homogeneous in the PEN. However, literature review shows that the temperature of higher aspect ratio PEN can vary significantly along the direction of the flow [START_REF] Hwang | Heat/mass transfer in porous electrodes of fuel cells[END_REF][START_REF] Jeon | A comprehensive CFD model of anode-supported solid oxide fuel cells[END_REF][START_REF] Ho | Effects of heat sources on the performance of a planar solid oxide fuel cell[END_REF][START_REF] Wang | 3-D model of thermo-fluid and electrochemical for planar SOFC[END_REF]. Thus, the previous observation seems valid only for a PEN with a low aspect ratio. The temperature profile in the cathode channel slightly differs from the anode: one can notice that the thermal boundary layer is thinner in the anodic channel than in the cathodic channel. This is because the Prandtl number (Eq. 25), which compares the relative thickness of the momentum and the thermal boundary layer, is not the same for the gases at the anode and at the cathode.

Pr = μ g Cp (g) /λ g ( 2 5 )
As a result, the interaction between the fluid and the heat flux is higher for the gas mixture at the anode (Pr=1.24) than for the hydrogen at the cathode (Pr=0.6).

Thermal resistances modelling

The objective of this model is to represent the thermal behaviour of the cell as a network of thermal resistances and heat sources. The behaviour of a component of the cell is given by the following equation:

∆T = R th Φ th ( 2 6 )
∆T is the difference between the mean temperature of the two boundaries of the component (in K), Φ th is the thermal flux that crosses the element (in W) and R th is the thermal resistance of the component (in W K -1 ). This model relies on different assumptions:

 Heat is transported through the component of the PEN and the gas channels only.

 Heat sources are uniform along the directions of the gas flows.

 The heat sources inside the electrodes are generated at the electrode/electrolyte interface.

 The properties of the electrodes can only vary in their depth.

 Heat sources are not influenced by the velocity of gases inside the channel  Steady state

Thermal resistances

This thermal resistance R th can be related to the parameters of the medium. In the electrodes and in the electrolyte, heat transfer occurs through conduction and the resistance is given by

R th = e/λS ( 2 7 )
e is the width of the medium, λ its thermal conductivity and S the area crossed by the flux (or the length in 2D models).

The thermal resistance in the gas channels depends on the characteristics of the diffusion and thermal layers at the interface with the electrode and is usually given by a correlation. So far, no correlation has been found to match the results of this study, however data reduction from the COMSOL calculations gives the expression of a thermal resistance for the channels in this model: R th = 1/(0.24 Re 0.53 Pr 0.23 ) ( 2 8 )

Heat sources

We have shown in the parametric study that heat is mostly generated by ohmic losses inside the electrodes. We have also shown that the heat source in the electrodes is a proportion of the heat source in the electrolyte. Therefore all the heat sources depend only on the value of the inlet current, electrode external surface and thickness, and ionic conductivity.

For the electrolyte:

(

The expression for total power source inside the electrodes is given by:

( 3 0 )

According to the results shown in Table 4, we calculate K =0.1. From the observations made in section 4, we have shown that the value of K is independent of the inlet current density and inlet velocity.

Supprimé : from

However it probably depends on the electrode microstructure, on the electrical conductivities and other parameters of the electrode.

Equivalent circuit

Since the sources are not dependent on the phenomenon of the process and that the calculations are performed in steady state, the thermal transfer inside the cell can be modelled in the form of a thermal resistances network (Fig. 7). The resistance of the different components are in series while the heat sources are symbolised as current generators in the equivalent circuit. R ace represents the thermal resistance of the anode, R e represents the thermal resistance of the electrolyte and R acc represents the thermal resistance of the cathode. The expression of the thermal conductivities of the electrodes is given by Eq. ( 22). The thermal resistances of the anodic and cathodic gas channels anodic are respectively called R ach and R cch . T 0 is the temperature of inlet gases while T ach , T ace , T e , T cce and T cch are the maximum temperature respectively reached inside the anodic gas channel, the anode, the electrolyte, the cathode and the cathodic gas channel. R ∞ is an infinite resistance needed for the resolution of the system and has no physical meaning.

The objective of this network system is to predict the maximum temperatures inside each component.

By applying the Kirchoff's circuit law on each node of the network, we obtain the representative matrix of the system.

For each node k, Kirchoff's circuit law is given by:

(T k-1 -T k )/R k-1 + (T k+1 -T k )/R k + Φ k = 0 ( 3 1 )
The matricial representation of the system is given by:

[R]×[T]=[P] ( 3 2
Solving this network allows the prediction of the maximum temperature reached inside each component of the PEN. The temperatures in the channels are then obtained with the potential divider law. The expression of maximum temperature is given by Eq. ( 36) whose coefficients are given in Table 6.

(36)

Inverse physics and isothermal assumption checking

The resistance network model analytically relates the maximum temperature attained in the components to the parameter of the process. This model also allows determining for which values of those parameters a given temperature increase is attained.

For a sample calculation, let us consider that the heat sources inside the electrodes are negligible. Let us also consider that the properties of the two electrodes are equal and that all the elements of the PEN have the same external surface. The maximum temperature is therefore attained in the middle of the electrolyte and is given by: (

Using the expressions of the resistance and heat sources given previously, we obtain:

(

This equation shows how the parameters of the process influence the value of the T e . Given a set of parameters, it is therefore possible to determine whether the process is isothermal or not before running an experiment or a calculation. Changing the form of this equation gives the value of a parameter for a given T e . Thus, resistance modelling makes inverse physics modelling is possible, with or without the assumptions given for this sample calculation.

Conclusion

This paper has presented different models for the simulation of thermal transfer inside the electrode assembly of a high temperature electrolyzer. The local thermal equilibrium assumption in the cermet was demonstrated through the construction of a criterion, and the continuity equation for energy was built via a volume averaging method. This model was implemented in COMSOL Multiphysics along with a set of equations demonstrated in previous work for the simulation of a planar proton conducting solid oxide electrolysis cell in galvanostatic mode. A parametric analysis was carried out by studying the influence of inlet current density and inlet velocity of gases in the channels on temperature and heat released in the cell. This analysis showed that nearly all the heat power source is released in the ceramic phase by Joule's effect and that the temperature varies mostly inside the channels. Moreover, the repartition of heat power sources inside the electrode and the electrolyte is mostly independent from inlet current density or from the inlet velocity of gases. This means that for a given problem, the heat release in the cell can be predicted with only few parameters: the dimensions of the electrolyte, the value of the ionic conductivity, the inlet current density and the proportion of power released in the electrodes. Following that observation, it was demonstrated that it is possible to model the thermal behaviour of the cell as a network of heat sources and thermal resistance in order to analytically evaluate the temperature inside each component of the cell. A comparison with the results given by the parametric analysis showed that this model predict those temperatures accurately. This modelling strategy strongly reduces computation time and allows reverse physical analysis in order to determine for which values of those parameters a given temperature increase is attained.

λ

  Abbreviations

6 )

 6 By coupling the previous equations together, and by summing the energy equation of each phase, the local conservation equation for energy inside the control volume is written as:

Fig. 4 -Fig. 5 -

 45 Fig. 4 -Oxygen fraction in the anode channel for different operating temperature.

Fig. 6 -

 6 Fig. 6 -Temperature elevation in the cell (in K) for different inlet velocities.

Fig. 7 -

 7 Fig. 7 -Thermal network modelling of the electrolyzer.

Table 2 -

 2 Reactor dimensions

	Parameter	Value
	Anodic channel width	4×10 -4 m
	Anode width	1.5×10 -4 m
	Electrolyte width	1×10 -4 m
	Cathode width	1.5×10 -4 m
	Cathodic channel width	4×10 -4 m
	Anode length	10 -3 m
	Electrolyte length	10 -3 m
	Cathode length	10 -3 m
	Anodic channel length	2×10 -3 m
	Cathodic channel length	2×10 -3 m

Table 3 -

 3 Values of variable parameters and inlet parameters for the reference calculation

	Parameter	Reference value	Unit	Other values tested
	Inlet current density Inlet velocity	10000 0.1	A m -2 m s -1	1, 10 3 , 2×10 3 , 5×10 3 , 10 4 , 1.5×10 4 0.2, 0.5, 1, 1.5, 2
	Operating temperature Operating pressure	873 10 6	K Pa	773, 973 -
	Particle radius (metal)	10 -6	m	-
	Particle radius (ceramic)	5×10 -7	m	-
	Initial oxygen molar fraction	0,01
	Ionic conductivity (ceramic)	1.2302×e -1827.6/T	S m -1	-
	Ceramic thermal conductivity	1 [152]	W m -1 K -1 -
	Metal thermal conductivity	1000	W m -1 K -1 -
	Proton transport number	1	-

Fig. 3 -Reactor geometry and boundary conditions.

Table 4 -

 4 Boundary conditions

	Boundary	Heat
	1	Zero flux
	2	Fixed temperature (T 0 =873 K)
	3	Outlet , no temperature gradient
	4	Zero flux
	5	Heat flux continuity
	6,8,10,11,13,14 Zero flux
	7, 15	Zero flux

Table 5

 5 Repartition of the heat release inside the PEN

	Parameter	Anode Electrolyte Cathode Total
	Joule's effect (ceramic)	8%	82%	8%	98%
	Joule's effect (metal)	1%	-	1%	2%
	Electrochemical heat source	0.1%	-	0.1%	0.1%
	Total	9%	82%	9%	100%

Table 6

 6 Total heat released in the PEN

	Inlet current density (A m -2 )	Total power (W)
	1	9.55×10 -8
	1000	9.54×10 -2
	10000	8.91

Table 7

 7 Coefficients for the calculation of temperature T ach R ach (R e +R cch +R cce ) R ach (R e /2+R cch +R cce ) R ach (R cch +R cce ) T ace (R ach +R ace )(R e +R cch +R cce ) (R ach +R ace )(R e /2+R cch +R cce ) (R ach +R ace )(R cch +R cce ) T e (R e /2+R cch +R cce )(R ach +R ace ) (R ach +R ace +R e /2)(R e /2+R cch +R cce ) (R ach +R ace +R e /2)(R cch +R cce ) T cce (R cch +R cce )(R ach +R ace ) (R cch +R cce )(R e /2+R ach +R ace ) (R cch +R cce )(R e +R ach +R ace ) T cch R cch (R ach +R ace ) R cch (R e /2+R ach +R ace ) R cch (R e +R ach +R ace )

	a a	a e	a c

Acknowledgements

The authors acknowledge the French ANR (Agence Nationale pour la Recherche, France) for the financial support (Project ANR-09HPAC-007-01). The authors thank Kamal Rahmouni for the measurements of the ionic conductivity.

Supprimé : o

Code de champ modifié Code de champ modifié Code de champ modifié Code de champ modifié Code de champ modifié Code de champ modifié Code de champ modifié Code de champ modifié Code de champ modifié Code de champ modifié Code de champ modifié Code de champ modifié Code de champ modifié Code de champ modifié Code de champ modifié Code de champ modifié Code de champ modifié Code de champ modifié Code de champ modifié Code de champ modifié Code de champ modifié Code de champ modifié Code de champ modifié Code de champ modifié Code de champ modifié Code de champ modifié Code de champ modifié Code de champ modifié Code de champ modifié Code de champ modifié Code de champ modifié Code de champ modifié Code de champ modifié Code de champ modifié Code de champ modifié Code de champ modifié Code de champ modifié 2.3 Energy conservation Code de champ modifié Code de champ modifié Supprimé : y Code de champ modifié Supprimé : variations Supprimé : side Code de champ modifié Supprimé : from Code de champ modifié Code de champ modifié Code de champ modifié Supprimé : they

Comparison with the results of the parametric analysis

The model has been tested for every case that has been performed in the parametric analysis with COMSOL. The results given in Table 8 and Table 9 show that the temperatures given by the network model are very close to the temperatures given by the parametric analysis, with maximum difference of 1K between the two models. The resistance network model is therefore valid to estimate the maximum temperatures given by the calculations in COMSOL. It must be noted however that the thermal resistance modelling is valid if the variation of temperature inside a component does not cause a large variation in the value of its thermal resistance. x i Molar fraction for species i (dimensionless)