Degradation forecast for PEMFC cathode-catalysts under cyclic loads

Journal of Power Sources 359 (2017) 611e625

Contents lists available at ScienceDirect

Journal of Power Sources journal homepage: www.elsevier.com/locate/jpowsour

Degradation forecast for PEMFC cathode-catalysts under cyclic loads M. Moein-Jahromi a, b, M.J. Kermani a, b, *, S. Movahed a a

Department of Mechanical Engineering, 424 Hafez Avenue, Amirkabir University of Technology (Tehran Polytechnic), Tehran 15875-4413, Iran Fuel Cell Laboratory, Department of Mechanical Engineering, 424 Hafez Avenue, Amirkabir University of Technology (Tehran Polytechnic), Tehran 158754413, Iran

b

h i g h l i g h t s
Cathode Catalyst Layer (CL) degradation is investigated under cyclic load.
Electro-Chemical Surface Area (ECSA) loss of CL is predicted during cyclic load.
Coarsening of agglomerate particle in the CL is modeled by Ostwald ripening.
ECSA degradation is analogized with carbon steel fatigue phenomenon.
A set of parametric study is done in order to rank the most influential parameters.

a r t i c l e i n f o

a b s t r a c t

Article history: Received 10 May 2017 Accepted 29 May 2017

Degradation of Fuel Cell (FC) components under cyclic loads is one of the biggest bottlenecks in FC commercialization. In this paper, a novel experimental based algorithm is presented to predict the Catalyst Layer (CL) performance loss during cyclic load. The algorithm consists of two models namely Models 1 and 2. The Model 1 calculates the Electro-Chemical Surface Area (ECSA) and agglomerate size (e.g. agglomerate radius, rt,agg) for the catalyst layer under cyclic load. The Model 2 is the already-existing model from our earlier studies that computes catalyst performance with fixed structural parameters. Combinations of these two Models predict the CL performance under an arbitrary cyclic load. A set of parametric/sensitivity studies is performed to investigate the effects of operating parameters on the percentage of Voltage Degradation Rate (VDR%) with rank 1 for the most influential one. Amongst the considered parameters (such as: temperature, relative humidity, pressure, minimum and maximum voltage of the cyclic load), the results show that temperature and pressure have the most and the least influences on the VDR%, respectively. So that, increase of temperature from 60  C to 80  C leads to over 20% VDR intensification, the VDR will also reduce 1.41% by increasing pressure from 2 atm to 4 atm. © 2017 Published by Elsevier B.V.

Keywords: Catalyst layer degradation Cyclic load Electro-chemical surface area Ostwald ripening Catalyst layer performance

1. Introduction Cost-effectiveness and sustainability are the two most significant factors for power generation systems such as Fuel Cells (FC). Durability and cost are the most important challenges for commercialization of the FC systems [1]. In fact, there should be a tradeoff between three key factors: performance, cost and endurance. PEM Fuel Cells (PEMFCs), as a serious alternative of internal combustion engines, has higher efficiency and power density and

* Corresponding author. Department of Mechanical Engineering, 424 Hafez Avenue, Amirkabir University of Technology (Tehran Polytechnic), Tehran 158754413, Iran. E-mail addresses: [email protected] (M. Moein-Jahromi), mkermai@aut. ac.ir (M.J. Kermani), [email protected] (S. Movahed). http://dx.doi.org/10.1016/j.jpowsour.2017.05.102 0378-7753/© 2017 Published by Elsevier B.V.

no emission of greenhouse gases. However, according to the 2020 roadmap of the Department Of Energy (DOE) of USA, the lifetime of the PEMFCs should be up to 5000 h to be commercialized in the automotive industries [1]. Achieving this target needs accurate and sophisticated models of predicting the degradation rate of PEMFCs. Catalyst Layers (CL), as the sites of electrochemical reaction, are a vital part of the PEMFCs that play an important role on their durability. Transient operation procedure and poisoning with some gases such as sulfur and CO are the two main reasons of ElectroChemical Surface Area (ECSA) decline. The transient operation of the PEMFCs has considerable effects on Platinum (Pt) particle agglomeration, dissolution, and consequently, ECSA decline [2]. The transient operations procedures of PEMFC for automotive application can be categorized to idling, Open Circuit Voltage (OCV), cyclic load and startup/shutdown operation [2]:

612

M. Moein-Jahromi et al. / Journal of Power Sources 359 (2017) 611e625

Idling operation: a little but constant power is drawn from the PEMFC as an auxiliary load (constant voltage about ~0.9 V vs. RHE), which may take for a long time. The idling operation without any fluctuating in voltage is a more durable operating procedure than the cyclic one (longer lifetime). Its reason will be discussed in the following sections. Wilson et al. tested a PEMFCs with 0.5 and 0.75 V constant cell voltages to examine the effect of idling operation on the ECSA decline both in anode and cathode sides [3]. The results showed that if the idling potential is 0.5 V vs. RHE, the ECSA decline percentage of anode and cathode will be 33% and 57% of its initial value after 3750 h, respectively. However, by increasing idling voltage to 0.75 V, the ECSA decline percentage of anode and cathode will reach to 34% and 58% after 4230 h, respectively. OCV operation: under this condition, current density of the PEMFC is zero. The corresponding cell voltage of OCV operation is ~0.95 V vs. RHE. It usually occurs immediately after the PEMFC starts up or exactly before shutting down. During the fuel cell vehicle lifetime, the OCV operation happens about 100 h [2]. Ferreira et al. performed experimental study to test idling and also OCV operations of the PEMFC [4]; in their study, the FC was run for 2000 h under 0.75 V (corresponding to 0.2 A cm2) and 0.95 V in idling and OCV operations, respectively. Their obtained results indicated that the ECSA decline rate for OCV operation is faster than idling one, thus after 2000 h of operation, the ECSA decreases from 2 1 2 1 2 1 75 m2 g1 Pt to 40 m gPt and from 75 m gPt to 15 m gPt for idling and OCV operation, respectively. Cyclic load: it is the most prevalent and probable kind of transient operation for the FC vehicles. Different kinds of Accelerated Degradation Tests (ADT) have been defined to investigate the CL durability and lifetime during cyclic load [5,6]. Many researches were implemented about this issue [4,6e19]. Bi et al. investigated the influence of temperature and Relative Humidity (RH) of the reactant on the degradation rate and ECSA decline [7,8]. In their experiments, the ECSA is measured after some cyclic load; they revealed that an increase in the temperature will result a higher rate of ECSA decline [7]. Furthermore, their results showed that the relative humidity also plays an important role on the Pt degradation mechanism; an increase of the relative humidity causes enrichment of the ECSA decline, and consequently, the percentage of Voltage Degradation Rate (VDR%) [8]. Ferreira et al. also studied effects of cyclic load, (cyclic voltage from 0.6 V to 1 V) on the TKK 2 1 catalyst layer [4]. The ECSA decreased from 63 m2 g1 Pt to 23 m gPt after 10000 cycles (up to 63% reduction). Their obtained X-Ray Diffraction (XRD) results showed that the volume average of Pt particles radius increased from 2.3 nm for pristine crystal to 10.5 nm for degraded one. Experiments on the effects of different types of cyclic load protocol were performed [6,13,20], e.g. Paik et al. and Ohma et al. applied triangle and square cyclic voltage waves to the FC and showed that for the triangle/square protocol with the same minimum, maximum voltage and time period, the ECSA decline is approximately the same, and will be closer for the larger time period [6,13]. Ohma et al. also investigated the impact of amplitude and average of cyclic load protocol on the ECSA decline [6]. Their results show that the higher amplitude of cyclic load the higher rate of ECSA decline. Some different techniques were proposed for reducing the ECSA decline rate such as: synthesis of catalyst layer with some other noble metal [15,16,19], controlling operating parameters [7,8,17] and using Nano-Structure Thin Film (NSTF) [9]. Debe et al. compared two different catalyst layers (Pt/C and NSTF) during the cyclic load protocol of 0.6Ve1.2 V vs. RHE [9]. The normalized ECSA decreased to 0.7 and 0.1 for NSTF and Pt/C after about 2000 cycles, respectively. Their obtained results indicates that NSTF is more resistant against the ECSA decline caused by the cyclic loads respect to the Pt/C.

