A multiscale physical model of a polymer electrolyte membrane water electrolyzer

Electrochimica Acta 110 (2013) 363–374

Contents lists available at ScienceDirect

Electrochimica Acta journal homepage: www.elsevier.com/locate/electacta

A multiscale physical model of a polymer electrolyte membrane water electrolyzer Luiz Fernando L. Oliveira a,2 , Christian Jallut b , Alejandro A. Franco a,∗,1 a b

Alternative and Atomic Energies Commission of France CEA/LITEN, 17 Rue des Martyrs, 38054 Grenoble cedex 9, France Université Claude Bernard Lyon 1, CNRS, UMR 5007, LAGEP, F-69622 Villeurbanne, France

a r t i c l e

i n f o

Article history: Received 17 December 2012 Received in revised form 29 July 2013 Accepted 29 July 2013 Available online 13 August 2013 Keywords: Water electrolysis Hydrogen production PEMWEs Transition metal oxides Multiscale modeling

a b s t r a c t In this paper we report a multiscale physical and transient model describing the operation of a polymer electrolyte membrane water electrolyzer single cell. This model includes a detailed description of the elementary electrode kinetics, a description of the behavior of the nanoscale catalyst–electrolyte interface, and a microstructural description of the transport of chemical species and charges at the microscale along the whole membrane electrodes assembly (MEA). We present an impact study of different catalyst materials on the performance of the PEMWEs and a sensitivity study to the operation conditions, both evaluated from numerical simulations and with results discussed in comparison with available experimental data. © 2013 Elsevier Ltd. All rights reserved.

1. Introduction Electrochemical dissociation of water has been known for many years. Rudolph and van Trootswijk are credited to be the first researchers who split water into oxygen and hydrogen, back in 1788. Later, in 1800, Nicholas and Carlisle performed electrolysis using the Voltaic cell which was invented in the same year [1]. Nowadays, different approaches have been proposed in order to have an industrial production of hydrogen, using water electrolysis. These differences depend upon the type of electrolyte and electrodes and therefore on their performance. One of the technologies used is the so-called PEM water electrolyzers (PEMWEs). PEMWEs can be considered as a reverse device to the well-known polymer electrolyte membrane fuel cells (PEMFCs). PEMWEs use a perfluorosulfonic acid (PFSA) proton conducting polymer as electrolyte. Water is oxidized in the anode, through the so-called oxygen evolution reaction (OER), with the production of molecular oxygen, O2 , and protons, H+ . Protons flow to the cathode side through the

∗ Corresponding author. Present address: a) Laboratoire de Réactivité et de Chimie des Solides (LRCS) Université de Picardie Jules Verne & CNRS, UMR 7314 Amiens, France. b) Réseau sur le Stockage Electrochimique de l’Energie (RS2E), FR CNRS 3459, France. Tel.: +33 3 22 82 75 86. E-mail address: [email protected] (A.A. Franco). 1 ISE member. 2 Present address: Laboratoire de Chimie et Physique Quantiques LCPQ/IRSAMC, Université de Toulouse (UPS) and CNRS, 118 route deNarbonne, F-31062 Toulouse, France (present affiliation of Luiz Fernando Oliveira Lopes). 0013-4686/$ – see front matter © 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.electacta.2013.07.214

membrane, and O2 is either released to the atmosphere or collected. The electrons produced through the reaction circulate through an external circuit. Finally, in the cathode side H+ are reduced to H2 , through the so-called hydrogen evolution reaction (HER). Such devices, offer several advantages, including ecological cleanness, better efficiency, from both current density and energy, compactness and low temperature operation [2–4]. Although PEM technology was first used in the 1960s [5], PEMWEs started to be conceived in the 1980s [6], and it received more attention from the scientific community only from the middle of the 1990s and the beginning of the 2000s. Different studies have discussed the promising features of PEMWEs, for the production of pure hydrogen from renewable energy sources [2–4,7]. Nowadays PEMWEs are already commercially available. Some of the PEMWEs manufactures are: Giner [8] and Proton OnSite [9] from USA; CETH [10] from France; and Hydrogenics [11] from Canada. As far as the price is concerned, it can reach several hundreds of euros. The price depends mainly on the required efficiency. Nevertheless, to turn this technology widely available, there are still technical and scientific issues to be overcome. One of the main drawbacks is the price of the material used in the electrodes. Concerning the catalyst material, platinum nanoparticles supported on percolated carbon nanoparticles are used for the HER in the cathode electrode. Whereas in the anode side, metal oxides, like IrO2 and RuO2 are currently used, with the double role of OER catalyst and electronic conductor. Other materials, mainly alloys, have also been

364

L.F.L. Oliveira et al. / Electrochimica Acta 110 (2013) 363–374

proposed for both electrodes [44,46]. The ionomer, which ensures the proton conduction within and between the electrodes, is usually made of PFSA proton conducting polymers, and it also has an elevated cost. The ionomer, developed by Dupont, is the most used one in PEMWEs, since its proton conductivity properties are proven extremely good [12]. In addition to the issues related with cost, there are many phenomena behind the operation principles still not well understood. One example is the mechanism of the electrochemical reactions taking place in a PEMWEs. The possible interactions between the ionomer and the catalyst materials are also not very well described in the literature. It is known, however, that these types of interactions cannot be neglected for the case of PEMFCs [13,14]. Moreover, PEMWEs present substantial technical challenges related to their efficiency and lifetime. For example, the catalyst dissolution, that leads to an evolution of the catalyst layer microstructure properties, remains a problem [15]. A further problem is the degradation of the membrane, a well known issue for the case of PEMFCs but poorly described in the literature for PEMWEs [16]. To enhance the performance and durability of the PEMWEs, a deep understanding of the physicochemical processes, related to the nano and microstructural properties of the electrodes, is crucial. With this aim, the complex mechanisms occurring at the multiple scales, under PEMWEs operation conditions have been a subject for numerous experimental works. Clearly there are phenomena that cannot be completely analyzed from experiments, mainly due to the strong coupling between those phenomena which take place simultaneously under the PEMWEs operation. In addition, experimental efforts concentrate in dealing only with one of the above cited challenges. In contrast to this, physical models allow a detailed analysis of processes according to various temporal and spatial scales. In each scale different approaches should be used, and sub-models can be built up. These sub-models can be connected giving rise to the socalled multiscale or multiphysics models. Concerning the modeling of PEMWEs, differently of PEMFCs, not many physical models have been reported. In the majority of the PEMWEs cell modeling studies, empirically fitted equations and/or Butler–Volmer equations are used for the electrochemical reactions description. Such equations do not consider atomistic processes, since the reactions are described using effective global steps. Thus, this type of approach does not allow connection of the materials chemical and structural properties with the experimental observables. Models that use the neural network approach [17], or equivalent electrical circuits [18] have been developed. Such investigations concern mainly macroscopical description of the cell operation providing insights on the temperature profiles across the cell or on water flow management aspects. Busquet et al. [19] proposed an empirical model with the aim to simulate the i–V curves from its typical behavior. The model can be used for fuel cells, electrolyzers and regenerative fuel cell. Choi et al. [20] presented an analytical model based on material balance equations and an equivalent circuit to the electrical analysis of the electrodes. Onda et al. [21] built a model similar to the one of [20], using empirical equations, they concern mainly in the ideal open-circuit voltage. Few modeling studies show sensitivity analysis of the i–V curves to important operation parameters, for instance, results have been shown only for temperature and/or pressure [18,19]. Although some preliminary multiphysical model of an electrolyzer has been reported [22], published results showed only the dependency of the cell operation on the temperature, and the study