Cyclic Start-Stop: it is the frequent startup and shutdown of PEMFC. Frequent start-stop cycling causes an instantaneous condition which lead to fill the anode with H2 and air at once, and thus increasing the local interfacial voltage reaches as high as 1.5 V which is very dangerous for catalyst layer lifetime. However, the PEMFC of vehicle is exposed to this condition approximately 38500 cycles in its total lifetime [2]. There are some simple modeling of ECSA degradation for cathode catalyst layer under cyclic load which have been presented by Bi et al. and Debe et al. in their experimental investigations [7e9]. In that simple model a relationship between the ECSA and the number of cycles, N, was presented which predicts the ECSA for a specific cyclic load (not any given arbitrary cyclic protocol) at definite relative humidity. Their model also did not consider any procedure to use the predicted ECSA after N cycles to forecast the cell performance under any cyclic load. In this study, an experimental based algorithm is presented to predict the cell performance for any given arbitrary cyclic load. The algorithm consists of two models, namely Models 1 and 2 as shown in the I-O diagram of Fig. 1. The Model 1 calculates the ElectroChemical Surface Area (ECSA) and agglomerate size (e.g. agglomerate radius, rt,agg) under any given cyclic load at different relative humidity and cell temperature. This task is performed based on cathode catalyst layer degradation and the Ostwald ripening, which relies on numerous experimental data obtained from the literature. Model 2 computes the catalyst layer performance with fixed structural parameters, say, ECSA and agglomerate size, rt,agg [21]. Combination of Models 1 and 2 predicts the catalyst layer performance under an arbitrary cyclic load at different relative humidity and cell temperature. The accuracy of the algorithm (combination of Models 1 and 2) is assessed versus extensive experimental data (over 25 cases), as reported in Sections 5.1 to 5.3. A sensitivity study on the percentage of Voltage Degradation Rate (VDR %) was performed to investigate the effects of operating parameters (such as temperature, relative humidity, pressure, minimum and maximum voltage of the cyclic load) at various numbers of cycles with rank 1 for the most influential parameter. 2. Descriptions of models 1 and 2 Model 1 calculates the Electro-Chemical Surface Area (ECSA) and agglomerate size (i.e. agglomerate radius, rt,agg) under a given number of cycles N, the type (also called as protocol) of the cyclic load, and the cell operating conditions. Numerous experimental observation attests that the effects of these parameters are

Fig. 1. I-O diagram repressing combinations of Models 1 and 2 to predict the catalyst layer performance for a given arbitrary cyclic load (a) Model 1 and (b) Model 2 based on Refs [21,22].

M. Moein-Jahromi et al. / Journal of Power Sources 359 (2017) 611e625

appeared in ECSA and rt,agg [6e9,13,23]. The parameters ECSA and rt,agg act very important role in kinetics of electrochemical reaction within the catalyst layer (CL performance). According to our previous investigation [21], the oxygen reaction rate and subsequently the CL performance depend on these two parameters. In a nutshell, the catalyst layer performance may be predicted for any given structural parameters, e.g. ECSA and rt;agg by means of our previous simulation, see in Fig. 1(b), which is repeated here for convenience by name of Model 2. It forms a Boundary Value Problem (BVP) consisting three ODEs and three Boundary Conditions (BCs) at the sides. They are reported in the following equations in which the two mentioned parameters (ECSA and rt,agg) are highlighted: ODEs:

(1)

ODE­2 : O2 concentration gradient

ODE­3 : activation over  potential

VCO2 ¼

Vhact

dCO2 i  Id ¼ eff ; dz 4FDO2

(2)

where kC is: (4) The subscript t in rt;agg represents the total time that catalyst has been under the cyclic load. For example, if a catalyst layer is under N cycles with a fixed repeating unit Dt for a cycle, t becomes t ¼ N  Dt: So, a given N infers to a corresponding fixed t. Eqs. (1) and (4) show the direct effect of ECSA and rt,agg on cell performance in Model 2. BCs:

catalyst layer=membrane interface


dCO2

¼ 0; dz
CL=membrane (5)

BC­2 : ¼ 0; i catalyst layer=GDL interface CL=GDL

BC­3 : catalyst layer=GDL interface

 
PO2 CO2
CL=GDL ¼ ð1  sÞ: RT |fflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflffl} CO

eff

ionomer, Id is the total current density of the cell, DO2 is the total effective diffusion coefficient of oxygen into the catalyst layer which is described in Appendix A. εCL is the catalyst layer porosity, keff is bulk effective protonic conductivity, seff is bulk effective electronic conductivity, COref2 is reference oxygen concentration, iref 0 is reference exchange current density, hact is activation over potential, F is the Faraday constant, aa and ac are cathodic and anodic transfer coefficient, respectively. Some of these parameters are reported in Table 1 and the others are explained with full details in Ref. [21]. On the other side, Model 1 calculates the ECSA and rt,agg during the cyclic load after any given number of cycles N, the details are described in the following:

2.1. Electro-Chemical Surface Area (ESCA)

SN ¼

dh i iI ¼ act ¼ eff þ eff d : dz k s

VCO2 ¼

the catalyst layer, dagg is the thin ionomer film surrounding each agglomerate, mPt is platinum loading of catalyst layer, LCL is the catalyst layer thickness, DO2 ;i is oxygen diffusion coefficient in

The ðECSAÞN is normalized by its primary value ðECSAÞN¼0 using:

(3)

BC­1 :

613

(6)

2 ;wv

ðECSAÞN : ðECSAÞN¼0

So, For a pristine catalyst it is SN¼0 ¼ 1. SN , is considered to be a function of the number of cycles N, where the rate of ECSA decline could be shown by dSN =dN. Assuming a first order ECSA decline rate, dSN =dN would be linearly proportional to its SN [9,12], so:

dSN dS fSN 0 N ¼ KSN ; dN dN

(9)

where K is the proportionality constant. The minus sign in Eq. (9) corresponds to the decline of ECSA. Integration from Eq. (9) within the limits N2 ½0; ∞Þ provides SN 2½0; 1, which is a linear plot between lnðSN Þ and N. However, several experimental data do not fit to the noted linear plot [9]. In fact, they support the existence of an endurance limit for the number cycles Ne , beyond which ECSA remains stable at a minimum value [7,9,24]. That is Smin ¼ SN > Ne : Then:

Smin ¼

ðECSAÞmin : ðECSAÞN¼0

(10)

That means there exists a maximum for the normalized decayed ECSA that corresponds to 1  Smin . Hence the integrating limit is taken as S2 ½1  Smin ; SN  Smin . So one can write:

  PO2 þ s: HO2 ; W |fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl} CO

In Eqs. (1)e(4), PO2 is oxygen partial pressure in cathode, CO2 is oxygen concentration in catalyst layer, HO2 ;i is the Henry's constant of oxygen into the ionomer, Er is agglomerate effectiveness factor which is explained in Ref. [21], rt;agg is agglomerate particle radius, aagg is the total surface of all agglomerate particles (usable to dissolve oxygen into the agglomerate particles) per unit volume of

(8)

SðNÞS Z min

1Smin

;

(7)

2 ;lw

  ZN dS SðNÞ  Smin ¼  KdN0ln ¼ KN; S 1  Smin 0

(11)

SðNÞ ¼ Smin þ ð1  Smin ÞexpðKNÞ: Lots of experimental studies performed on the effect of cyclic

614

M. Moein-Jahromi et al. / Journal of Power Sources 359 (2017) 611e625

Table 1 Operating, structural and kinetic parameters (base case condition). Parameters

Quantity

Value/Unit

Reference

T P RH Id XO2 XH2 LCL

Temperature, Pressure, Relative Humidity, Total current density, Mole fraction of oxygen in cathode side, Mole fraction of hydrogen in anode side, CL thickness, CL tortuosity, Initial ECSA of pristine CL, GDL porosity, Volume fraction of GDL in CL, Ohmic resistance, Mass loading of platinum per unit area of the cathode, Mass loading of carbon per unit area of the cathode, Platinum density, Carbon density, Reference oxygen concentration,

60  C 1 atm 100% 8000 A m2 21% 100% 50 mm 1:5  0:5  εCL 72000 m2 kg1 0.78 0.1 7.2  106 U m2 0.005 kg m2 0.045 kg m2 21400 kg m3 1800 kg m3 1.2 mol m3

[7] [7] [7] e [7] [7] [36] [7] [35] Fitted [7] [7] [7] e e [21]

Cathodic transfer coefficient, Anodic transfer coefficient, Bulk protonic conductivity, Bulk electronic conductivity, Initial radius of agglomerate, Thickness of agglomerate, Initial radius of Pt particles, Reference exchange current density

1.0 0.5 17 U1 m1 72700 U1 m1 0.3 mm 40 nm 3.3 nm 10ð0:0341T16:96Þ A m2

[21] [21] [21] [21] Fitted Fitted [7] [21]

Henry's constant of oxygen in ionomer, Henry's constant of oxygen in water, Oxygen diffusion coefficient in ionomer phase, Pt diffusion coefficient in ionomer phase, Interfacial energy surface density of Pt, Pt molar volume, Proportionality constant under standard condition,

0.3125 atm m3 mol1 1:33 expð666=TÞatm m3 mol1 1:4276  1011  T  4:2185  109 m2 s1 1  1010 m2 s1 2.37 J m2 9.09  106 m3 mol1 4.537  104

[21] [21] [21] [37] [28] e [9]

t

ðECSAÞN¼0 εGDL LGDL Rohmic mPt mC

rPt rC

ref

CO 2 aC aa

k s

r0;agg

dagg R0;Pt iref 0 HO2 ;i HO2 ;w DO2 ;i DPt;i

gPt UPt KSC

load on the ECSA decline, showed the consistency of Eq. (11) with experimental data [7,9,10,24]. It is worthwhile to mention that there is a very similar phenomenon in solid mechanic science in which the strength of carbon steel under fatigue load (a cyclic stress load) decreases in a way that it is linearly proportional to the logarithm of cycle number, see Fig. 2(a) [25,26]. If the cyclic stress load on carbon steel is small enough, its fatigue strength reduces to a minimum value after a certain number of cycles beyond which fatigue strength remains stable (not broken area in Fig. 2(a) which is called endurance limit) [26]. By comparing the loss in fatigue strength, Fig. 2(a), and the ECSA decline of the catalyst layer, Fig. 2(b), one can conclude the analogy between them. The ECSA/ fatigue strength reduces logarithmically with increasing number of cycles and finally reaches to its minimum value after certain number of cycles. In other word, it seems that after cyclic load imposition, the catalyst layer/carbon steel may transmute into a new material with lower ECSA/fatigue strength. By using this analogy, it is possible to adopt the same solution procedure as in fatigue analysis of carbon steel to define standard conditions and then modify other different operating conditions by considering some correction factors. The standard condition is defined as the following:
Temperature: TSC ¼ 60o Cð333KÞ;
Cathode Relative Humidity: RHSC ¼ 100%;
Cyclic load protocol is a cyclic square-shaped voltage with maximum and minimum voltage of Vmax ¼ 1:2 V and Vmin ¼ 0:6 V: The repeating unit for a cycle is considered as Dt ¼ 60 s: Therefore, the proportionality constant in Eq. (11), K, is considered to be revised by three correction factors which corresponding

to the main effective operating parameters affected on the ECSA decline rate such as: temperature [7,9], relative humidity [8] and cyclic load protocol [6,13,20,23]. If the operating conditions differs from the standard condition, K will be modified to:

K ¼ KSC  KT  KRH  KV ;

(12)

in which KSC is proportionality constant under standard condition. KT , KRH and KV are: KT : Temperature correction factor which is equal to KT ¼ 1 if T ¼ TSC , KRH : Relative humidity correction factor which is equal to KRH ¼ 1 if RH ¼ RHSC , KV : Cyclic load protocol correction factor which is equal to KV ¼ 1 for the standard cyclic load protocol. 2.1.1. Temperature correction factor If temperature is the only violating parameter from the standard condition, then K would be:

K ¼ KSC  KT :

(13)

K is a rate constant having the typical Arrhenius dependence on the temperature, similar to lots of other chemical equations, so it is:

  Ea Ea 0lnðKÞ ¼ lnðGÞ  : K ¼ G exp RT RT

(14)

Eq. (14) is an empirical relation, in which G is the preexponential factor, Ea is the activation energy for ECSA decline reaction and R is the universal gas constant. Debe et al. experimentally obtained the proportionality constant at standard condition, KSC, it is reported in Table 1 [9]. They also

M. Moein-Jahromi et al. / Journal of Power Sources 359 (2017) 611e625

615

obtained K for different temperature under standard cyclic load protocol and relative humidity. Curve fitting of their experimental results between K and T regarding to Eq. (14) yields that K depends on T according to:

  2720 : K ¼ exp 0:47  T

(15)

Substitution of Eq. (15) into Eq. (13) yields the temperature correction factor, it is:

KT ¼

K ¼ KSC

  exp 0:47  2720 T KSC

:

(16)

2.1.2. Cyclic load protocol correction factor With considering cyclic load protocol deviation as the only parameter which differs from the standard condition, Eq. (12) reduces to:

K ¼ KSC  KV :

(17)

There are three important factors influence on KV ; they are (i) amplitude, Vamp ; (ii) average, Vave and (iii) the period time, Dt; of the cyclic load protocol. The mechanism of platinum particles (reaction sites) solution into ionomer phase of the catalyst layer and also the Ostwald ripening of platinum particles are intensified by increasing Vamp : The experimental data demonstrates that lnðSÞ is inversely dependent to Vamp [6]. On the other hand, Eq. (11) shows that lnðSÞ is also inversely proportional to K. Thus one can conclude that:

lnðSÞ inverselyfVamp lnðSÞ inverselyfK

 /KfVamp :

(18)

Wang et al. [27], Darling et al. [28] and Shao-Horn et al. [29] showed that the platinum solute concentration in the ionomer phase is exponentially proportional to Vave ; it is described in SubSection 3.1 in detail. However, the more platinum solute concentration in the ionomer phase, the lower value of normalized ECSA (means that: lnðSÞ inversely f C∞;Pt ). Therefore, K, would be proportional to expðVave Þ: The mentioned description is summarized to:

9 C∞;Pt fexpðVave Þ = lnðSÞ inverselyfC∞;Pt /KfexpðVave Þ: ; lnðSÞ inverselyfK

(19)

C∞;Pt f t ¼ N Dt !

C∞;Pt fDt

K

(21)

Considering the defined standard conditions and Eq. (21), it can be expressed that:

K3

(22)

Therefore, regarding to Eqs. (17) and (22) and the standard conditions, KV is: K3