of the impact of other important operating conditions have not been presented. Only few works take into consideration important physicochemical phenomena that occur under real conditions. One example is the work done by Marangio et al. [23], in which they have combined a theoretical model with experimental measurements. Based on material balance equations and electrical analysis of the electrodes and plates, the open-circuit voltage was calculated. They presented i–V curves and sensitivity studies to different parameters, as for example the concentrations and diffusivity of H+ in the membrane. Finally a comparison with the experimental results were made. A second example is the work by Grigoriev et al. [24]. Using a multiphase flow modeling, different macroscopical phenomena were analyzed which include cross-permeation of the gas through the membrane, formation of gas bubbles, flow rate of water and heat production. To the best of our knowledge, any physicochemical model allowing the prediction of the performance and durability of a PEMWEs as a function of its operation conditions and the materials chemical and structural properties – catalyst loading and microporosities at multiple scales – has been reported in the literature. In one previous paper [25] we have proposed a first model within this sense to address the simulation of IrO2 -based PEMWEs anodes. In this paper we propose an extension of this model to the whole MEA. First, we present the model and the coupled scales accounted for. Then, we discuss simulation results regarding different parameter sensitivity studies, and finally we conclude and highlight further work necessary to improve our model. 2. Methodology 2.1. Multiscale model: general presentation The multiscale model presented here was first developed for PEMFCs by Franco et al. [26–29]. This model is based on different sub-models, each of them dealing with a specific mechanism. Several sub-models can be included in order to simulate the different processes. This modeling approach also takes into account the feedback between the transient performance and the intrinsic material aging process thus allowing the prediction of the cell durability. For the case of PEMFCs, this approach has provided new insights on the interaction between different aging phenomena. It is also possible to analyze the cell response to operational parameters such as temperature, pressure and relative humidity. Many other parameters such as the loadings of the catalyst material, the carbon phase and the ionomer can also be easily studied. It is also possible to study the behavior of time-dependent variables at the different scales, as the time evolution of the electrocatalytic activity. Within this context, the model has provided information on the competition of aging phenomena [30,47,48]. This model can be extended not only for PEMWEs, as it has been done in the present work, but for other electrochemical devices like Li-ion batteries and supercapacitors [31,32]. Similar parameters as those mentioned in the previous paragraph can be studied as well. In the following sections the sub-models for PEMWEs anode, cathode and membrane are detailed. In what concerns the multiscale model we consider the following assumptions: • the ionomer is fully saturated by water, so that water transport is not described; • the cell operates under isothermal and isobaric conditions; • the model neglects mechanical constraints (e.g. at the membrane/electrode interface);

L.F.L. Oliveira et al. / Electrochimica Acta 110 (2013) 363–374

365

• the membrane is impermeable to O2 and H2 ; • the complex geometry of the electrodes is simplified by considering a cylindrical structure.

namely the difference of potential in the compact layer, is defined, for both, anode and cathode, as

Concerning the last assumption some clarifications are needed. The real geometry of the electrodes is complex due to the fact that the three phases, the liquid water, the ionomer and the solid catalyst, are intimately brought into contact forming a three dimensional structure. Each of the three phases under consideration has to be continuous so that the transport of all the species involved in the electrochemical process is possible. Such a complex geometry is simplified by considering a cylindrical structure, which satisfies the continuity of the three phases. In Sections 2.2 and 2.3 models dealing with the phenomena occurring, respectively, at the nanoscale and microscale are described. In Fig. 1 it is shown a flowchart showing the interconnection of the main equations used to model the phenomena which are discussed in the following sections. Further details and the derivation of the equations presented therein are given in the several previous works by Franco et al. [26–28,33,34].

where N (L, y, t) (N stands for nanoscale) is the electrostatic potential in the electrolyte at the boundary between the diffuse and compact layers, and (y, t) is the catalyst electrostatic potential. The potential (y, t) is used in the calculation of the reaction rates, described in the next section. Furthermore, given the potential of the catalyst phase, (y, t), of each electrode, one can calculate the potential of the cell (Section 2.5). The approach used to calculate the potential in the electrolyte phase, N , is detailed in Section 2.2.3. Hereafter it is presented the equations used to calculate (y, t). It has already been postulated by Franco et al. [27,28] that, in addition to Eq. (1), (y, t) can also be calculated by the following equation:

2.2. Nanoscale model The elementary kinetics part of the model accounts for an explicit description of the electrochemical double layer formation in the vicinity of the catalyst under conditions far from the equilibrium. This is one of the key parts of our model and in order to account for the physicochemical phenomena at the nanoscale, ab initio data can be used. The present section is dedicated to explain our approach in modeling the non equilibrium electrochemical double layer of both the cathode and the anode side. The compact layer is modeled by a 0-D approach, which is presented in Sections 2.2.1 and in Section 2.2.2. Section 2.2.3 is dedicated to the description of the 1-D approach used to study the transport of protons and chemical species through the diffuse layer. In this part of the model, the sulfonate groups related with the Nafion, are taken into account. However, since they are considered as fixed negatively charges, only their concentration is considered. An illustration of our anode nanoscale model can be seen in our previous publication [25]. In this case the species considered in the compact layer are H2 O, O2 , OH, OOH and H+ , that is, the intermediate species produced and consumed in the OER. In the diffuse layer we consider mainly the transport of H+ and O2 . Similar considerations hold for the cathode nanoscale, Fig. 2. The main differences are the directions of the charge and mass transport and obviously the chemical species involved. Concerning the transport of chemical species through the diffuse layer, here we consider H+ coming from the anode through the membrane, and the H2 produced. As a result, in the compact layer, the species treated are H+ , H2 , which are related to the HER, and the H2 O. The nanoscale model is connected with the microscale model (Section 2.3) through the expressions of its boundary conditions and source terms (Eqs. (45) and (47)). 2.2.1. Compact layer model This part of the nanoscale model describes the interaction between the chemical species and the catalyst surface. Adsorption and desorption of molecules and reaction intermediates as well as charge transfer processes take place at this level. The electrochemical reactions, related with those phenomena, are driven by an electrostatic potential difference between the electrostatic potential of the catalyst and the electrolyte at the reaction plane assumed to be situated at x = L, see Fig. 2. Such a quantity,

(y, t) − N (L, y, t),

(y, t) =

(1)

(y, t) = ϕ1 + ϕ2 ,

(2)

where ϕ1 is the potential drop related to the thickness of the compact layer and ϕ2 is the potential drop related with the dipolar nature of water molecules. The first term, ϕ1 , is expressed by means of Gauss law. The electric field created by the surface charge, (y, t), of the catalyst is given by

∂ (y, t) . =− CL ∂x

(3)

However as the thickness of the compact layer, d, is small, the previous equation can be written as ϕ1 = −