zfflffl}|fflffl{ zfflfflfflfflffl}|fflfflfflfflffl{ K2 Vamp ðVÞ zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{ DtðsÞ KV ¼ $expðVave ðVÞ  0:9Þ$ ; 60 0:6

(23)

which splits into three parts:

KV ¼ K1  K2  K3 ;

9 > > = lnðSÞ inverselyfC∞;Pt lnðSÞ inverselyfK /KfDt: > > ;

Summing up Eqs. (18e20) results that:

KfVamp  expðVave Þ  Dt:

K

2 1 zfflfflfflfflfflffl}|fflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflffl{ zfflfflffl}|fflfflffl{ Vamp K expðVave Þ Dt   $ $ ¼ : KSC Vamp SC exp ðVave ÞSC ðDtÞSC

K1

Wang et al. showed that the solubility of platinum particles in the ionomer phase is linearly proportional to the working time duration, t. Therefore, for the predefined number of cycles, N, it would be expressed that Platinum solute concentration in the ionomer is proportional to time period, Dt [27], so it can be concluded that:

Fixed N

Fig. 2. Analogy between (a) fatigue phenomenon of carbon steel (it is reprinted with permission from Taylor and Francis Group LLC book [26]) and (b) ECSA decline during cyclic load according to the experiments of Ref. [7].

(24)

(20)

in which K1, K2 and K3 are the correction factor for deviation of Vamp ; Vave and Dt from the standard condition, respectively. 2.1.3. Relative humidity correction factor Bi et al. investigated the influence of relative humidity on the

616

M. Moein-Jahromi et al. / Journal of Power Sources 359 (2017) 611e625

ECSA decline rate [8]. They experimentally obtained K at different relative humidity. Their results indicate that K is approximately proportional to the relative humidity, i.e. by decreasing relative humidity from 100% to 50%, K will be reduced by around half [8], as a result the relative humidity correction factor would be:

KfRH0KRH ¼

K RH RH% : ¼ ¼ KSC RHSC 100%

(25)

ripening and form a bigger one with radius of Rt;Pt ; the interfacial energy will be:

 ERt ;Pt ¼ gPt 4pR2t;Pt :

(27)

From mass conservation, one can conclude:

massPt ¼ rPt

    1 4 3 4 pRt;Pt ¼ 2rPt pR30;Pt 0Rt;Pt ¼ 23 R0;Pt : 3 3 (28)

2.2. Agglomerate radius during the Ostwald ripening of catalyst layer

Replacing Eq. (28) into Eq. (27) and comparing the result with Eq. (26) yields that:

During the Ostwald ripening, platinum particles have tendency to dissolve into the ionomer phase by departing two electrons and forming Pt2þ (Pt / Pt2þþ2e, dissolution). Pt2þ could be dissolved in ionomer more easily. Then, in the next reaction, the solute Pt2þ ions in the ionomer phase aggregate together and precipitate into the carbon support. It forms bigger platinum particles (Pt2þþ2e/ bigger Pt, precipitation) [30]. Different stages of this phenomenon, which is called Ostwald ripening, are shown in Fig. 3(a). Small particles shrink and the other big particles grow, i.e. the large particles swallow the small ones. In this effect, the mean radius of platinum particles grow but their numbers decrease. Aggregating small platinum particles into larger one leads to minimize their interfacial energy, it is the main motivation for Ostwald ripening, e.g. for two different Pt particles with initial radius R0;Pt and the interfacial energy surface density gPt ; the total interfacial energy will be:





ER0 ;Pt ¼ 2gPt 4pR20;Pt :

(26)

(29)

The Ostwald ripening model simulates dissolution of the dilute system of dispersed particles in a media or matrix [31,32]. So, the dissolution of Pt particles in ionomer phase is modeled as a dilute system of Pt particles in the ionomer media by using Ostwald ripening. Hence, the solute concentration of Pt particle with radius Rt;Pt in ionomer phase, may be described by well-known GibbsThomson relation, it is [31,32]:

  n ; CRt;Pt ¼ C∞;Pt exp Rt;Pt

where, C∞;Pt is the solute concentration in ionomer phase for a flat Pt particle (infinite radius). n is a Capillary length defined as [32]:

n¼

2gPt UPt ; RT

(31)

(ii)

(iii)

(a)

Ostwald Ripening Timeline Ionomer Media Pt Particle

aagg

O2

O2

O2

O2 O2 O2

O2

Carbon Support Ionomer Phase Pt Particle

(b)

Penetration Surface of Agglomerate against O2 Electrochemical Surface Area (ECSA)

ECSA

(30)

in which UPt is Platinum molar volume. Eq. (30) could be approximated to:

If these two particles are joined together by the Ostwald

(i)

ERt ;Pt < ER0 ;Pt :

O2 Flux into the Agglomerate

Fig. 3. (a) Ostwald ripening procedure of Pt particles in ionomer during the time, (b) oxygen penetration flux into an agglomerate to achieve to the ECSA.

M. Moein-Jahromi et al. / Journal of Power Sources 359 (2017) 611e625

 CRt;Pt zC∞;Pt 1 þ



n Rt;Pt

:

(32)

Eqs. (39) and (41) show that the Pt particles with smaller/bigger radius (Rt;Pt ) than the Rt;Pt;critical shrink/coarsen, as it has been shown in Fig. 3(a) from (i) to (iii), in other words:

if Rt;Pt < Rt;Pt;critical ;

The Pt concentration conservation is:

V2 CPt ¼ 0:

(33)

if Rt;Pt > Rt;Pt;critical ;

Eq. (33) is solved in spherical coordinate subjecting to two following boundary conditions [32]:  (a) Cr¼Rt;Pt



n ; Gibbs-Thomson equation for solute ¼ C∞;Pt 1 þ Rt;Pt

concentration of Pt particle with radius Rt,Pt, see Eq. (32). (b) Cr/∞ ¼ C Pt : C Pt is the mean concentration of Pt particle in the bulk ionomer far from the reaction sites. The solute concentration in ionomer at distance r from Pt particle with radius Rt;Pt is obtained by solving Eq. (33), it is:

iR h C∞;Pt n t;Pt þ : CPt ðrÞ ¼ C Pt  C Pt  C∞;Pt r r

(34)

On the other hand, the platinum mass balance in the catalyst layer is:



d ðnPt Þ ¼  4pR2t;Pt N_ Pt ; dt

(35)

where, N_ Pt is the Pt molar flux into the ionomer. The minus sign in Eq. (35) corresponds to the Pt particles separation from the catalyst layer. nPt which is the platinum moles inside the catalyst layer would be expressed as:

nPt ¼

rpt 4p

m R3 : ¼ MPt MPt 3 t;Pt

(36)

MPt and rPt are the molecular weight and density of Pt, respectively. Pt molar flux is expressed according to the Fick's law, it is:



N_ Pt ¼  DPt;i VCPt r¼R :

(37)

t;Pt

(38)

in which, DPt;i is Pt diffusion coefficient in ionomer phase. Time variation of the Pt particle radius is obtained by replacing Eq. (34) in Eq. (38), it is [32]:

 ! dRt;Pt UPt DPt;i C Pt  C∞;Pt C∞;Pt ¼ n : Rt;Pt  dt R2t;Pt C Pt  C∞;Pt

(39)

Lifshitz, Slyozov and Wagner solved Eq. (39) with considering the assumption that the distance between the Pt particles compared to their radius is large enough in order to have no interaction on each other. Also, the volume fraction of the dispersed Pt in the ionomer phase is very small (dilute dissolution) [31e33]. Their analytic and asymptotic solution of Eq. (39) is:

R3t;Pt  R30;Pt ¼

8gPt C∞;Pt U2Pt DPt;i 9RT

t ¼ kd t;

kd ≡

8gPt C∞;Pt U2Pt DPt;i 9RT

:

(40) Critical radius of Pt particle is obtained by solving dRt;Pt =dt ¼ 0, so it is [32]:

Rt;Pt;critical ¼

C∞;Pt C Pt  C∞;Pt

n:

dRt;Pt < 0; dt dRt;Pt > 0; then dt then

so Pt Particles shrink; so Pt Particles coarsen: (42)

Coarsening of the Pt particles according to Eq. (40), and the Pt solute concentration, Eq. (32), indicate that the driving moment for the Ostwald ripening is D ¼ CRt;Pt  C∞;Pt [30]. As the Pt particles grow during Ostwald ripening (Rt;Pt increases), CRt;Pt tends to C∞;Pt (meaning D/0), see Eq. (32), which leads to reduce the rate of the Ostwald ripening. Consequently, as long as Ds0, the Ostwald ripening will occur to coarsen the catalyst particles. When the PEMFC is under cyclic load, D is not zero (Ds0), since C∞;Pt depends on the average voltage of cyclic protocol, see Eq. (42), hence it will not be constant during the potential cycling. Thus, the driving moment of the Pt ripening is more severe in cyclic load respect to the idle load in which D tends to zero with time goes on. As illustrated in Fig. 3(b), the total surface of all agglomerate particles (usable to dissolve oxygen into the agglomerate particles) per unit volume of the catalyst layer in terms of [m2 m3].  unit  volume of the catalyst aagg ; is proportional to the ECSA per Pt layer in terms of [m2 m3] which is m  ECSA; therefore [34]: LCL



aagg f

mPt LCL



 ECSA:

(43)

The surface area is proportional to square of its corresponding radius, meaning areaf R2 , so for a fixed value of the platinum loading, mPt ; and thickness of the catalyst layer, LCL ; the similarity Eq. (43) is summarized to:

rt;agg r0;agg

!2

 ¼

Rt;Pt R0;Pt

2

Rt;Pt rt;agg 0 ¼ : r0;agg R0;Pt

(44)

The mean radius of the agglomerate particles at any time, rt;agg ; is obtained by combining Eqs. (40) and (44), it is:

Substitution of Eqs. (36) and (37) into Eq. (35) yields:

   d rPt 4p 3 dC Rt;Pt ¼ 4pR2t;Pt DPt;i Pt dt MPt 3 dr r¼Rt;Pt

617

(41)

rt;agg ¼ r0;agg 

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k t 3 1 þ 3d : R0;Pt

(45)

3. Model input parameters The present model which predicts the catalyst layer performance during cyclic load, depends on many operating/structural parameters which considerably influence on the results. Therefore, many efforts has been done to precisely calculate and collect them from the literature. Some of them which are used in Model 2 are explained in Appendix A and also are reported in Table 1 and the others are described in the following [7,9,21,28,35e37]: 3.1. Solute concentration of flat Pt particle There has not proposed a comprehensive theory about the solute concentration of Pt in the ionomer yet [29]. It seems that the assumption of thermodynamic/electrochemical equilibrium during the Pt dissolution and precipitation in the ionomer is not enough accurate. There are different studies which investigate on Pt solute concentration in different aqueous acidic solutions [4,28,29,38]. They are good pattern to put forward a relation for solute concentration of flat Pt in the ionomer. Ferreira et al. measured the Pt

618

M. Moein-Jahromi et al. / Journal of Power Sources 359 (2017) 611e625

solute concentration in 0.5 M H2SO4 at 80  C [4]. They have shown that the Pt2þ solubility increases with the average voltage of cyclic load. Pourbaix also investigated this effect and indicated that in low potential, the solubility of platinum is negligible [39], however it cannot be neglected if Vave > 0:85 V: Platinum solubility in 25  C was reported in his work against the average voltage of cyclic load. The Pt dissolution in ionomer phase occurs by decomposition of Pt atoms to Pt2þ ions according to Eq. (46). Since, Pt2þ could easily dissolve into the ionomer [28,29].

PtðsÞ/Pt 2þ ðaqÞ þ 2e :

(46)

Some different factors play important role in this issue such as the average voltage of cyclic load and chemical potential of Pt2þ ions [29]. Therefore, Shao-Horn et al. presented an experimental/ analytical model to calculate the solute concentration of flat Pt particle in ionomer phase according to [29]:

    h i 2F Vave  Erev;Pt 4g U þ Pt Pt : C∞;Pt mol m3 ¼ exp Rt;Pt RT RT

(47)

Vave and Erev;Pt are the average voltage of cyclic load protocol, and reversible voltage of Pt decomposition reaction, Eq. (46). ShaoHorn et al. described the calculation of Erev;Pt in detail, it is [29]:

Erev;Pt ðTÞ ¼ 1:188 þ 3:2376  104 ðT  298Þ; T½K:

(48)

step, the degraded values of ðECSAÞN and the agglomerates radius, rt;agg ; after N cycles of the protocol, must be obtained by using Eqs. (11) and (45), respectively (Model 1). In the final step, the ODEs (Eqs. (1)e(3)) are solved by the new value of these parameters (ðECSAÞN and rt;agg ) in order to obtain the degraded catalyst layer performance after N cycles. Shooting method with forth order Runge-Kutta procedure is used to solve the ODEs by using bvp4c function of MATLAB software package. The catalyst layer is divided by 500 nodes which seems enough for independent grid solution. The operating voltage of the cell, VCell ; is obtained by subtracting the losses from the thermodynamic reversible voltage, Vrev ; it is:

VCell ¼ Vrev  hact  Rohmic Id :

(51)

Rohmic is the ohmic loss in term of [U.m2]. It should be pointed out that the ohmic loss of MEA with Pt/C particles, is roughly constant during the cyclic load especially at low cell temperatures [7,40]. In the current study, the used values of Rohmic are extracted in any temperature and cycle number from Bi and fuller's experimental data [7]. hact is attained from the output of Model 2. Vrev is obtained from thermodynamic reaction of producing H2O from O2 and H2 [21], it is:

Vrev ¼ 1:229  0:85  103 ðT  298:15Þ þ 4:31      105 T ln xH2 P þ 0:5ln xO2 P :

(52)

It should be mentioned that, using Eqs. (47) and (40) to evaluate Pt particle radius during the Ostwald ripening leads to form a set of recursive formulation for Rt;Pt : It is solved by trial and error procedure.

T and P, are the operating temperature and pressure in terms of [K] and ½atm, respectively. xH2 and xO2 are hydrogen and oxygen mole fraction at the inlet of anode and cathode, they are:

3.2. Minimum normalized ECSA

xH2 ¼

The ECSA is proportional to APt : APt is defined as the ratio of the outer surface of all Pt particles per its unit mass: APt ¼

4pR2r;Pt . rPt ð4=3pÞR3r;Pt

APt is inversely proportional to the Pt radius, Rt;Pt ;[7]. Therefore, it

sat P  RH  PH 2O

P

;

3:3 z0:1: 30:5

:

(53)

 105 ðT  273:15Þ2 þ 1:4454  107 ðT  273:15Þ3 ;

T½K:

(49)

The DOE roadmap target by 2020 for PEMFC catalyst layer under cyclic load protocol is at least 5000 h, time life durability [1]. If the initial radius of Pt is supposed to be R0;Pt ¼ 3:3 nm[7], according to Eq. (40) after about 5000 h cyclic load under standard conditions, the Pt radius will be increased to Rt;Pt ¼ 30:5 nm; thus, the minimum ECSA is considered as:

Smin ¼

4:76P

 sat log10 PH ½atm ¼ 2:1794 þ 0:02953ðT  273:15Þ  9:1837 2O

ECSA, Smin ; is estimated as:

Smin

sat P  RH  PH 2O

sat is saturation pressure of water at the cell operating temperaPH 2O ture in term of ½atm, it is [41]:

cloud be stated that, ECSAfR1 t;Pt : So, the minimum normalized

R0;Pt ðECSAÞmin  ¼ ¼ : ðECSAÞN¼0 Rt;Pt min

xO2 ¼

(50)

In this study, the value of Smin is considered to be 0.1. Similar values are remarked in other investigations such as Refs. [9,24,40]. 4. Solution procedure The complete solution procedure was shown in Fig. 4. The first step is guessing the initial value of agglomerate radius, r0;agg ; and also extracting the initial value of ðECSAÞN¼0 for the pristine catalyst layer from experiments. In the second step, the performance of the catalyst layer is calculated using Model 2 (Eqs. (1)e(3)) and the values of r0;agg and ðECSAÞN¼0 : If the performance of the pristine catalyst layer does not comply with the result of experiments, the solution procedure should be repeated from the first step by guessing new value for r0;agg : In the next

(54) 5. Results and discussion The results are validated by experimental data in three steps. The first is about the validation of the Model 1 output, ECSA, against a large number of experiments. The second is verification of the growth of Pt particles radius during the Ostwald ripening, and the last is the validation of the catalyst layer performance calculator during cyclic load, Model 2. After validation purpose, the influence of five operating parameters (T, RH, P, Vmin and Vmax) are studied on the degradation rate of the catalyst layer performance and the ECSA decline. All of the parameters which are selected here, are under control of the PEMFC operators. Therefore, the results of this study are very useful to minimize the percentage of the Voltage Degradation Rate (VDR%). VDR% is defined as:

VDR% ¼

N¼0  V N Vcell cell N¼0 Vcell

 100 ¼

1

N Vcell N¼0 Vcell

!

 100:

(55)

The parametric study is performed at current density Id ¼ 8000 A m2 for 1000, 2500, 4000 and 5500 numbers of cycles under the standard condition. At last, a sensitivity analysis is performed to rank the most influential parameters on the VDR%.

M. Moein-Jahromi et al. / Journal of Power Sources 359 (2017) 611e625

619

Fig. 4. Solution procedure flow chart.

5.1. ECSA validation The ECSA calculator during cyclic load under different conditions is validated with a large number of experiments. Fig. 5 shows that the predicted ECSA significantly complies with a lot of experimental data of different operating conditions [4,6e13,16,17,23,24,34,42e53]. The Model 1 is theoretical and empirical based model which predicts the ECSA decline very well during different cyclic load protocol and operating condition. It should be pointed out that the principle of ECSA prediction in present study, is using an analogy between the ECSA degradation and the fatigue of carbon steel in solid mechanic science. As discussed in full detail in the mechanical engineering design books such as Refs. [25,26], the fatigue is very complicated issue exactly like complex nature of ECSA degradation. So, it seems that there are maybe some other reasons which leads to fatigue phenomenon in carbon steel [25,26], and ECSA decline in CL of fuel cell through the transient operation especially

cyclic load [29]. As a result, it is natural that the algorithm of this paper which is based on this analogy and does not consider probable unknown reasons, predicts the ECSA degradation with an acceptable approximation. 5.2. Ostwald ripening of Pt particles validation The Ostwald ripening of Pt particles during cyclic load protocol, Eq. (40), is validated with the experiments of Bi and Fuller [7]. Table 2 indicates a good agreement between their X-Ray Diffraction (XRD) test results and the predicted radius of Pt particle after 7000 cycles under the protocol Vmax ¼ 1:2 V; Vmin ¼ 0:87 V, Dt ¼ 60s at different operating temperatures. It is worthwhile to notice that the tendency of Pt particle to grow up through Ostwald ripening depends on its radius, see Eqs (40) and (47). In fact, the smaller Pt particle, the higher value of Cf,Pt and consequently higher value of kd. On the other side, in a real catalyst layer, Pt particles are spread with a distribution of different radius inside the catalyst. So, during

620

M. Moein-Jahromi et al. / Journal of Power Sources 359 (2017) 611e625

Fig. 5. Verification of predicted normalized ECSA with experimental data of different Refs. [4,6e13,16,17,23,24,34,42e53].

an experimental test, the small Pt particles have a greater growth rate than the larger ones, which leads to create a new distribution of Pt particles with a large variance in radius which result to different value of Ostwald ripening rate, kd [7]. Therefore, considering an average value of radius for all of Pt particles to calculate the Ostwald ripening rate, kd, leads to make some error in the estimation. Considering an average value of cyclic load protocol (Vave) instead of the instant voltage value may add some error, too. So, the presented model of Ostwald ripening overestimates the experimental value by maximum 5% error. 5.3. Performance validation The cathode catalyst layer performance during cyclic load, output of Model 2, is validated by two different cyclic load protocols such as: Case (i): Vmax ¼ 1:2 V; Vmin ¼ 0:87 V and Dt ¼ 60s, at T ¼ 60o C and T ¼ 80o C[7], Case (ii): Vmax ¼ 1:2 V; Vmin ¼ 0:6 V and Dt ¼ 60s, at T ¼ 75o C[9]. Fig. 6(a) and (b) indicate the cathode catalyst layer performance of case (i) for different number of cycles at T ¼ 60 c and 80  C,

respectively. There is a good agreement between the results of present study and the experiments of Ref. [7]. Validation case (ii) is also shown in Fig. 6(c). It shows that the predicted catalyst layer performance is tracing the experimental results for both pristine and degraded catalyst layer. Fig. 6 illustrates the cell performance loss due to the ECSA degradation and Pt particles agglomeration after N number of cyclic load. The voltage reduction after N number of cyclic load is more severe by increasing the cell temperature i.e. if the cell current density is considered as 8000 A m2, then the cell voltage reduces more than 14% and 34% after 2500 number of cyclic load at cell temperature of 60  C and 80  C, respectively. So, the voltage reduction enhances with increasing the number of cyclic load, N. The reason is the ECSA degradation/Pt particle agglomeration with increasing number of cyclic load, N. 5.4. Temperature The influence of operating temperature is investigated on both the ECSA decline rate and the VDR%. Fig. 7(a) shows that the logarithm of normalized ECSA, is linearly proportional to the number of cycles. It also illustrates that the ECSA decline for a given number of cycles N, is more severe at higher temperature. As an example, the normalized ECSA losses to SN¼5500 ¼ 0:17 and SN¼5500 ¼ 0:45 after 5500 cycles at T ¼ 80  C and T ¼ 60  C, respectively. The reason

Table 2 Comparing the results of Pt particles growth during the Ostwald ripening after 7000 voltage cycling from Vmax ¼ 1.2V to Vmin ¼ 0.87V and Dt ¼ 60s with the experimental data of XRD test Ref. [7]. T [K]

C∞,Pt [mol m3]

kd [nm3 s1]

t ¼ N  Dt [s]

R0,Pt [nm]

Rt,Pt [nm] Present study

Rt,Pt [nm] Experimental data from XRD test [7]

Error [%]

313 333 353

3.629  108 1.673  107 2.400  107

2.427  104 1.051  103 1.423  103

7000  60 ¼ 420000 420000 420000

3.3 3.3 3.3

5.17 7.81 8.59

4.9 7.6 8.2

5.51 2.76 4.75

M. Moein-Jahromi et al. / Journal of Power Sources 359 (2017) 611e625

621

Fig. 6. Validation of cathode catalyst layer performance during cyclic load protocol Vmax ¼ 1.2V, Vmin ¼ 0.87V and Dt ¼ 60s, with experiments of Ref. [7], at (a) T ¼ 60oC, (b) T ¼ 80oC, and also during cyclic load protocol Vmax ¼ 1.2V, Vmin ¼ 0.6V and Dt ¼ 60s with experiments of Ref. [9], at (c) T ¼ 75  C.

is that, K enhances exponentially with increasing T, see Eq. (16). On the other hand, according to Eq. (11), the higher value of K, the lower value of S. Fig. 6 shows that the temperature has a significant effect on the loss of I-V curve during the cyclic load. Compared to the case of T ¼ 60  C, the loss in I-V curve for T ¼ 80oC is considerable. The VDR% at T ¼ 60  C and current density Id ¼ 8000 A m2, is about 19.12% after 5500 cycles. However, it will be about 34.53% at Id ¼ 8000 A m2 and T ¼ 80oC after only 2500 cycles. Therefore, it can concluded that the VDR% enhances noticeably with increasing the temperature. To summarize, one can conclude that an increase in cell temperature makes a rise in the ECSA degradation rate, K, see Eq. (14). So, the ECSA of the catalyst layer decreases more severely at higher temperature which leads to reduce the cell performance (meaning increasing VDR%). 5.5. Relative humidity Relative humidity is the second parameter which is studied. Fig. 7(b) reveals that the relative humidity is linearly proportional to the ECSA degradation rate, K. It is mentioned earlier in Eq. (25). This result was also concluded by Bi and Fuller [8]. It seems that presence of water in catalyst layer provides more opportunity for Pt ion diffusion into the ionomer phase, so Pt ion diffusion into the ionomer phase is facilitated by increasing the relative humidity which causes to rise the ECSA degradation rate and finally reduce the cell performance [8]. As an example, for N ¼ 5500, the VDR% decreases from 14.27% to 8.57%, and the N =V N¼0 ) increases from 0.86 to 0.915, normalized cell voltage (Vcell cell by reducing the relative humidity from 100% to 50%, see Fig. 7(c)

and (d). Therefore, relative humidity is one of operating parameters which can control the VDR%. 5.6. Pressure The results indicate that operating pressure does not have significant influence on the VDR%. Fig. 7(e) and (f) shows that as the pressure increases from 2 atm to 4 atm, the normalized cell voltage rises from 0.894 to 0.908 (for N ¼ 5500 cycles and Id ¼ 8000 A m2). It corresponds to the VDR% of 10.51% and 9.10%, respectively. It means that the VDR% is approximately constant for different operating pressure. Therefore, pressure is not a key parameter for controlling the VDR%. Other investigations on catalyst layer degradation also concluded that pressure have a minor effect on the VDR% compared to the other operating parameters such as T and RH [8]. 5.7. Minimum voltage of cyclic load protocol The effect of minimum voltage of cyclic load protocol, Vmin, is investigated on the VDR% and normalized cell voltage in the range of Vmin ¼ 0.2V to Vmin ¼ 1V. If the maximum voltage of cyclic load protocol, Vmax, is fixed, the amplitude of cyclic load, Vamp, will be increased by decreasing Vmin, (meaning that Vamp [ ¼ Vmax  Vmin Y). A rise in Vamp corresponds to an increase in K, and consequently, more ECSA decline, see Eq. (18). The reason is that the driving moment of Ostwald ripening which is discussed in Section 2.2, D, is not zero in cyclic load operation and is intensified by increasing the Vamp. However, the value of D tends to zero

622

M. Moein-Jahromi et al. / Journal of Power Sources 359 (2017) 611e625

Fig. 7. Effects of operating parameters such as (a) temperature, (b, c, d) relative humidity, (e, f) pressure, (g, h) minimum voltage of cyclic load and (i, j) maximum voltage of cyclic load on ECSA decline rate, on normalized cell voltage, and VDR % during cyclic load.

under a constant/almost fixed voltage operation. So as a results, an increase in Vamp ends to more dissolution of Pt particles and also more ECSA degradation rate which means a higher value of VDR%. Fig. 7 (g) and (h) show that if Vmin ¼ 0.2V, after N ¼ 5500 cycles, the normalized cell voltage and also the VDR% are 0.728 and 27.12% respectively. However, they are changed to 0.875 and 12.41% if Vmin ¼ 1V.

VDR%, respectively. If Vmin is constant, the Vamp will be increased with enhancing Vmax, thus, the ECSA decline rate is reinforced and leads to CL performance decreasing. The reason is the same as mentioned the Section 5.7. The normalized cell voltage at Id ¼ 8000 A m2 will be reduced from 0.897 to 0.857 after N ¼ 5500 cycles if Vmax rises from 1V to 1.2V. The VDR also increases from 10.27% to 14.27% in correspond.

5.8. Maximum voltage of cyclic load protocol

6. Concluding remarks

Fig. 7(i) and (j) illustrate the effect of the maximum voltage of cyclic load protocol, Vmax, on the normalized cell voltage and the

In the current study, an empirical based model was presented for prediction the ECSA decline and also the growth of Pt particles

M. Moein-Jahromi et al. / Journal of Power Sources 359 (2017) 611e625 Table 3 Ranking of the most influential parameters on the VDR%. Rank Parameter Range studied 1 2 3 4 5

T [oC] Vmin [V] RH [%] Vmax [V] P ½atm

60 to 80 0.2 to 1 50 to 100 1 to 1.2 2 to 4

Optimum value 60 1 50 1 4

Sensitivity (jDVDRj, %)

j34.53e14.27j ¼ 20.26, after 2500 cycles j27.12e12.41j ¼ 14.71, after 5500 cycles j14.27e8.57 j ¼ 5.70, after 5500 cycles j14.27e10.27j ¼ 4.00, after 5500 cycles j10.51e9.10 j ¼ 1.41, after 5500 cycles

radius under any given cyclic load protocol. The outputs of this modeling were used in an agglomerate model of catalyst layer to forecast the CL performance loss after some cycles. The results show a very good agreement between the present study and the experimental data for ECSA decline, Pt particle size and catalyst layer performance. A set of parametric studies was performed to determine the influence of each parameter on the VDR%. To rank the most effective parameters which were studied, a sensitivity analysis with VDR% as the objective function was also implemented. The result showed that the VDR% is 14.27% under the standard condition, at Id ¼ 8000 A m2 after N ¼ 5500 cycles. It changes if the operating parameters deviate from their defined standard conditions. The results are reported in Table 3 as well as following:
The most influential parameter on the VDR% is temperature. Increase of the temperature from 60  C to 80  C causes a significant intensification in the VDR% from 14.27% to 34.53% after N ¼ 2500 cycles; it means 20.26% variation in the VDR%, see Table 3. The results is more severe for higher number of cycles such as N ¼ 5500. However, the cathode catalyst layer performance rises with increasing the temperature [21]; therefore, it should be an optimum point in which the two cost function, VDR% and cathode catalyst layer performance, are minimized and maximized, respectively. In conclusion, the operating temperature plays a key role on monitoring the VDR% and the performance which is under control of the PEMFC operators.
The second effective parameter is the minimum voltage of cyclic load protocol. The parametric study has been accomplished to change the Vmin in the range of 0.2Ve1V; the results show that the VDR% decreases with increasing Vmin. The corresponding VDR% is 27.12% and 12.41%, respectively. Therefore, the operator behavior for imposing cyclic load to the PEMFC has direct effect on the catalyst layer degradation.
The relative humidity of the reactant is the third parameter. The result of sensitivity analysis indicates that the VDR% will decrease from 14.27% to 8.57% if relative humidity changes from 100% to 50%.
The next effective parameter is maximum voltage of cyclic protocol, Vmax. The influence of Vmax is investigated in the range of 1 Ve1.2 V. The result shows that the VDR% varies from 10.27% to 14.27% respectively.
The pressure is the last parameter which has negligible effect on the VDR% as mentioned earlier. The VDR% reduces from 10.51% at P ¼ 2 atm to 9.1% at P ¼ 4 atm. Nomenclature aagg Crt ;Pt C∞;Pt

Effective outer surface of agglomerate per volume of catalyst layer (m2 m3) Solute concentration of Pt particle (mol m3) Solute concentration of flat Pt with infinity radius (mol m3)

623

C Pt Di,j DKn

Mean Pt concentration in the bulk ionomer (mol m3) Binary diffusion coefficient (m2 s1) Knudsen diffusion coefficient (m2 s1)

dP Er Erev,Pt F HO2 ;i HO2 ;w Id K LCL LS Li lO2 mPt mC N Nagg nagg

Mean diameter of pores in the porous catalyst layer (m) Agglomerate effectiveness factor Reversible potential of decomposition reaction of Pt (V) Faraday's constant (96485.3365 C mol1) Henry's constant of oxygen in ionomer (atm m3 mol1) Henry's constant of oxygen in water (atm m3 mol1) Total current density (A m2) Proportionality constant of ECSA decline rate, Eq. (11) Thickness of cathode catalyst layer (m) The solid phase volume fraction of catalyst layer The ionomer volume fraction of catalyst layer Mean free path of oxygen inside the catalyst layer (m) Platinum loading (kg m2) Carbon loading (kg m2) Number of cycles Total number of agglomerates Total number of agglomerates per catalyst layer volume (m3)

N_ Pt nPt sat PH 2O PO2 r0,agg rt,agg R Rohmic S Vrev

Pt mole flux (mol m2 s1) Mole of Pt (mol) Saturation pressure of water (atm) Local oxygen pressure (atm) Initial radius of agglomerate for pristine catalyst layer (m) Radius of agglomerate for degraded catalyst layer (m) Universal gas constant (8.314 J mol1 K1) Specific ohmic resistance (U m2) Normalized ECSA Reversible thermodynamic PEM fuel cell voltage (V)

Greek letter aa Anodic transfer coefficient ac Cathodic transfer coefficient dagg Thickness of the electrolyte film covering an agglomerate (m) εagg Ionomer phase volume fraction inside an agglomerate particle εCL Catalyst layer porosity εGDL Gas diffusion layer porosity k Bulk protonic conductivity (U1 m1) gPt Interfacial energy surface density of Pt UPt Pt molar volume, hact Local activation over potential (V) r Density (kg m3) s Bulk electronic conductivity (U1 m1) t Tortuosity of catalyst layer mwater Liquid water viscosity (cP) Subscripts and superscripts amp Amplitude agg Agglomerate ave Average eff Effective i Ionomer max Maximum min Minimum P Pore Pt Platinum ref Reference rev Reversible RH Relative humidity S Solid

624

M. Moein-Jahromi et al. / Journal of Power Sources 359 (2017) 611e625

other word if Knudsen number is in the range of 0:1 < Kn ¼ l=dp < 10, both molecular and Knudsen diffusion will be dominant. So, DO2 ;vapor e will be calculated as [54e56]:

Abbreviation ADT Accelerated Degradation Test CL Catalyst Layer ECSA Electro-Chemical Surface Reaction FC Fuel Cell GDL Gas Diffusion Layer OCV Open Circuit Voltage ODE Ordinary Differential Equation ORR Oxygen Reduction Reaction PEMFC Polymer Exchange Membrane Fuel Cell RH Relative Humidity RHE Reference Hydrogen Electrode VDR Voltage Degradation Reaction

1 DO2 ;vaper

e

¼

1 1 þ : DO2 ;vaper DKn

(A.3)

DO2 ;vapor and DKn are binary diffusion coefficient of oxygen into water vapor and the Knudsen diffusion coefficient. Binary diffusion coefficient of i species into j is [57]:

h

Di;j m2 s1

i

0 1b 1=3  5=12 aB T C  ¼ @qffiffiffiffiffiffiffiffiffiffiffiffiA PCi PCj TCi TCj P T T Ci Cj

1 1 þ Mi Mj

Appendix A. Effective diffusion coefficient of oxygen into the catalyst layer, Deff O2 The effective diffusion coefficient of oxygen into the catalyst eff layer, DO2 , is determined here. The transport routes of the reacting species towards the reaction sites are schematically shown in Fig. A.1.

!1=2 ;

T½K;

(A.4)

P½atm:

The a and b parameters in Eq (A.4) are 3.64  108 and 2.334 for water vapor and oxygen, respectively. Mi, TCi and PCi are molecular weight, critical temperature and pressure of i species. Oxygen Knudsen diffusion coefficient in porous media with pore

Fig. A.1. Reactants diffusion routes inside the catalyst layer to reach the reaction sites.

Fig. A.1 illustrates the cathode catalyst layer partitions in the agglomerate model and the oxygen diffusion routes towards the reaction sites. Oxygen to reach the reaction sites within the catalyst layer has two options [54]: (1) penetrating into liquid water, and (2) diffusion in water vapor. Thus, the overall effective oxygen diffusion in the cathode catalyst layer containing: (1) effective oxygen diffusion into the liquid water, diffusion into water vapor, eff

eff

eff DO2 ;water ;

Deff O2 ;vapor e , eff

DO2 ¼ bs DO2 ;water þ b1s DO2 ;vapor e :

(2) effective oxygen

[54]:

(A.1)

Oxygen diffusion coefficient in liquid water, Deff O2 ;water , is explained in detailed in our former work [21]. The effective oxygen diffusion into water vapor is obtained by Bruggeman's correction eff

DO2 ;vaper

e

¼ ½ð1  sÞεCL t DO2 ;vaper e ;

(A.2)

where, DO2 ;vapor e is equivalent oxygen diffusion in water vapor. If mean free path of oxygen molecule (lO2 ) is in the same order of magnitude with mean diameter of porous media pores (dp ), or in

diameter, dp is [54,56,58]:

h i d  8RT 1=2 p DKn m2 s1 ¼ : 3 pMO2

(A.5)

dP is estimated by hydraulic diameter of the porous media (catalyst layer). The hydraulic diameter could be defined as the quadruple ratio of volume available for oxygen flow, Vflow ; to total wetted surface, Awet ; so it could be express as [36]:

dp ¼ dhyd ¼ 4

Vflow : Awet

(A.6)

The total available volume for the oxygen flow is Vflow ¼ εCL VCL ; and the wetted surface is product of number of agglomerates and its outer surface area which is Nagg  4pðrt;agg þ dagg Þ2 : Therefore, Eq. (A.6) is simplified to:

dp ¼ 4

Nagg

εCL  VCL εCL  2 ¼ 4  2 :  4p rt;agg þ dagg nagg  4p rt;agg þ dagg (A.7)

Nagg is the total number of the all agglomerate particles and nagg is

M. Moein-Jahromi et al. / Journal of Power Sources 359 (2017) 611e625

equal to Nagg =VCL :nagg is calculated in our previous work [21], substitution nagg into Eq. (A.7) yields that [36,56]:

dp ¼

  4 εCL rt;agg þ dagg : 3 1  εCL  Ls;GDL

(A.8)

Appendix B. Supplementary data Supplementary data related to this article can be found at http:// dx.doi.org/10.1016/j.jpowsour.2017.05.102. References [1] U. Drive, US DRIVE Partnership, 2013. [2] J. Garche, C.K. Dyer, P.T. Moseley, Z. Ogumi, D.A. Rand, B. Scrosati, Encyclopedia of Electrochemical Power Sources, Newnes, 2013. [3] M.S. Wilson, F.H. Garzon, K.E. Sickafus, S. Gottesfeld, J. Electrochem. Soc. 140 (1993) 2872e2877. [4] P. Ferreira, Y. Shao-Horn, D. Morgan, R. Makharia, S. Kocha, H. Gasteiger, J. Electrochem. Soc. 152 (2005) A2256eA2271. [5] N. Garland, T. Benjamin, J. Kopasz, ECS Trans. 11 (2007) 923e931. [6] A. Ohma, K. Shinohara, A. Iiyama, T. Yoshida, A. Daimaru, ECS Trans. 41 (2011) 775e784. [7] W. Bi, T.F. Fuller, J. Electrochem. Soc. 155 (2008) B215eB221. [8] W. Bi, Q. Sun, Y. Deng, T.F. Fuller, Electrochimica Acta 54 (2009) 1826e1833. [9] M.K. Debe, A.K. Schmoeckel, G.D. Vernstrom, R. Atanasoski, J. Power Sources 161 (2006) 1002e1011. [10] F. Hiraoka, K. Matsuzawa, S. Mitsushima, Electrocatalysis 4 (2013) 10e16. [11] M.F. Mathias, R. Makharia, H.A. Gasteiger, J.J. Conley, T.J. Fuller, C.J. Gittleman, S.S. Kocha, D.P. Miller, C.K. Mittelsteadt, T. Xie, Electrochem. Soc. Interface 14 (2005) 24e35. [12] B. Merzougui, S. Swathirajan, J. Electrochem. Soc. 153 (2006) A2220eA2226. [13] C. Paik, G. Saloka, G. Graham, Electrochem. Solid-State Lett. 10 (2007) B39eB42. [14] A. Riese, D. Banham, S. Ye, X. Sun, J. Electrochem. Soc. 162 (2015) F783eF788. [15] S.V. Selvaganesh, G. Selvarani, P. Sridhar, S. Pitchumani, A. Shukla, Phys. Chem. Chem. Phys. 13 (2011) 12623e12634. [16] S.V. Selvaganesh, G. Selvarani, P. Sridhar, S. Pitchumani, A. Shukla, J. Electrochem. Soc. 159 (2012) B463eB470. [17] M. Uchimura, S. Sugawara, Y. Suzuki, J. Zhang, S.S. Kocha, ECS Trans. 16 (2008) 225e234. [18] Y. Zhang, S. Chen, Y. Wang, W. Ding, R. Wu, L. Li, X. Qi, Z. Wei, J. Power Sources 273 (2015) 62e69. [19] Z.M. Zhou, Z.G. Shao, X.P. Qin, X.G. Chen, Z.D. Wei, B.L. Yi, Int. J. Hydrogen Energy 35 (2010) 1719e1726. [20] A. Marcu, G. Toth, S. Kundu, L. Colmenares, R. Behm, J. Power Sources 215 (2012) 266e273. [21] M. Moein-Jahromi, M. Kermani, Int. J. Hydrogen Energy 37 (2012) 17954e17966. [22] M. Moein-Jahromi, S. Movahed, M. Kermani, J. Power Sources 277 (2015)

625

180e192. [23] P. Urchaga, S.M. Goli, C.A. Rice, in: Meeting Abstracts, the Electrochemical Society, 2012, p. 1316. [24] V. Avakov, A. Aliev, L. Beketaeva, V. Bogdanovskaya, E. Burkovskii, A. Datskevich, B. Ivanitskii, L. Kazanskii, A. Kapustin, O. Korchagin, Russ. J. Electrochem. 50 (2014) 773e788. [25] J.E. Shigley, C.R. Mischke, R.G. Budynas, Mechanical Engineering Design, McGraw-Hill, 2004. [26] S.R. Schmid, B.J. Hamrock, B.O. Jacobson, Fundamentals of Machine Elements: SI Version, CRC Press, 2014. [27] X. Wang, R. Kumar, D.J. Myers, Electrochem. Solid-State Lett. 9 (2006) A225eA227. [28] R.M. Darling, J.P. Meyers, J. Electrochem. Soc. 150 (2003) A1523eA1527. [29] Y. Shao-Horn, W. Sheng, S. Chen, P. Ferreira, E. Holby, D. Morgan, Top. Catal. 46 (2007) 285e305. [30] A.V. Virkar, Y. Zhou, J. Electrochem. Soc. 154 (2007) B540eB547. [31] A. Baldan, J. Mater. Sci. 37 (2002) 2379e2405. [32] J.H. Yao, K. Elder, H. Guo, M. Grant, Phys. Rev. B 47 (1993) 14110. [33] I.M. Lifshitz, V.V. Slyozov, J. Phys. Chem. Solids 19 (1961) 35e50. [34] S. Jomori, N. Nonoyama, T. Yoshida, J. Power Sources 215 (2012) 18e27. [35] FuelCellStore, Toray Carbon Paper 060 Data Sheet, 2016. http://www. fuelcellstore.com/toray-carbon-paper-060/ (Accessed 29 December 2016). [36] M.N. Abbas, J. Eng. Dev. 15 (2011). [37] S. Burlatsky, M. Gummalla, V. Atrazhev, D. Dmitriev, N. Kuzminyh, N. Erikhman, J. Electrochem. Soc. 158 (2011) B322eB330. [38] P. Bindra, S.J. Clouser, E. Yeager, J. Electrochem. Soc. 126 (1979) 1631e1632. [39] M. Pourbaix, Atlas of Electrochemical Equilibria in Aqueous Solutions, Pergamon Press, Oxford, 1966. [40] F.S. Saleh, E.B. Easton, J. Electrochem. Soc. 159 (2012) B546eB553. [41] T.E. Springer, T. Zawodzinski, S. Gottesfeld, J. Electrochem. Soc. 138 (1991) 2334e2342. [42] S.G. Rinaldo, in: Science: Department of Chemistry, 2013. , M. Oezaslan, P. Strasser, J. Electrochem. Soc. 159 (2011) B24eB33. [43] F. Hasche [44] D. He, S. Mu, M. Pan, Carbon 49 (2011) 82e88. [45] J. Speder, A. Zana, I. Spanos, J.J. Kirkensgaard, K. Mortensen, M. Hanzlik, M. Arenz, J. Power Sources 261 (2014) 14e22. [46] A. Marcu, G. Toth, R. Srivastava, P. Strasser, J. Power Sources 208 (2012) 288e295. €t Ulm, 2014. [47] A. Marcu, in, Universita [48] R. Makharia, S. Kocha, P. Yu, M.A. Sweikart, W. Gu, F. Wagner, H.A. Gasteiger, ECS Trans. 1 (2006) 3e18. [49] Z. Xu, H. Zhang, H. Zhong, Q. Lu, Y. Wang, D. Su, Appl. Catal. B Environ. 111 (2012) 264e270. [50] S. Ye, M. Hall, H. Cao, P. He, ECS Trans. 3 (2006) 657e666. [51] Y.-C. Park, K. Kakinuma, M. Uchida, H. Uchida, M. Watanabe, Electrochimica Acta 123 (2014) 84e92. [52] A.C. Fernandes, V.A. Paganin, E.A. Ticianelli, J. Electroanal. Chem. 648 (2010) 156e162. [53] R.L. Borup, R. Mukundan, ECS Trans. 33 (2010) 17e26. [54] N. Khajeh-Hosseini-Dalasm, M. Fesanghary, K. Fushinobu, K. Okazaki, Electrochimica Acta 60 (2012) 55e65. [55] C. Sattefield, in, MIT Press, 1970. [56] D. Huizenga, D.M. Smith, AIChE J. 32 (1986) 1e6. [57] R.B. Bird, W.E. Stewart, E.N. Lightfoot, New York, (2002). [58] M. Abbasi, J. Evans, I. Abramson, AIChE J. 29 (1983) 617e624.