(y, t)

CL

d.

(4)

The calculation of ϕ2 is more complex. It has already been shown in previous work [26–28] that ϕ2 (y, t) is given by: ϕ2 (y, t) = −

anmax s (t) H2 O  sinh[X(y, t)]

DL

,

(5)

where a = 2e−(Eads /RT) in which Eads is the water chemical adsorption energy, R is the ideal gas constant and T the absolute temperature. DL and CL are the electric permittivities of the diffuse and the compact layer, respectively; s (t) is the fraction of the free sites;  is the dipolar moment of a water molecule and H2 O is the water activity; nmax is the maximum number of free sites per unit of area of the catalyst which is calculated by the surface area of the unit cell of MO2 (1 1 0) (M = Ir or Ru) for the anode side. For the cathode side it assumed to be 1/d2 as in Ref. [27]. X(y, t) is given by the solution of the following transcendental equation: a sinh[X] =  − X, (1/ H2 O ) + (1/ H2 O )( inter / s ) + a cosh[X]

(6)

with = (d3 /CL )A and = KTd3 /A2 ; inter is the surface coverage by the reaction intermediates, calculated by Eq. (29). Eq. (6) and its parameters are detailed in Ref. [27]. Finally from Eqs. (3) and (5) we have for (y, t) (y, t) = −

(y, t)

CL

d−

anmax s (t) H2 O  sinh[X(y, t)]

DL

.

(7)

The different quantities necessary to solve Eqs. (6) and (7) are described in the following paragraphs. The time evolution of the charge density, (y, t), is given by the following balance equation for the anode JeN− (y, t) − JFar A (y, t) =

∂ , ∂t

(8)

366

L.F.L. Oliveira et al. / Electrochimica Acta 110 (2013) 363–374

Fig. 1. Flowchart representation of the interconnection of the main equations used in the model.

and similarly for the cathode (9)

the electrochemical reaction kinetics described in the next section. JeN− (y, t) is the local current density at the catalytic surface and it is calculated through Eq. (48) (Section 2.3).

where JFar A (y, t) and JFar C (y, t) are respectively the anodic and cathodic faradaic current densities. The faradaic current is due to the charge transfer between the electrode and the electrolyte because of the reduction or oxidation of chemical species. In our model the faradaic current densities are calculated by means of

2.2.2. Kinetic models The overall electrochemical reactions of a PEMWEs are due to the oxidation and reduction half cell reactions which occur, respectively, at the anode and at the cathode.

JeN− (y, t) − JFar C (y, t) = −

∂ , ∂t

Fig. 2. Schematic view of our nanoscale model for the cathode PEMWEs electrode. The chemical species, considered are the H+ and H2 . The former is produced in the anode layer, which reacts at the catalyst surface to produce H2 . As a result, in the compact layer the species treated are H+ and H2 , both related to the HER.

L.F.L. Oliveira et al. / Electrochimica Acta 110 (2013) 363–374

For the case of PEMFCs, this methodology has been already Franco et al. [35]. In order to have the same type of data for the present model, we used DFT calculations to analyze the stability and the adsorption properties of water on the IrO2 (1 1 0) and RuO2 (1 1 0) surfaces which are two of the most used catalyst materials in PEMWEs devices. The details on our DFT calculations are being published elsewhere [36]. At present there is no consensus on the OER mechanism. From our DFT calculations, and based on other references [37], we propose a kinetic model for the OER for both materials. Concerning the HER, the mechanism used is simply the inverse of the well-known hydrogen oxidation reaction (HOR). We propose the following mechanism for the OER on the IrO2 surface: H2 O + s + Os → OHs + OHads

(10)

H2 O + s → H2 Oads

(11)

H2 Oads → OHads + e− + H+

(12)

OHads → Oads + e− + H+

(13)

H2 Oads + Oads → OOHads + e− + H+

(14)

OOHads → O2 + e− + H+

(15)

Reaction (10) comes from our finding of the spontaneous dissociation of water on IrO2 (1 1 0) described in [36]. We have also found that the oxygen bridge species, Obr [38], which are part of surface, participates in the reaction mechanism. In such reaction, Os refers to Obr species of the surface, and “s” represents a metallic active site. Through the dissociation of water two types of OH species will be produced. That is the reason why, in reaction (10) we distinguish two species, OHs and OHads . The former refers to the OH species formed through one hydrogen atom of the water molecule and the Obr species, and the latter refers to the OH bounded on top of the iridium atom and produced during the water dissociation. Therefore their behavior is different. The OHs species will not diffuse on the surface due to its strong bond with two iridium atoms. On the other hand the OHads species is more susceptible to diffuse on the surface. However, in the calculation of the evolution of the OH coverage, OH , presented below, we do not distinguish the previously mentioned species. Reaction (11) is related with the molecular adsorption of water, also a possible process encountered in the DFT calculations. The other elementary reactions, (11)–(15), included in our mechanism, are based on elementary reactions proposed in the literature [37]. For those reactions we cannot precise if the adsorption of the species takes place at the top of the metallic active site. This is because we have not performed yet detailed DFT calculations for the adsorption of OH and OOH. However, it seems to be the most probable, hence we keep the subscript “ads” for the intermediates participating in such reactions. We found that the numerical sensitivity of the backward reactions is very low, although physically important. We have then used the assumption of irreversibility. This has been done in order to reduce the number of parameters to fit with experimental data. Thus, we are dealing only with the forward reactions. The rates, of the above described reactions, are respectively as follows:

367

v5 = k5 H2 O O e(1−˛3 )f(y,t)

(20)

v6 = k6 OOH e(1−˛4 )f(y,t)

(21)

where H2 O is the water activity; ˛ is the dimensionless elementary electronic transfer coefficient; f = F/RT; j is the coverage of free sites or adsorbed species; H2 O is calculated as the total number of water dipoles adsorbed toward or opposed to the catalytic surface. The indexes 1–6 refer, respectively, to the reactions (10)–(15); ki denotes the rate parameters. Based on our findings, we have used in our simulations a similar mechanism to RuO2 . The main difference is the fact that the reaction which describes the spontaneous dissociation of water, reaction (10), is not present since this mechanism was not observed in our DFT calculations. Considering k1 = 0 in Eq. (16), the above described mechanism remains valid for the present case. The adsorption of water was found to be less exothermic on RuO2 (1 1 0) than on IrO2 (1 1 0) [36]. Hence this process should be faster on the latter material. Therefore, the value of k2 in Eq. (17) for RuO2 (1 1 0) should be lower than for IrO2 (1 1 0). For the HER in the cathode electrode, we assume the following three elementary reactions steps: H+ + s + e− → Hads

(22)

Hads + H+ + e− → H2 + s

(23)

2Hads → H2 + 2s

(24)

The following equations respectively describe the rates of the reactions (22)–(24):

v1 = K1 s CH+ e−˛1 f(y,t) ,

(25)

v2 = K2 Hads CH+ e˛2 f(y,t) ,

(26)

v3 = K3 H2 ads .

(27)

The surface coverages of each intermediate species, j , on the catalytic surface are calculated from the following balance equation d j

NAv = max n dt

 

vp −

p





vc

,

(28)

c

where NAv is the Avogadro number and vp and vc are, respectively, the production and consumption source terms of the species j due to the chemical reactions. By solving these equations, we can calculate the total coverage of intermediates species, inter , inter =



j .

(29)

j

For the case of the water dipoles directed toward the catalyst, the coverage is given by [26]  =

(a/2)eX , (1/ H2 O ) + (1/ H2 O )( inter / s ) + a cosh[X]

(30)

and for the opposed ones, ←

=

(a/2)e−X , (1/ H2 O ) + (1/ H2 O )( inter / s ) + a cosh[X]

(31)

v1 = k1 H2 O s O

(16)

v2 = k2 H2 O s

(17)

← H2 O =  + .

v3 = k3 H2 O e(1−˛1 )f(y,t)

(18)

The free sites, the coverage of the intermediate species and the water coverage on the catalyst surface are related through

v4 = k4 OH e(1−˛2 )f(y,t)

(19)

s + inter + H2 O = 1,

the previous equations allow us to calculate the water coverage, (32)

(33)

368

L.F.L. Oliveira et al. / Electrochimica Acta 110 (2013) 363–374

thus, we have for s s =

1 − inter . 1 + H2 O a cosh[X]

(34)

The chemical rates also allow us to calculate the faradaic currents, given in Eqs. (8) and (9). This is done through: JFar = F



vh

(35)

h

where F is the Faraday constant, and the sum is made over all the rates, vh , which involve a proton transfer. 2.2.3. Diffuse layer model The diffuse layer model handles the problem of the transport of chemical species to or from the catalyst surfaces. The fluxes of the chemical species, O2 and H2 , are calculated from the mass balances, by assuming a Fickian relationship for the flux expression, with a constant diffusion coefficient

⎧ N ⎨ ∂Ci = −∇ J N x i ∂t ⎩ N N

0≤x≤L

(36)

Ji = −Di ∇ x Ci ,

where N stands for nanoscale and i for H2 or O2 . The boundary conditions of the previous equation are calculated using the following equation: JiN (x = L, y, t) =



vp −



p

vc ,

(37)

c

in which vp and vc are the production and consumption source terms due to chemical reactions in which the species i is involved. For the case of H+ , the balance and the flux equations are given respectively by

⎧ N ⎨ ∂CH+ = −∇ J N x H+ ∂t ⎩ N N J

H+

0≤x≤L

= −DH+ (∇ x C

H+

(38)

N

− fC H+ ∇ x N ),

with boundary condition N JH + (x = L, y, t) =



vh ,

(39)

h

similarly to the Eq. (35), here vh refers to the rates which involve a proton transfer. The electrostatic potentials, N (x, y, t), for each electrode are calculated using the Poisson equation: F

DL

N (CH +

− CSO3− ) =

−∇ 2x N ,

(40)

where CSO3− is the concentration of sulfonic acid groups of the ionomer, which is assumed to be homogeneous across the diffuse layer thickness. Further model versions, should describe heterogeneous sulfonic acid group spatial distributions, accordingly to recent molecular dynamics simulation results[14]. N , and the potential, N (x, y, The concentrations, C N+ , CON and CH H 2 2 t), are connected with Eqs. (41) and (46) of the microscale models (Section 2.3). These values will also give us the necessary data to be used for the calculation of the reaction rates presented in the previous section. 2.3. Microscale model

Fig. 3. Schematic view of our microscale model for the cathode electrode.

according to the y and z directions. These directions reflect the main directions of the species transport. The reactants, O2 and H2 , are assumed to diffuse mainly in the perpendicular direction to the catalyst surface toward the water liquid phase within the pores of the catalyst material. The protons are assumed to be transported in parallel from the catalyst surface to the membrane or vice versa, depending of the electrode under consideration. Finally, electrons are also transported parallel from the catalyst surface toward the external electrical circuit or vice versa, depending again on the electrode considered. To account for the O2 (or H2 ) transport through the ionomer layer in the z direction, we use, as in the nanoscale model, the following material balance equations:

⎧ M ⎨ ∂Ci = −∇ J M 0 ≤ z ≤ e z i NI ∂t ⎩ M M

where M stands for microscale and eNI is the thickness of the ionomer layer. The molar flux of O2 (or H2 ) is assumed to be well represented by the Fick’s model with a constant diffusion coefficient. Two boundary conditions have to be posed to solve Eq. (41). In order to account for a possible mass transfer limitation within the liquid phase, a simplified view of the situation is to use the classical film model [39]: W,eq

JiM (y, z = 0, t) = ˇi (Ci

− CiW,∞ )

(42)

W,eq

where Ci and CiW,∞ are, respectively, the concentration of the species i in the liquid phase at the ionomer surface and far from it and ˇi is the transfer coefficient. As far as the ionomer surface is concerned, it is assumed to be at equilibrium with the liquid phase in its immediate vicinity. And, as we assume that O2 (or H2 ) in water is in equilibrium with O2 (or H2 ) in gas phase, COW,∞ is calculated 2 from the Henry law: CiW,∞ =

p , kHi

(43)

where p is the partial pressure in the gas phase and kHi is the Henry constant for O2 or H2 . The oxygen and hydrogen concentrations in the ionomer at z = 0 are assumed to be at thermodynamic equilibrium with the oxygen and hydrogen concentration in the water phase eq

The microscale structure is simplified in a cylindrical-like structure for the purpose of the simulation reported in this work. Fig. 3 illustrates a cut view of this approach for the cathode electrode. The model is a combination of two interconnected 1-D models,

(41)

Ji = −Di ∇ z Ci

W,eq

CiM (z = 0, t) = Ci (T, Ci eq

W,eq

)

(44)

) is the O2 (or H2 ) solubility within the ionomer where Ci (T, Ci when it is in contact and at equilibrium with water containing O2 W,eq (or H2 ) according to the concentration Ci . At z = eNI , the flux of O2

L.F.L. Oliveira et al. / Electrochimica Acta 110 (2013) 363–374

(or H2 ) at the microscale level, JiM (y, z = eNI , t), is given by the flux at the nanoscale model, JiN (x = 0, y, t), by assuming the continuity of the fluxes at the interface between the two modeling scales: JiM (y, z

= eNI , t) =

JiN (x

= 0, y, t).

(45)

The proton transport within the microscale model is assumed to occur in the y direction and at steady-state. The corresponding balance equation is

∇y

iM+ H

Selectrode

M = SH 0 ≤ y ≤ eCA + (y, t)

(46)

M (y, t) is the local where eCA is the thickness of the catalyst and SH+ volumetric source term of protons of the microscale model, which is given by the nanoscale model through the following equation: M N SH + (y, t) = CA · J + (x = 0, y, t), H

(47)

where  CA is a geometrical parameter representing the specific area of the active surface (units m2 m−3 ) allowing to express a volumetric source term for the balance equation, Eq. (46), from the value of the flux of protons given by the nanoscale model. Similarly to the transport of protons, the electronic transport is given by

∇y

ieM− Selectrode

= SeM− (y, t) 0 ≤ y ≤ eCA ,

(48)

where SeM− (y, t) = CA · JeN− (y, t). JeN− (x

(49) J N+ (x H

We consider that = 0, y, t) = = 0, y, t), no electronic accumulation at x = 0. Eq. (48) has the following boundary conditions:

369

Table 1 List of parameters related with the cathode electrode. Parametersa

Unit

References

K1 = 10−5 K2 = 10−6 K3 = 10−1 ˛1 = 0.5 ˛2 = 0.5 WPt = 4 Eads = 1000

mol s−1 m−2 mol s−1 m−2 mol s−1 m−2 Dimensionless Dimensionless mg cm−2 J mol−1

[27] [27] [27] [27] [27] Assumed [27]

a

Unless stated otherwise.

2.4. Membrane As aforesaid, we assume that the membrane is impermeable to O2 and H2 . That is, only protons and water cross it. The proton transport is represented by Ohm’s law, as we are assuming the electroneutrality in the electrolyte phase: Cbulk (t) = Abulk (t) − RM · iH+ (t)

(55)

where Cbulk (t) and Abulk (t) are, respectively, the electrostatic potentials of the cathode and of the anode at the electrolyte bulk, (x = 0); iH+ , the protonic current, which is considered equal to the external electric current, I(t); and RM is the membrane resistance, which is calculated using RM =

EM Selectrode gH+

(56)

where EM is the membrane thickness and gH+ , the protonic conductivity, is given by [41]:

ieM− (y = eCA , t) = 0,

(50)

gH+ = (0.46 − 0.25)e[−1190((1/T )−(1/298.15))]

ieM− (y = 0, t) = I(t)

(51)

where  is the local number of water molecules per sulfonate site and T the temperature.

where I(t) is the external current applied to the system. The value of the structural parameter,  CA , can be estimated either from experimental data or by models such as Monte Carlo or coarse grained molecular dynamics [40] and it is related with the porosity of the metal oxide materials. As our model is written in a single mesh approach no discretiN (y, t) can be zation along the electrode thickness is performed, Je− calculated through: N Je− (y, t) =

I(t) SCA

(52)

where I(t) is the global current applied to system and SCA is calculated for the cathode electrode manufacture parameters using SC =

3 WPt rp Pt Ral

2.5. Cell potential Finally, the cell potential of the PEMWEs, UWE , is calculated from the catalyst electrostatic potential, (y, t), of each electrode, through the following equation UWE =

C (y, t) −

A (y, t),

(58)

where C (y, t) and A (y, t) are, respectively, the catalyst electrostatic potential of the cathode and of the anode. 3. Results and discussion

(53)

where rp is the average radius of the Pt/C particle, Pt it is the platinum density, WPt is the platinum loading in the electrode, and Ral is the thickness of the active layer. For the anode electrode we use values of electrochemical surface area (ECSA), obtained by Brunauer–Emmett–Teller (BET) experiments. For the cathode side, such type of measurements will not be useful because the Pt nanoparticles are deposed on larger carbon nanoparticles. Therefore we use the following equation SA = massCA · ECSA

(57)

(54)

where ECSA value is given in Table 2 and massCA is the total mass of catalyst on the anode, which is calculated from its loading and the total area of the electrode.

The results of the simulations presented here were numerically solved using the MEMEPhys simulation package which is coded in C language and Matlab/Simulink [26,27]. The C codes are embedded into the Matlab/Simulink environment thanks to the use of the socalled S-functions. The simulations have been carried out using the Simulink solver ode23s with a relative tolerance of 10−3 . For the purposes of this paper, only one mesh was used for the simulation of each electrode, and a quasi-stationary formulation has been used for the calculation of the Poisson Nernst Planck equations in the diffuse layer regions [49]. This section is organized as follows. Firstly we present results concerning the coverage of intermediates species at the surface of the anode catalyst materials. Secondly a sensitivity study for the polarization curves of the whole MEA is presented.

370

L.F.L. Oliveira et al. / Electrochimica Acta 110 (2013) 363–374

Table 2 List of general parameters used for the simulation of both materials. Parametersa

Unit

References

CL = 60 DL = 100

C2 J−1 m−1 C2 J−1 m−1 m m m

[42] [42] Assumed Assumed Assumed

T

m2 s−1

[27]

T

atm m3 mol−1

[43]

atm m3 mol−1

[27]

m2 g−1 cm2 S m−1 ␮m K Pa

[44] Assumed Assumed Assumed Assumed Assumed

eNI = 10−9 eCA = 5 ×10−6 d = 2 ×10−10 DO2 = 3.1 × 10−7 × exp eq

CO = 1.34 × 102 × exp

 −2768  −1540

2

KHO = 5.08 × 106 × exp

 −498

2

T

ECSA = 105 Selectrode = 25 gH+ = 9.3 EM = 127 T = 333 p = 105 a

Fig. 5. Calculated water, OH and O coverages as a function of the applied current.

Unless stated otherwise.

Table 3 List of parameters related to the IrO2 electrode. Parametersa

Unit

References

k1 = 10−1 k2 = 10−2 k3 = 10−2 k4 = 10−2 k5 = 2 ×10−2 k6 = 5 ×102 ˛1 = 0.5 ˛2 = 0.5 ˛3 = 0.5 ˛4 = 0.5 WIrO2 = 0.64 Eads = −137

mol s−1 m−2 mol s−1 m−2 mol s−1 m−2 mol s−1 m−2 mol s−1 m−2 mol s−1 m−2 Dimensionless Dimensionless Dimensionless Dimensionless mg cm−2 K J mol−1

Assumed Assumed Assumed Assumed Assumed Assumed Assumed Assumed Assumed Assumed Assumed [36]

a

j/mA cm-²

Unless stated otherwise.

3.1. Coverage study for the anode electrode In Table 1 the conditions concerning the cathode electrode are summarized. In Table 2 the values of the parameters used in the results presented here and in the following sections are presented. In Table 3 we list the parameters related to IrO2 and in Fig. 4 we show the corresponding calculated polarization curve. In Fig. 5 it is presented a comparison of the calculated water, OH, and O coverages as a function of the applied current and corresponding to the conditions used for the calculation of the

polarization curve shown in Fig. 4. This result shows that water and O coverage decrease as the current density increases. For the coverage of OH the trend is the opposite. That is, OH coverage is higher at higher current densities. This is due the fact that the production of OH is directly connected with the dissociation of water and its interaction with the Obr species. In other words, the coverage of water and O decreases with the current density because the amount of produced OH increases. Similarly to the coverage of OH, the coverage of OOH also increases with the current densities, as shown in Fig. 6 but the sensitivity is lesser as can be seen by the values of the coverage. This suggests that the formation of OOH species on the catalyst surface, through reaction (14), is much lower than the formation of OH. DFT calculations will be a useful tool to deeply analyze such interaction. In order to study the temperature impact on the coverage, we outline in Figs. 7–9, respectively the coverage of water, O and OH species on the catalyst. For the coverage evolution of water and O, desorption is favored at higher temperatures. On the other hand, higher coverage of adsorbed OH is expected at higher temperature. This, although counterintuitive, can be explained by the fact that the adsorption of OH is, as discussed above, directly connected to the adsorption of O and water. In this case the increase in the OH coverage is also connected with the temperature dependency of the reaction rates, i.e. at higher temperatures higher is the production of OH.

j/mA cm-² Fig. 4. Calculated polarization curve for the IrO2 electrode.

Fig. 6. Calculated OOH coverage as function of the applied current.

L.F.L. Oliveira et al. / Electrochimica Acta 110 (2013) 363–374

371

Table 4 List of parameters related with RuO2 electrode. Parametersa

Unit

References

k1 = 0 k2 = 10−3 k3 = 10−2 k4 = 10−2 k5 = 2 ×10−2 k6 = 5 ×102 ˛1 = 0.5 ˛2 = 0.5 ˛3 = 0.5 ˛4 = 0.5 WRuO2 = 6 Eads = −99

mol s−1 m−2 mol s−1 m−2 mol s−1 m−2 mol s−1 m−2 mol s−1 m−2 mol s−1 m−2 Dimensionless Dimensionless Dimensionless Dimensionless mg cm−2 K J mol−1

Assumed Assumed Assumed Assumed Assumed Assumed Assumed Assumed Assumed Assumed Assumed [36]

a

Unless stated otherwise.

j/mA cm-² Fig. 7. Calculated water coverage at different temperatures.

j/mA cm-² Fig. 8. Calculated O coverage at different temperatures.

3.2. Polarization curves for RuO2 and IrO2 based catalysts In this section we show the results concerning the simulation of a PEMWEs anode with RuO2 as catalyst material. In Table 4 it is summarized the values of the parameters used for the simulation of the RuO2 -based catalyst performances. The coverage of the intermediate species evolves in a similar way as the ones presented for IrO2 , so that they are not presented

j/mA cm-² Fig. 9. Calculated OH coverage at different temperatures.

here. In Fig. 10 it is presented the calculated polarization curve for each material. In contrast to the experimental knowledge, these curves show a higher performance for the IrO2 than for the RuO2 based catalyst. Although the adsorption of water on the oxide surfaces is an important issue, is not sufficient to describe the whole mechanism. Therefore, further DFT calculations should be performed to address this issue. From our model we can suggest which of the steps of the proposed elementary reaction mechanism for OER should be first investigated. With this aim, an extensive sensitivity study with respect to the values of the rate parameters was made. We have found out that reaction (14) is the most sensitive one. Fig. 11 shows a second comparison between RuO2 and IrO2 -based catalyst, using the same values of k s used in the previous presented result but with a different value of k5 for RuO2 . This result suggests that the production of OOH through the interaction between water and adsorbed O species on the RuO2 surface must play a important role in the mechanism of the OER. 3.2.1. Sensitivity studies Our model allows us to perform a variety of studies. Such phenomena are, for example, the study of the concentration and the flows of species along the thickness of the electrodes, the potentials drop across the compact layers, (y, t), or on the diffuse layer, N (x, y, t) and the coverages of intermediate species. The capabilities to simulate one of those phenomena – the coverage of the intermediate species – have already been demonstrated in the

Fig. 10. Comparison of the evolution of the calculated anode potential for the IrO2 and RuO2 based catalysts. Once the electrochemical mechanisms were obtained from DFT calculations, we use k1 = 0 and k2 = 10−3 mol s−1 m−2 for the corresponding kinetic constants.

372

L.F.L. Oliveira et al. / Electrochimica Acta 110 (2013) 363–374

Fig. 11. Comparison of the evolution of the calculated anode potential for the IrO2 and RuO2 based catalysts. This simulation was performed with the aim to investigate which of the proposed elementary steps for the OER most influence in the higher performance of RuO2 when compared with IrO2 .

previous section. In this section we have chosen to concentrate our results on the sensitivity to different parameters of the PEMWEs MEA performance. This is done by the analysis of the calculated polarization curve. The i–V curves presented in this section were generated by taking the values of the potential, for a given current density, after the steady state is reached. 3.2.1.1. Sensitivity to temperature. The performance of the PEMWEs increases as the temperature increases (i.e. at a given current, potential decreases with temperature, meaning that less input electric power is necessary to perform the electrolysis). It has been already shown experimentally [45] that the operating temperature influences substantially the PEMWEs performance, the experimental trend is successfully reproduced by our model, Fig. 12. Such behavior is directly related to the kinetic rates mainly of the anode electrode. Once again the importance of the understanding of nanoscale phenomena is highlighted here. In contrast with the experimental curves, our simulated i–V curves do not show an increase of the potential ohmic slope with temperature. This should be connected with the fact that we have temperature dependencies in different parts of our model, such

Fig. 12. Calculated polarization curves for different temperatures.

Fig. 13. Calculated polarization curves for different IrO2 loadings.

as in the electrochemical sub-models and the material transport sub-models. Hence, the prediction of the ohmic loss for different temperatures, is compensated by other phenomena which is also dependent on temperature. 3.2.1.2. Sensitivity to electrode properties. In Fig. 13 and 14 we report a sensitivity study of the i–V curves to the catalyst loading in the anode and in the cathode, respectively. As expected, the value of the cell potential increases as the catalyst loading value decreases. This behavior is verified in both electrodes. Due to the high cost of the catalyst materials, many efforts have been made for the design of new catalysts that have both high performance and lower cost, compared with the usual catalysts. If the appropriate elementary kinetic model and associated ab initio data are known, our model could open the possibility of simulation for the testing of new materials under real conditions of water electrolysis. Interestingly the increase of the potential ohmic slope, commented above, is present for the anode case but not for the cathode. This difference, in real systems, is probably due to the electronic support. In other words, when varying the platinum loading, the amount of carbon is assumed constant, hence the electronic percolation through carbon nanoparticules is unchanged. Whereas, when varying the IrO2 loading the electronic percolation is also changed, as IrO2 is not deposited in a further material. Therefore, for the

Fig. 14. Calculated polarization curves for different Pt loadings.

L.F.L. Oliveira et al. / Electrochimica Acta 110 (2013) 363–374

Fig. 15. Calculation of polarization curves for different anode ECSA values.

anode case, once the amount of catalyst material is decreased, the electronic percolation will also decrease and in turn a higher ohmic loss will be detected. Another feature of the electrode that we can analyze is the ECSA, Fig. 15, that is related with the morphology of the catalyst material. In this case the analysis is made only for the anode catalyst. For the cathode case we are not using the values of the ECSA. 4. Conclusions and perspectives In this paper we have proposed a multiscale model of a PEMWEs MEA. This model includes a detailed description of the elementary kinetics, a description of the non-equilibrium behavior of the nanoscale catalyst–electrolyte interface, and a microstructuralresolved description of the transport of charges and chemical species at the micro and mesoscales along the electrodes. This new multiscale physical transient model gives the possibility to analyze the sensitivity of the PEMWEs behavior under a large set of operating conditions and material parameters (e.g. catalyst loading). This model does not use the empirical Butler–Volmer equation but is based on an elementary kinetic approach. The model has been built up in such a manner that it allows us to compare macroscopic observables, such as polarization curves, as well as other interesting parameters at other scales such as the coverage of the species at the catalyst surface. Encouraging results were obtained which agree well with some experimental trends reported in the published literature [45]. As previously mentioned our model relies on different assumptions leading to some limitations. In our opinion, our model exhibits three main limitations. The first one is the fact that we have assumed an uniform and constant temperature of the MEA. Within the framework of our geometrical description of the system, energy balances could be easily included in order to calculate the temperature. The second one is due to the simple geometrical representation of the electrodes by the three cylinders that we have adopted. Such a simplified representation guaranties that the three phases (ionomer, water, catalyst and electronic conductor) are continuous but leads to an oversimplification of the catalyst active surface description. For example, this description does not allow to take into account the oxides porosity. A precise catalyst active surface area modeling is important since the link between the elementary kinetics description and the macroscale performances is based on this parameter. The third limitation is due to the simplification that we have made to describe the mass transfer between the ionomer and the water. For high current densities, bubble formation occurs that can be hardly described by a global mass transfer coefficient. Despite these simplifications, the model allows

373

providing trends in agreement with experimental knowledge, and thus provides an interesting basis for the development of more complex models. With relation to the multiscale model, another important phenomenon to be studied is the membrane degradation due to H2 and O2 cross-permeation through the membrane. A second important phenomenon, which was not described in this paper, is the formation of gas bubbles. Such phenomenon could impact the PEMWEs performance and durability at high current densities. Hence, a model to study the formation of bubbles and its further inclusion in the multiscale model, would be of an important value and will be the subject of one of our future publications. One example in which our model could contribute to the improvement of the MEAs in PEMWEs concerns the possibility of existing catalyst screening by simulation under real conditions of water electrolysis. This type of simulation is of paramount importance because, as already discussed the high cost of the catalytic materials is one of the main drawbacks of PEMWEs technology. Therefore simulations could suggest specific routes for the experiments avoiding then the need of extensive benchmark tests. One can also plan to sketch the ideal catalyst properties with respect to the other properties of the MEA. Thanks to the modular structure of these types of models it is also possible to develop a system level model. Coupling the PEMWEs with a hydrogen storage system and the PEMFCs. Acknowledgements Valuable discussions with Dr. D. Loffreda (ENS-Lyon) are gratefully acknowledged. We also thank Dr. P. Giunta, who worked as engineer at Dr. Franco’s group, for proof reading. Appendix A. List of symbols

a f Ci CH+ CSO−

a = 2 exp(− Eads /RT) f = F/RT hydrogen or oxygen concentration (mol m−3 ) proton concentration (mol m−3 ) sulfonic acid groups concentration (mol m−3 )

d Di eNI eCA eDL EM ECSA F gH+ I(t) JO2 JH+ JFar ki nmax

thickness of the compact layer (m) diffusion coefficient of hydrogen or oxygen (m2 s−1 ) microscale impregnated ionomer thickness (m) thickness of the electrode (m) thickness of the diffuse layer (m) membrane thickness (m) eletrochemical surface area (m2 g−1 ) Faraday constant (96485) (C mol−1 ) proton conductivity (S m−1 ) external electric current (A) oxygen molar flux in the ionomer (mol s−1 m−2 ) protons molar flux in the ionomer (mol s−1 m−2 ) total faradaic current (A m−2 ) standard reaction rate parameter (mol s−1 m−2 ) maximal quantity of free sites per unit of area of catalyst phase (m−2 ) Avogrado number (6.02 × 1023 ) (mol−1 ) partial pressure of hydrogen or oxygen (Pa) active surface area of the electrode (m2 ) ideal gas constant (8.34) (J K−1 mol−1 ) membrane resistance () geometrical surface area of the electrode (m2 ) absolute temperature (K) water electrolyzer potential (V) reaction rate (mol s−1 m−2 )

3

NAv pi S R RM Selectrode T UWE

vi

374

W ˛i

0 CL DL  i    

L.F.L. Oliveira et al. / Electrochimica Acta 110 (2013) 363–374

catalyst loading (mg cm−2 ) elementary reaction step electronic transfer coefficients (dimensionless) electric permittivity in the vacuum (8.85 × 10−12 ) (C2 J−1 m−1 ) electric permittivity in the compact layer (C2 J−1 m−1 ) electric permittivity in the diffuse layer (C2 J−1 m−1 ) electrostatic potential difference between the catalyst and the electrolyte phases (V) coverage fraction of free site or adsorbed species (dimensionless) local number of water molecules per sulfonate site (dimensionless) dipolar moment of a water molecule (D) surface charge density on the catalyst surface (C m−2 ) electrostatic potential in the electrolyte phase (V) electrostatic potential in the catalyst phase (V)

References [1] R. de Levie, The electrolysis of water, J. Electroanal. Chem. 476 (1999) 92. [2] A. Marshall, B. Borresen, G. Hagen, M. Tsypkin, R. Tunold, Hydrogen production by advanced proton exchange membrane (PEM) water electrolyser – reduced energy consumption by improved electrocatalysis, Energy 32 (2007) 431–436. [3] E. Rasten, G. Hagen, R. Tunold, Electrocatalysis in water electrolysis with solid polymer electrolyte, Electrochim. Acta 48 (2003) 3945–3952. [4] P. Millet, R. Ngameni, S.A. Grigoriev, N. Mbemba, F. Brisset, A. Ranjbari, C. Etiévant, PEM water electrolyzers: from electrocatalysis to stack development, Int. J. Hydrogen Energy 35 (2010) 5043–5052. [5] W.T. Grubb, Batteries with solid ion-exchange electrolytes, J. Electrochem. Soc. 106 (1959) 275–278. [6] U.D. of Energy, Solid polymer electrolyte water electrolysis technology development for large-scale hydrogen production, 1981, DOE/ET/26 202-1. [7] F. Barbir, PEM electrolysis for production of hydrogen from renewable energy sources, Solar Energy 78 (2005) 661–669. [8] Giner, 2013. Available on http://www.ginerinc.com [9] ProntoOnSite, 2013. Available on http://www.protonsite.com [10] Ceth, 2013. Available online on www.ceth.fr [11] Hydrogenics, 2013. Available on http://www.hydrogenics.com [12] P. Millet, M. Pineri, R. Durand, New solid polymer electrolyte composite for water electrolysis, J. Appl. Electrochem. 19 (1989) 162–166. [13] A.A. Franco, Y. Suzue, R.F. de Morais, A. Kachmar, T. Mashio, A. Ohma, K. Shinohara, Multiscale modeling of electrochemical double layer formation on Pt electrodes covered by Nafion, in: Abstract 226, 221st ECS Meeting, Seattle, WA, 2012. [14] D.D. Borges, S. Mossa, K. Malek, G. Gebel, A.A. Franco, Effect of surface hydrophilicity on the formation of Nafion thin films inside PEMFCs catalyst layers: a computational study, ECS Trans. 45 (2013). [15] S. Song, H. Zhang, X. Ma, Z. Shao, R.T. Baker, B. Yi, Electrochemical investigation of electrocatalysts for the oxygen evolution reaction in PEM water electrolyzers, Int. J. Hydrogen Energy 33 (2008) 4955–4961. [16] H. Ito, T. Maeda, A. Nakano, H. Takenaka, Properties of nafion membranes under PEM water electrolysis conditions, Int. J. Hydrogen Energy 1 (2011) 1–14. [17] S. Becker, V. Karri, Predictive models for PEM-electrolyzer performance using adaptive neuro-fuzzy inference systems, Int. J. Hydrogen Energy 35 (2010) 9963–9972. [18] O. Atlam, M. Kolhe, Equivalent electrical model for a proton exchange membrane (PEM) electrolyser, Energy Convers. Manage. 52 (2011) 2952–2957. [19] S. Busquet, C. Hubert, J. Labbe, D. Mayer, R. Metkemeijer, A new approach to empirical electrical modelling of a fuel cell, an electrolyser or a regenerative fuel cell, J. Power Sources 134 (2004) 41–48. [20] P. Choi, D.G. Bessarabov, R. Datta, A simple model for solid polymer electrolyte (SPE) water electrolysis, Solid State Ionics 175 (2004) 535–539. [21] K. Onda, T. Murakami, T. Hikosaka, M. Kobayashi, R. Notu, K. Ito, Performance analysis of polymer-electrolyte water electrolysis cell at a small-unit test cell and performance prediction of large stacked cell, J. Electrochem. Soc. 149 (2002) A1069–A1078.

[22] K.S. Agbli, M.C. Péra, D. Hissel, O. Rallières, C. Turpin, I. Doumbia, Multiphysics simulation of a PEM electrolyser: energetic macroscopic representation approach, Int. J. Hydrogen Energy 36 (2011) 1382–1398. [23] F. Marangio, M. Santarelli, M. Calì, Theoretical model and experimental analysis of a high pressure PEM water electrolyser for hydrogen production, Int. J. Hydrogen Energy 34 (2009) 1143–1158. [24] S.A. Grigoriev, A.A. Kalinnikov, P. Millet, V.I. Porembsky, V.N. Fateev, Mathematical modeling of high-pressure PEM water electrolysis, J. Appl. Electrochem. 40 (2010) 921–932. [25] L.F.L. Oliveira, S. Laref, E. Mayousse, C. Jallut, A.A. Franco, A multiscale physical model for the transient analysis of PEM water electrolyzer anodes, Phys. Chem. Chem. Phys. 14 (2012) 10215–10224. [26] A.A. Franco, P. Schott, C. Jallut, B. Maschke, A dynamic mechanistic model of an electrochemical interface, J. Electrochem. Soc. 153 (2006) A1053–A1061. [27] A.A. Franco, P. Schott, C. Jallut, B. Maschke, A multi-scale dynamic mechanistic model for the transient analysis of PEFCs, Fuel Cells 07 (2007) 99–117. [28] A.A. Franco, Un modèle physique multiéchelle de la dynamique électrochimique dans une pile à combustible à électrolyte polymère, Université Claude Bernard Lyon 1, 2005 (Ph.D. Thesis). [29] A.A. Franco, A multiscale modeling framework for the transient analysis of PEM fuel cells – from the fundamentals to the engineering practice, Université Claude Bernard Lyon 1, 2010, Habilitation (H.D.R.) manuscript. [30] A.A. Franco, PEMFCs degradation modeling and analysis, in: C. Harting, C. Roth (Eds.), Polymer Electrolyte Membrane and Direct Methanol Fuel Cell Technology (PEMFCs and DMFCs) – vol. 1: Fundamentals and Performance, Woodhead, Cambridge, UK, 2012. [31] A.A. Franco, A multiscale physical model of electrochemical storage systems, ECS Trans. 45 (2013) 11–19. [32] A.A. Franco, Multiscale modeling and numerical simulation of rechargeable lithium ion batteries: concepts, methods and challenges, RSC Adv. 3 (32) (2013) 13027–13058. [33] A.A. Franco, M. Tembely, Transient multiscale model of aging mechanisms in a PEFC cathode, J. Electrochem. Soc. 154 (2007) B712–B723. [34] A.A. Franco, S. Passot, P. Fugier, C. Anglade, E. Billy, L. Guétaz, N. Guilet, E. de Vito, S. Mailley, Ptx Coy catalysts degradation in PEFC environments: mechanistic insights, J. Electrochem. Soc. 156 (2009) B410–B424. [35] R.F. Morais, P. Sautet, D. Loffreda, A.A. Franco, A multiscale theoretical methodology for the calculation of electrochemical observables from ab initio data, Electrochim. Acta (2011) 10842–10856. [36] L.F.L. Oliveira, A.A. Franco, C. Jallut, D. Loffreda, Paper in preparation, 2013. [37] J. Rossmeisl, Z.-W. Qu, H. Zhu, G.J. Kroes, J.K. Norskov, Electrolysis of water on oxide surfaces, J. Electroanal. Chem. 607 (2007) 83–89. [38] K.R.Q. Sun, M. Sheffler, Hydrogen adsorption on RuO2 density-functional calculations, Phys. Rev. B 70 (2004) 235402. [39] R.B. Bird, W.E. Stewart, E.N. Lightfoot, Transport Phenomena, John Wiley and Sons, New Yark, 2002. [40] A.A. Franco, K. Malek, Microstrucure-based modeling of aging mechanisms in catalyst layers of polymer electrolyte fuel cells, J. Phys. Chem. B 115 (2011) 8088–8101. [41] F. Meier, G. Eigenberger, Transport parameters for the modelling of water transport in ionomer membranes for PEM-fuel cells, Electrochim. Acta 49 (2004) 1731–1742. [42] J.O. Bockris, S. Khan, Surface Electrochemistry: A Molecular Level Approach, Plenum Press, New Yark; London, 1993. [43] H. Ito, T. Maeda, A. Nakano, Y. Hasegawa, N. Yokoi, C.M. Hwang, M. Ishida, A. Kato, T. Yoshida, Effect of flow regime of circulating water on a proton exchange membrane electrolyzer, Int. J. Hydrogen Energy 35 (2010) 9550–9560. [44] E. Mayousse, F. Maillard, F. Fouda-Onana, O. Sicardy, N. Guillet, Synthesis and characterization of electrocatalysts for the oxygen evolution in PEM water electrolysis, Int. J. Hydrogen Energy 36 (2011) 10474–10481. [45] L. Ma, S. Sui, Y. Zhai, Investigations on high performance proton exchange membrane water electrolyzer, Int. J. Hydrogen Energy 34 (2009) 678–684. [46] A. Marshall, B. Børresen, G. Hagen, M. Tsypkin, R. Tunold, Electrochemical characterisation of Irx Sn1−x O2 powders as oxygen evolution electrocatalysts, Electrochimica Acta 51 (2006) 3161–3167. [47] A.A. Franco (Ed.), Polymer Electrolyte Fuel Cells: Science, Applications and Challenges, Pan Stanford Publishing, Singapore, 2013. [48] A.A. Franco, Multiscale modeling of electrochemical devices for energy conversion and storage, in: R. Savinell, K.I. Ota, G. Kreysa (Eds.), Encyclopedia of Applied Electrochemistry, Springer, UK, 2013. [49] A.A. Franco, M. Gerard, Transient model of carbon catalyst-support corrosion in a PEFC: multi-scale coupling with electro-catalysis and impact on performance degradation, J. Electrochem. Soc. 155 (4) (2008) B367.