1. technology and computing

Journal of The Electrochemical Society, 162 (7) E73-E83 (2015)
E73
A Multi-Paradigm Computational Model of Materials
Electrochemical Reactivity for Energy Conversion and Storage
Matias A. Quirogaa,b and Alejandro A. Francoa,b,c,∗,z
a Laboratoire
de R´eactivit´e et Chimie des Solides (LRCS), Universit´e de Picardie Jules Verne and CNRS, UMR 7314,
80039 Amiens Cedex, France
b R´
eseau sur le Stockage Electrochimique de l’Energie (RS2E), FR CNRS 3459, France
c ALISTORE European Research Institute, 80039 Amiens Cedex, France
In this paper we report a new multi-paradigm modeling approach devoted to the investigation of the electrochemical reactivity of
materials in electrodes for energy conversion or storage applications. The approach couples an atomistically-resolved Kinetic Monte
Carlo (KMC) modeling module describing the electrochemical kinetics in an active material, with continuum modeling modules
describing reactants transport at the active material/electrolyte nanoscopic interface (electrochemical double layer region) and along
the mesoscale electrode thickness. The KMC module is developed by extending the so-called Variable Step Size Method (VSSM)
algorithm (called here Electrochemical-VSSM) and constitutes the first VSSM extension reported so far which allows calculating
the electrode potential as function of the imposed current density. The KMC module can be parameterized with activation energies
calculated from Density Functional Theory (DFT), and thanks to the coupling with the transport modules, it describes the materials
reactivity in electrochemical conditions. This approach allows us to study how the surface morphology (e.g. distribution of inactive
sites, size of the active material particle, etc.) impacts the performance of the electrode. As an application example, we report here
a computational investigation of the Oxygen Reduction Reaction (ORR) kinetics in a Pt(111)-based Polymer Electrolyte Membrane
Fuel Cell (PEMFC) cathode.
© The Author(s) 2015. Published by ECS. This is an open access article distributed under the terms of the Creative Commons
Attribution 4.0 License (CC BY, http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse of the work in any
medium, provided the original work is properly cited. [DOI: 10.1149/2.1011506jes] All rights reserved.
Manuscript submitted December 26, 2014; revised manuscript received February 20, 2015. Published April 1, 2015.
Because of the global warming and the fossil fuels depletion, zeroemission electrochemical devices for energy conversion and storage,
such as fuel cells and secondary batteries, are called to play a significant role in the sustainable development of the Humanity. The operation principles of these devices involve complex competitions and
synergies between electrochemical, transport and thermo-mechanical
mechanisms occurring at multiple spatial and temporal scales. Multiscale modeling techniques can help on establishing correlations between their materials chemical and microstructural properties, operation conditions, performance and durability.1–3 Such correlations are
important in order to suggest enhanced materials and cells designs.4–7
Several multiscale modeling approaches have been reported so far.
The most popular one consists in extracting data from a lower-scale
and using it as input for the upper-scale via their parameters. For
example, activation barriers for elementary reaction steps can be extracted from Density Functional Theory (DFT) calculations8 and then
injected into Transition State Theory (TST) Eyring’s expressions to
estimate the parameters which are used to solve Mean Field (MF)
kinetic models simulating the evolution of the reactants, intermediate
reaction species and products.9,10
Alternatively, direct multiparadigm multiscale models consist in
coupling “on-the-fly” models developed in the frame of different
paradigms. For instance, continuum equations describing the transport
phenomena of multiple reactants in a porous electrode may be coupled with atomistically-resolved simulations describing electrochemical reactions among these reactants. Several numerical techniques
nowadays are well established to develop such a type of models,
e.g. coupling Continuum Fluid Dynamics (CFD) with Monte Carlo
(MC) approaches, applied to the simulation of catalytic and electrodeposition processes.11–13 Indeed, the MC approach consists in an
elegant simulation technique which allows closing the gap between
the atomistic and the continuum viewpoints of the mechanisms under
investigation.
In contrast to Einstein and his statement that “God does not play
dice”, Hawking believes that the Universe is being governed by nondeterminism (chance) as illustrated by his answer to Einstein “Not
only does God definitely play dice, but He sometimes confuses us
by throwing them where they cannot be seen”.14 The MC approach
∗
z
Electrochemical Society Active Member.
E-mail: [email protected]
exploits this fact by extensively repeating random executions to obtain
results that can mimic those physical systems where the random character is inherent. Several MC techniques have been developed so far,
the most popular one being the Metropolis MC algorithm consisting
on performing random swaps from a given arrangement (e.g. spatial
distribution of particles constituting a system) in order to search for
the global minimal energy configuration.15–17
The so-called Kinetic Monte Carlo (KMC) approach is a MC
method that helps to study the evolution of the systems where the
associated transition rates are known. The KMC approach treats the
events as a simple Markov process, in other words, a series of “memoryless” random swaps. The probability density function Pi , that is
the probability to find the system at state i, is described by the master
equation of Markovian form18
d Pi (t)
=−
[1]
ki j Pi (t) +
k ji P j (t)
dt
j=i
j=i
where kij is the average escape rate from state i to state j. The main
KMC postulate is that during thermal vibrational motion the system
loses memory and has the same probability of finding the escape path
from state i to state j for any short time increment t, leading to an
exponential decay statistics.12 This is mathematically described by the
Poisson distribution
[2]
Pi (t) = ki j exp −ki j t
In particular, the KMC approach is relevant for investigating electrochemical reactions in devices for energy conversion and storage
which show sensitivity to the microstructure of the active materials.19
In the frame of Lithium Ion Batteries (LIBs), Methekar et al.20 applied
a KMC simulation approach to investigate the formation of the Solid
Electrolyte Interphase (SEI) layer on an intercalation graphite anode.
Van der Ven and Ceder21 used KMC simulations to study lithium
diffusion through a divacancy mechanism in layered Lix CoO2 . Yu
et al.22 implemented KMC simulations parameterized by moleculardynamics-based activation energy barriers for Li+ and e− diffusion in
TiO2 , revealing the central role of the electrostatic coupling between
Li+ and e− on their collective drift diffusion.
In the frame of the Polymer Electrolyte Membrane Fuel Cells
(PEMFCs),23 the description of the relationship between the catalyst/electrolyte interfacial charge distribution and its evolution represents a significant issue for the understanding of the redox kinetic
Downloaded on 2015-04-01 to IP 31.33.112.81 address. Redistribution subject to ECS terms of use (see ecsdl.org/site/terms_use) unless CC License in place (see abstract).
E74
Journal of The Electrochemical Society, 162 (7) E73-E83 (2015)
Figure 1. Schematics of the classical VSSM
principles.
processes and the prediction of the associated effective catalyst activity. The so-called first-principles methods, like the DFT approach,
constitute interesting tools to study isolated redox events and can offer
powerful capabilities for predicting materials properties when coupled
with thermodynamics.24 However, the amount of competing events
and involved species in complex reactions such as the Oxygen Reduction Reaction (ORR) makes these methods unable to provide a wide
description of the redox events in realistic electrochemical operation
conditions where both the electric field and the temperature are not
zero. In this sense, in order to describe the electrode surface coverage
dynamics in relation to processes like adsorption, desorption and reactions, the implementation of a KMC method arises as a highly relevant
technique to get a proper description of surface electrocatalytic events.
Within this sense, Zhdanov23 implemented the KMC approach
to study the O2 reduction on Pt(100) finding that the impact in the
overall kinetics of the lateral interaction between adsorbates is not
so significant as the O2 adsorption rate. By using KMC simulations,
Abramova et al.25 investigated the O2 adsorption on different metallic
surfaces with fcc(100) termination as function of the metal surface
charge density and the molecular oxygen sticking probability.
The KMC approach was also used to study the electric field fluctuations in simple electrochemical reactions kinetic rates due to electrochemical double layer (EDL) effects.26 Furthermore, Harrington27
implemented the KMC approach to study the Pt oxide growth.
Among the KMC algorithms, one of the most popular is the
so-called Variable Step Size Method (VSSM),12,28 which has been
used to describe surface reactions in catalysis (Figure 1). The VSSM
algorithm
(a)
(b)
starts with a given species surface spatial distribution (called
“configuration”);
then, it creates a list of all possible M processes (M = 6 in the
example of Figure 1) with the corresponding transition rates, and
M
calculates ktot = ki ;
(c)
then it generates two random numbers ρ1 and ρ2 and selects a
single process N that fulfills
i
N
i=1
(d)
(e)
(f)
ki ≥ ρ1 ktot ≥
N −1
ki
[3]
i=1
then, it executes the selected process;
consequently, it updates the clock with a time step given by
t = ln(ρ2 )/ktot ;
if the time t < the desired total simulation time, then it backs to
(b), otherwise, it stops.
One common expression for the kinetic rate ki in equation 3 has
the Arrhenius form
Ei
[4]
ki = κi exp −
kB T
where κi is a pre-exponential factor depending on the vibrational
frequency of the transition state, Ei is the activation barrier between
the two states i and j while kB and T are the Boltzmann constant and
the absolute temperature, respectively.
Casalongue et al.29 described the dynamics of the ORR on Pt(111)
surfaces with the VSSM. Their simulations included an electrode
Butler-Volmer potential as an external parameter being controlled
(input parameter). A similar approach was implemented by other
authors.19,30,31
These previous works lack of predicting electrochemical conditions and full polarization curves (I-V) as they do not consider the
EDL effects, thus the surface charge density impact onto the electrochemical elementary kinetic steps, neither the reactants transport at
the vicinity of the catalyst.
In this paper we present a new multiparadigm modeling approach
to simulate electrochemical reactions in materials for energy conversion or storage. The approach combines a KMC algorithm, called
here “Electrochemical-VSSM” (E-VSSM), with continuum transport models. For instance, in this way, it allows simulating electrochemical conditions by accounting for reactants/products adsorption/
desorption dynamics, reaction intermediates surface diffusion,
chemical/electrochemical reactions, EDL effects and mesoscopic
reactants/products transport phenomena. Thus, this approach allows
capturing the dynamics of the reactants, reaction intermediates and
products at the atomistic/molecular level for an operating electrochemical cell. These processes cannot be addressed with state of the
art continuum kinetic models supported only on the MF approximation which neglects, by construction, adspecies surface diffusion
phenomena for example.32,33
This paper is organized as follows. First our model is presented
within the framework of PEMFCs. An illustrative application example
is provided for the description of the ORR kinetics in a Pt(111)-based
PEMFC cathode. Then results under different simulated scenarios are
discussed. Finally, we conclude and indicate further directions for our
model development and application.
Our Model
For the case of PEMFCs, our model is articulated into interconnected modules representative of three relevant scales in a PEMFC
cathode electrode (Figure 2a):
Downloaded on 2015-04-01 to IP 31.33.112.81 address. Redistribution subject to ECS terms of use (see ecsdl.org/site/terms_use) unless CC License in place (see abstract).
Journal of The Electrochemical Society, 162 (7) E73-E83 (2015)
E75
Figure 2. (a) Schematics of our general multiscale modeling framework; (b) schematics of the O2 mesoscopic transport model coupled with the EDL model;
(c) schematics of the EDL model with some of the adspecies represented.
r mesoscopic scale describing reactants transport in gas phase in
a parallel direction to the catalyst surface (Figure 2b);
r nanoscopic ionic transport scale in the electrolyte (constituted
of ionomer, water and protons) through a continuum model proposed
by us in Ref. 34 (Figure 2c) and describing the external layer (EL)
within the EDL region;
r the adlayer scale or inner layer (IL), where the redox reactions
take place, and described through a KMC modeling framework aiming
to rationalize the effects of the catalyst surface heterogeneity on the
net electrochemical reaction rate (Figure 3).
Mesoscopic scale.— For the simulation results reported in this
paper, as an illustrative example, O2 transport in gas phase at the
mesoscopic scale is described as a Darcy-like process in the parallel direction to the catalyst surface. The corresponding conservation
equation in the discretized form, for a system of M bins along the Y
Downloaded on 2015-04-01 to IP 31.33.112.81 address. Redistribution subject to ECS terms of use (see ecsdl.org/site/terms_use) unless CC License in place (see abstract).
E76
Journal of The Electrochemical Society, 162 (7) E73-E83 (2015)
Figure 3. Schematics of our E-VSSM principles.
direction (Figure 2b), is given by
pn,O2
=k
t
in
Jn,O
2
−
out
Jn,O
2
y
− ωkmc
n,O2
[5]
where the time evolution of the O2 pressure in each bin (n) is solved
by numerical integration.
In Eq. 5 pn,O2 is the oxygen pressure in bin (n), ωkmc
O2 is the oxygen
consumption rate calculated with KMC in the same bin and expressed
in atm/second, k is the O2 transport coefficient at the mesoscopic scale,
out
in
is the oxygen flow rate going out from the bin (n) and Jn,O
is
Jn,O
2
2
the oxygen flow rate going in the bin (n) defined respectively by
in
= ( pn−1,O2 − pn,O2 )
Jn,O
2
[6]
out
= −( pn+1,O2 − pn,O2 )
Jn,O
2
[7]
where the O2 partial pressures p0 and p N are the boundaries conditions
given in Table II. Once the gas O2 pressure is calculated for each mesh,
the corresponding absorbed O2 concentration in the electrolyte (EDL
region) is calculated from the Henry’s law c O2 = p O2 /α H , where α H
is the Henry’s constant.7
As a first approximation, we neglect H2 O transport in all the scales
as we suppose the EL (see below) to be fully hydrated. Furthermore,
H2 O2 transport is also neglected as its calculated production is very
small for the electrocatalytic system under investigation in this paper.
We underline however that the implementation of these transport models in the overall framework is quite straightforward and is reserved
for a future investigation.
External layer.— The EL model is a continuum model which
calculates the proton molar concentration at the reaction plane, i.e.
x = L, for each KMC time step. This EL model is based on the transport model in EDLs we have published in Ref. 34. According to this
model, the molar concentration ci in the electrolyte of a specie i is
given by the solution of
∂ci
2
ai j c j
= Bi RT ∇ ci + Bi ∇ · ci ∇
∂t
i
· ci ∇ϕ
− N A Bi ∇
· ci ∇
pi · ∇ϕ
+ Bi z i F ∇
+ f 1 (βi , ci ) − f 2 (βi , ci )
[8]
In the above expression f1 and f2 are terms capturing the finite size
and relative finite size entropic diffusion and their mathematical formulations are given in Ref. 34. In Eq. 8 Bi is the i specie diffusion
coefficient over RT where R is the ideal gas constant, NA is the Avogadro number, ϕ is the electrostatic potential, zi is the i specie charge,
pi is the i specie dipolar moment, aij is the quantum interaction energy
between the i and the j species, βi is the i specie relative size related
to the smallest specie.
For illustrative purposes, we consider here the steady-state solution
of the Eq. 8 ( ∂c∂ti = 0) with a single species size (f1 and f2 = 0) for H+ ,
for the water molecule and for counter ions concentration (SO− 3 ) of
the ionomer (e.g. Nafion), then Eq. 8 becomes for protons,
· c H + ∇ϕ
0 = B H + RT ∇ 2 c H + + B H + z H + F ∇
[9]
for water
and for
· c H2 O ∇
p H2 O · ∇ϕ
0 = −N A B H2 O ∇
[10]
· c S O − ∇ϕ
0 = BS O3− RT ∇ 2 c S O3− + BS O3− z S O3− F ∇
3
[11]
S O−
3
Because the SO− 3 and water concentrations are assumed to be uniform
along the EL thickness (L), expressions [10] and [11] arise to the
relation for the electrostatic potential
ϕ = a1 x + a2
[12]
S O−
3
and water concenwhere a1 and a2 are constants related to the
trations and calculated from the boundary conditions.
In one dimension Eq. 9 becomes
∂ϕ
∂c H +
− B H + z H + Fc H +
[13]
∂x
∂x
where JH + is an integration constant related to the proton flux at the
catalyst surface.
We can then rewrite Eq. 13 as
Fϕ ∂
Fϕ
c H + exp
[14]
JH + = B H + RT exp −
RT ∂ x
RT
JH + = −B H + RT
so
J
H+
exp
Fϕ
RT
=B
H+
∂
Fϕ
+
c H exp
RT
∂x
RT
[15]
Downloaded on 2015-04-01 to IP 31.33.112.81 address. Redistribution subject to ECS terms of use (see ecsdl.org/site/terms_use) unless CC License in place (see abstract).
Journal of The Electrochemical Society, 162 (7) E73-E83 (2015)
As we assume steady-state conditions (JH + uniform along the coordinate X), by integrating in the EL thickness L we obtain
c H + (L) exp Fϕ(L)
− c H + (0) exp Fϕ(0)
RT
RT
JH + = −B H + RT
[16]
L
d
x
exp Fϕ(x)
RT
Table I. List of activation energies used in the KMC module.
Activation
energy [kJ.mol−1 ]
Reaction
Forward
Backward
O2ads + s ←→ 2 Oads
30
149
0
In Eq. 12 we have an expression for the electrostatic potential profile
(L) = E
and by fixing the boundary conditions as ϕ(0) = 0 and dϕ
dx
(where E is the electric field that is obtained at the IL boundary from
the surface charge density σ in the same way as it has been done
in Ref. 34) we can explicitly solve Eq. 16 and calculate the proton
concentration at x = L as,
⎛
⎞
L
+
exp[Fϕ(x)]d
x
J
H
⎜
⎟
⎜
⎟
0
c H + (L) =⎜
+ c H + (0) exp[Fϕ(0)]⎟ exp[−Fϕ(L)]
⎝
⎠
B H + RT
H+
Internal layer.— Once again we follow the strategy reported in
our Ref. 34 where a correction in the activation energy is introduced
in Eq. 4 for all those kinetic parameters k where a charge transfer is
occurring, thus
kB T
−E act,l ∓ f (σ)
[18]
exp
kl = c H + (L) κ
h
RT
where c H + (L) is the proton concentration at the boundary between IL
and EL (x = L); Eact,l is the activation energy of an elementary reaction
l without influence of the interfacial electric field (extracted from DFT
calculations), h is the Planck constant. f(σ) is a function of the charge
density σ, the so-called surface potential, a quantity different from zero
only for electrochemical steps. f(σ) is also equal to the electrostatic
potential drop through the IL between the catalyst surface and the
electrolyte.34 We note that in the exponential argument Eact,l ± f(σ) is
an effective activation energy (the sign depends on whether reduction
or oxidation is occurring).
The charge density σ is calculated in terms of the imposed electronic current density J (input for our model) through
∂σ
[19]
∂t
where JFar is the faradaic current density which is equal to FJH+ and
calculated here as the total charge δQ transferred by the electrochemical reactions per unit of time δt and per unit of area A in the catalyst
surface:
1 ∂Q
JF A R =
[20]
A ∂t
Note that JFAR impacts the time evolution of σ which affects in turn
the kinetic rates in Eq. 18. The expression for f(σ) is given by34
J − JF A R = −
F
(σ − ) H
[21]
N A ε0
where F is the Faraday constant, ε0 the electric permittivity of vacuo
and is the adlayer dipolar charge density calculated by adding the
dipolar moment of all the adspecies divided by the catalyst surface
area (as a first approximation, only the dipolar moments of H2 O and
OH have been considered in this paper). The adlayer dipolar density
screens the effective charge density, and H is the IL thickness (the
distance between the catalyst surface and the assumed location of the
reaction plane).
f (σ) =
+ e−
k2 /k−2
+ O2ads ←→ HO2ads
38
43
2OHads ←→ H2 Oads + Oads
1
158
H+ + e− + Oads ←→ OHads
88
94
19
80
24
47
k3 /k−3
k4 /k−4
H+
+ e−
k5 /k−5
+ OHads ←→ H2 Oads
k6 /k−6
H+ + e− + HO2 ads ←→ H2 O2 ads
k7 /k−7
O2ads diffusion
Oads diffusion
OHads diffusion
H2 Oads diffusion
H2 O2ads diffusion
O2 Hads diffusion
[17]
Because of the small value of L (typically few nanometers), O2 concentration in the EL is assumed to be uniform.
Since H2 O transport is neglected an uniform and constant electric
permittivity is assumed for the EL. A description of the dependence
of this electric permittivity on the space can be included by following
our model in Ref. 34, something that we plan to do in close future.
E77
30
40
45
20
20
20
The algorithm behind the E-VSSM proposed in this paper is structurally similar to the VSSM but with some major differences:
(a
)
(b
)
(c
)
(d
)
(e
)
(f
)
first, the system starts with a given surface charge density (zero
in the present work) and initial configuration (free metal surface
in the present work);
then the total kinetic rate is calculated according to Eq. 18
and 21;
this step remains the same as in the classical VSSM (step (c));
when the algorithm executes the process it also calculates the
total charge transferred;
once the time step is obtained in the same way as in the classical
VSSM (step (e)), JFAR is calculated from Eq. 20, and for the given
imposed current density J, the new charge density is calculated
by numerical integration of Eq. 19;
this step remains the same as in the classical VSSM (step (f)).
The O2 and the H2 O adsorptions are modeled with the same
methodology that has been used in Ref. 32 and 35. The activation
barriers used in the KMC module for the considered ORR elementary
steps, are presented in Table I and the parameters values used in the
continuum modules are listed in Table II. All the reaction activation
barriers were taken from Refs. 35 and 36 and the diffusion activation
barriers were assumed. The initial configuration was set as a perfectly
clean Pt(111) surface but non-zero initial coverage can also be assumed for O2 and H2 O. Because the simulation starts with a clean
surface, the only possible steps at the beginning are the competitive
Table II. Parameter values implemented in the continuum
modules.
Parameter
Value
κ
O2 sticking coefficient
Proton diffusion coefficient
H2 O dipolar moment
OH dipolar moment
EL thickness
IL thickness
ε0
Temperature
Inlet pressure p0
Outlet pressure p N
7.29 × 10−6 s−1
0.045a,b
1.0 × 10−9 m2 s−1
5.6 × 10−10 C.m
6.7 × 10−10 C.m
5 × 10−9 m
2 × 10−10 m
8.854 × 10−12 F.m−1
350 K
1.2 atm
1.0 atm
a Reference 37.
b Reference 38.
Downloaded on 2015-04-01 to IP 31.33.112.81 address. Redistribution subject to ECS terms of use (see ecsdl.org/site/terms_use) unless CC License in place (see abstract).
E78
Journal of The Electrochemical Society, 162 (7) E73-E83 (2015)
adsorption of H2 O and O2 . The H2 O concentration was fixed constant
(55000 mol.m−3 ).
a
Computational Details
The KMC part of the model has been implemented within the
in-house LRCS software MESSI (Monte Carlo Electrochemistry
Simulation Software for Innovation), fully coded in Python programming language. The continuum parts of the model (EL and mesoscale
transport modules) are implemented within the MS LIBER-T framework which is also coded in Python (www.modeling-electrochemistry.
com). All the simulations were carried out in our University calculation platform Plateforme Mod´elisation et Calcul Scientifique.39 Typical calculation time for one single simulation on one single processor
is comprised between ten hours for a single constant current simulation and three days for the calculation of a full polarization curve. The
electrode potential of the cathode is obtained by the addition of the
electrostatic potential jump across the IL, equal to f (σ), and the electrostatic potential at the reaction plane calculated from the equation
12 in the EL.7,34
The Pt(111) surface was modelled with a supercell with 30 × 30
atoms unit cells (with a cell parameter a = 2.81 Å) with periodic
boundary conditions in order to guarantee simulation reproducibility
in a reasonable computational time. For all the results presented in
the following, each simulation was run at least five times, in order to
check that no significant deviation is occurring due to the probabilistic
distribution.
The EL thickness was fixed at 5 nm. The mesoscopic O2 transport
module was numerically solved with a single bin in most of the cases,
except for the polarization curve simulation where it was discretized
in three bins. The solved mesoscale transport length is 0.09 m. The
temperature for the simulations was fixed at 350 K unless otherwise
stated.
Open circuit condition.— In this section we report results with
zero imposed current density (J = 0 mA.cm−2 in Eq. 19) in order to
simulate the Open Circuit Voltage (OCV) condition.
The initial condition of the surface charge density, for all the
simulated potentials, is set to be zero. It has been demonstrated
that for uncharged surfaces adsorbed water molecules adopt a “flat”
configuration40 (their dipolar moments prefer being parallel to the
surface). Once the surface becomes charged through the protons
reduction the water molecules reorient and their dipolar moment
become perpendicular to the surface. As our model does not describe
the spatial configuration changes of each adspecie, it implements
an algorithm in order to hold at zero the polarization field of the
adsorbed water molecules at zero charge density until the first proton
reduction event occurs.
From these simulations we observe that, when some of the electrochemical steps start to occur, the surface becomes positively charged
and the argument of the exponential in Eq. 21 becomes negative;
in other words, electrons involved in the reduction reactions are attracted by the catalyst surface which increases the effective activation
barrier. Once the exponential argument becomes sufficiently negative
the electrochemical steps become inhibited and no more charge transfer occurs, leading to an equilibrium regime. The calculated species
coverage evolutions and the corresponding final configuration are presented in Figure 4.
As illustrated in Figure 5 for a simulation representative of OCV
conditions, our model is able to track the IL electric permittivity
evolution. Notice that for previous IL models (e.g. Ref. 32) treat
the electric permittivity as a constant value fitted from macroscopic
experimental data.
As one can see in Figure 5, the IL electric permittivity is not
constant throughout the time. Indeed, the electric permittivity depends
on the dipolar moment of the adsorbed species, their coverage and the
surface charge density, thus on the applied external current density.
The calculated OCV is reported in Figure 6.
b
Figure 4. a) Calculated coverage evolution of the ORR species under the
OCV condition; b) final obtained configuration.
In Figure 7, we investigate the impact of temperature onto the
OCV and the ORR intermediates surface coverage evolution (zero applied current density). One can see that the calculated OCV increases
with the temperature (Figure 7a) due to the increase of the effective
kinetic rates (cf. Eq. 18). Additionally we observe that the calculated
water coverage increases with temperature which leads to a slight O2
coverage decrease.
Because of its atomistic nature, the KMC approach allows investigating the impact of surface defects on the effective electrochemical
activity. This is due to the fact that it explicitly accounts for the
structural correlations between the adspecies locations and the spatial
configuration of active sites. This is in contrast to the MF approach
where it is assumed that all the adspecies are perfectly mixed (no intermediates surface diffusion limitations). We report here simulation
results at OCV conditions for a Pt(111) catalyst with surface inactive
sites randomly distributed to mimic the properties of a polycrystalline
or degraded surface (e.g. due to electrochemical dissolution), noted in
the following Pt∗ .
In Figure 8 we can see that the presence of inactive sites does not
only impacts the coverage evolution but also the final coverage of the
adspecies. The calculated final water coverage is 0.13 ML for Pt(111)
surface and 0.18 ML for Pt∗ . This result is due to the presence of the
inactive sites making less favorable the adspecies surface diffusion
and consequent reactions.
Calculated response to an applied current density step.— In this
section we study the influence of applied current density steps on the
dynamical response of the potential and the associated intermediate
coverage evolution. These current density steps are applied once the
OCV state regime is reached.
The calculated apparent reactivity results from the interplay between kinetic processes and surface transport of adspecies. Calculated ORR adspecies coverage evolutions for a 0 to 1 mA/cm2 current
Downloaded on 2015-04-01 to IP 31.33.112.81 address. Redistribution subject to ECS terms of use (see ecsdl.org/site/terms_use) unless CC License in place (see abstract).
Journal of The Electrochemical Society, 162 (7) E73-E83 (2015)
Figure 5. a) Calculated relative electric permittivity evolution at zero applied
current density (OCV condition); b) calculated surface charge density and
dipolar screening charge density evolutions.
E79
Figure 7. a) Calculated potential evolution and b) associated O, O2 and H2 O
coverage evolution, at 290 and 350 K at zero applied current density (OCV
condition).
density step applied at 10−2 seconds (once the stationary OCV is
reached) are presented in Figure 9. After the application of the current
density step, most of the ORR adspecies become significant in the
surface, in contrast to the OCV case (Figure 4), where only O, O2 and
H2 O are significantly present in the surface.
Figure 6. Calculated potential evolution at zero applied current density (OCV
condition).
Figure 8. (a) coverage evolution along the time in OCV condition; (b) calculated adlayer final configurations for Pt(111) and for Pt∗ .
Downloaded on 2015-04-01 to IP 31.33.112.81 address. Redistribution subject to ECS terms of use (see ecsdl.org/site/terms_use) unless CC License in place (see abstract).
E80
Journal of The Electrochemical Society, 162 (7) E73-E83 (2015)
Figure 10. Calculated potential evolutions following the current steps from
0 to 1 and 5 mA/cm2 applied at 10−2 seconds.
splitting into two O atoms as less free sites to receive the O atoms
are available. Other adspecies coverage show also lower coverage for
Pt∗ , H2 O (0.76 vs. 0.64), OH (0.17 vs. 0.13) and O2 H (0.11 vs. 0.08).
Indeed, the inactive sites make more difficult the adspecies surface
diffusion reducing the chances to have lateral reactions.
Figure 9. a) Calculated coverage evolution of the ORR adspecies with a 0 to
1 mA/cm2 current density step applied at 10−2 seconds; b) calculated final
configuration.
The application of the current density step leads to the increase
of the H2 O, O2 H and OH production and the O2 and O coverage
decrease vs. the OCV case. This is expected, since, H2 O, O2 H and
OH productions need a non-zero current density to form through H+
reduction. Another significant result is that the system shows a delay
in its response following the applied current density step, as it can be
seen for the coverage evolution. Indeed, when the current density step
is applied, the potential does not change instantaneously. Instead it
changes once the surface charge density is sufficiently high to induce
a sufficient number of electrochemical steps.
In Figure 10, we report the calculated potential evolution for two
current steps: from 0 to 1 and from 0 to 5 mA/cm2 , both applied at
10−2 seconds. The noisy aspect observed for the calculated potential
is due to the stochastic character of the faradaic current and the charge
density dynamics.
In Figure 11, we show the calculated ORR species coverage evolution for the case of a current density step from 0 to 5 mA/cm2 applied
at 10−2 seconds. Once again, other species become significant as O2 H
and OH and the H2 O production increase depleting the O2 and O
coverage.
In the following we investigate the system response to a current
step from 0 to 5 mA/cm2 applied at 10−2 seconds for a surface with inactive sites randomly distributed (called Pt∗ ) and mimicking partially
degraded Pt, and we compare the results with the ones obtained with
Pt(111) at the same conditions.
The calculated ORR species coverage evolution is illustrated in
Figure 12 in logarithmic scale. One can see that the inactive sites
impact significantly into the arising coverage especially for O, where
we find 0.17 ML coverage for Pt(111) vs. 0.06 ML coverage for Pt∗ .
Indeed the presence of the inactive sites makes more difficult the O2
b
Figure 11. Calculated coverage evolution of the ORR species following an
applied current step from 0 to 5 mA/cm2 at 10−2 seconds; b) calculated adlayer
final configuration.
Downloaded on 2015-04-01 to IP 31.33.112.81 address. Redistribution subject to ECS terms of use (see ecsdl.org/site/terms_use) unless CC License in place (see abstract).
Journal of The Electrochemical Society, 162 (7) E73-E83 (2015)
E81
Figure 13. calculated (in black) vs. experimental (in blue) I-V curve (experimental curve extracted from Ref. 41).
can see that the resulting O2 concentration distribution only impacts
significantly in the OH production, but not in the water production. We
also notice that the variation in the O2 pressure through the mesoscale
does not impact significantly into the ORR intermediates coverage
dynamics. The impact of larger pressure gradients and water transport
onto this dynamics merit further investigation.
Conclusions
Figure 12. Calculated coverage evolution in logarithmic scale of the ORR
adspecies following a current step from 0 to 5 mA/cm2 applied at 10−2 seconds
for the Pt∗ surface (b) and comparison with the Pt(111) results (a).
Polarization curve.— In the present section the mesoscopic transport model has been discretized in three identical bins. The external O2
input pressure was set in 1.2 atm and the external O2 output pressure
was set in 1 atm.
In order to simulate the polarization curve I-V, we scan the input
current density between 0 and 5 mA/cm2 with incremental current
density steps of 0.5 mA/cm2 . First we allow the potential to stabilize at
0 mA/cm2 for 0.1 seconds, and then the first current density increment
is applied. Then the system is allowed to relax again for another
0.1 seconds. Once the steady state is reached, a new current density
increment is applied. The procedure is repeated until the whole range
of current densities is explored. As shown in Figure 13, our calculated
polarization curve shows good agreement with the experimental data
obtained by Markovic et al.41 on Pt(111).
In Figure 13, the potential drop close to 4 mA/cm2 is due to water
saturation at the surface. In other words, since the reduction of H+ is
prevented because of the lack of O2 and O at the surface, the Faradaic
current no longer compensates the imposed current density J (Eq. 19),
and the surface charge density changes significantly.
The corresponding coverage of ORR species evolution in time all
along the applied current density steps for each mesh is presented in
Figure 14 in logarithmic scale.
The calculated steady state O2 pressure for each bin is 1.16 atm
(mesh #1), 1.05 atm (mesh #2) and 1.02 atm (mesh #3). In Figure 14 we
In this paper we presented a new multi-paradigm simulation framework of electrochemical systems coupling an atomistically resolved
model of electrochemical reactions with continuum models describing transport of charges and reactants in the electrolyte. To the best
of our knowledge, a model with these characteristics has not been
previously reported.
The atomistically resolved model is based on a new KMC algorithm developed by us (called MESSI), extending the VSSM method
by accounting explicitly of the electrochemical conditions. This arises
in particular through the coupling of the KMC algorithm describing the
reactions with a non-equilibrium electrolyte/active material interface
electrochemical double layer model based on statistical mechanics
and reported by us in Ref. 34. The continuum models are supported
on a set of coupled differential equations solved in the pre-existing
MS LIBER-T simulation package available at LRCS.
This simulation framework allows us to investigate how the macroscopic electrode performance responds to nanoscopic and atomistic
scale mechanisms by capturing details such as the surface morphology, inactive sites distribution, nanoparticle facets and sizes. In other
words, this approach allows capturing the dynamics of the reactants,
reaction intermediates and products at the atomistic/molecular level
for an operating electrochemical cell. This phenomena cannot be addressed with state of the art continuum kinetic models supported only
on the MF approximation which neglects, by construction, adspecies
surface diffusion phenomena for example. As an application example, a Pt(111)-based PEMFC cathode was studied. The model allows
predicting simultaneously electrochemical observables (potential vs.
time, polarization curves. . . ) and the associated surface intermediates reaction coverage. We have found that O, H2 O and OH are the
dominant adsorbed species in good agreement with the literature.29,42
Note that in comparison to previous works, in the present work more
elementary steps are considered (Table I) which allow the production
of other species like O2 H and H2 O2 .
This approach provides new capabilities for investigating the
interplaying between the EDL structure and electrochemical reactions since our simulations are carried out at galvanostatic
Downloaded on 2015-04-01 to IP 31.33.112.81 address. Redistribution subject to ECS terms of use (see ecsdl.org/site/terms_use) unless CC License in place (see abstract).
E82
Journal of The Electrochemical Society, 162 (7) E73-E83 (2015)
b
Figure 14. Calculated coverage evolution in logarithmic scale for the three bins (a) located at 2, 5 and 8 cm; (b) three meshes model scheme.
or galvanodynamic conditions, i.e. we do not impose cell potential
and constant capacities. Furthermore, the electric permittivity of the
IL is calculated on the fly from the coverage of the adsorbing and
forming species on the surface.
Because of the generality and the flexibility of the computational
modeling algorithm proposed in this paper, it is readily applicable
for other reactions where material heterogeneity plays an important
role. Indeed, we believe that the approach can offer interesting capabilities for the investigation of the ORR mechanism in lithium air
batteries or the electrochemical processes in lithium sulfur battery
cathodes.2,43,44
Acknowledgments
The authors deeply acknowledge the funding by the European
Project PUMA MIND, under the contract 303419. Prof. Dominique
Larcher (LRCS, Amiens, France), Dr. David Loffreda, Dr. Philippe
Sautet and Dr. Federico Calle Vallejo (ENS, Lyon, France) are also
acknowledged for fruitful discussions.
References
1. A. A. Franco, PEMFC degradation modeling and analysis, in Polymer electrolyte
membrane and direct methanol fuel cell technology (PEMFCs and DMFCs) - Volume
1: Fundamentals and performance, edited by C. Hartnig and C. Roth (Woodhead,
UK) (2012).
2. A. A. Franco and K.-H. Xue, ECS J. Solid State Sci. Technol., 2, M3084
(2013).
3. S. Strahl, A. Husar, and A. A. Franco, Int. J. Hydrog. Energy, 39, 9752 (2014).
4. A. A. Franco, Multiscale modeling, in Encyclopedia of Applied Electrochemistry,
edited by R. Savinell, K. I. Ota, G. Kreysa (Springer, UK) (2013).
5. A. A. Franco, Polymer Electrolyte Fuel Cells: Science, Applications, and Challenges,
CRC Press (2013).
6. A. A. Franco, RSC Adv., 3, 13027 (2013).
7. A. A. Franco, P. Schott, C. Jallut, and B. Maschke, Fuel Cells, 7, 99 (2007).
8. D. Sheppard, R. Terrell, and G. Henkelman, J. Chem. Phys., 128, 134106 (2008).
9. K. Malek and A. A. Franco, J. Phys. Chem. B, 115, 8088 (2011).
10. L. F. L. Oliveira, C. Jallut, and A. A. Franco, Electrochimica Acta, 110, 363
(2013).
11. L. Madec, L. Falk, and E. Plasari, Chem. Eng. Sci., 56, 1731 (2001).
12. K. Reuter, in Modeling and Simulation of Heterogeneous Catalytic Reactions: From
the Molecular Process to the Technical System, edited by O. Deutschmann, Weinberg:
Wiley-VCH (2011).
13. X. Li et al., J. Electrochem. Soc., 154, D230 (2007).
14. S. W. Hawking, A Brief Story of Time from the Big Bang to Black Holes. Bantam
Book, New York (1988).
15. J. P. Valleau and D. N. Card, J. Chem. Phys., 57, 5457 (1972).
16. W. Foulkes, L. Mitas, R. Needs, and G. Rajagopal, Rev. Mod. Phys., 73, 33 (2001).
17. B. Yang, M. Asta, O. N. Mryasov, T. J. Klemmer, and R. W. Chantrell, Scr. Mater.,
53, 417 (2005).
18. N. G. Van Kampen, Stochastic processes in physics and chemistry, Access Online
via Elsevier, (1992).
19. B. Andreaus and M. Eikerling, J. Electroanal. Chem., 607, 121 (2007).
20. R. N. Methekar, P. W. C. Northrop, K. Chen, R. D. Braatz, and V. R. Subramanian,
J. Electrochem. Soc., 158, A363 (2011).
21. A. V. der Ven and G. Ceder, Electrochem. Solid-State Lett., 3, 301 (2000).
22. J. Yu, M. L. Sushko, S. Kerisit, K. M. Rosso, and J. Liu, J. Phys. Chem. Lett., 3, 2076
(2012).
Downloaded on 2015-04-01 to IP 31.33.112.81 address. Redistribution subject to ECS terms of use (see ecsdl.org/site/terms_use) unless CC License in place (see abstract).
Journal of The Electrochemical Society, 162 (7) E73-E83 (2015)
23. V. P. Zhdanov, J. Electroanal. Chem., 607, 17 (2007).
24. J. K. Nørskov et al., J. Phys. Chem. B, 108, 17886 (2004).
25. L. A. Abramova, A. V Zeigarnik, S. P. Baranov, and E. Shustorovich, Surf. Sci., 565,
45 (2004).
26. V. P. Zhdanov and B. Kasemo, Appl. Surf. Sci., 219, 256 (2003).
27. D. A. Harrington, J. Electroanal. Chem., 420, 101 (1997).
28. A. P. J. Jansen, Comput. Phys. Commun., 86, 1 (1995).
29. H. S. Casalongue et al., Nat. Commun., 4 (2013).
30. V. Rai, M. Aryanpour, and H. Pitsch, J. Phys. Chem. C, 112, 9760
(2008).
31. J. J. Lukkien, J. P. L. Segers, P. A. J. Hilbers, R. J. Gelten, and A. P. J. Jansen, Phys.
Rev. E, 58, 2598 (1998).
32. R. F. de Morais, P. Sautet, D. Loffreda, and A. A. Franco, Electrochimica Acta., 56,
10842 (2011).
33. D. Eberle and B. Horstmann, Electrochimica Acta, 137, 714 (2014).
34. M. Quiroga, K. H. Xue, T. K. Nguyen, H. Huang, M. Tulodziecki, and A. A. Franco,
J. Electrochem. Soc., 161, E3302 (2014).
E83
35. R. Ferreira de Morais, Study of the stability and the reactivity of Pt and Pt3 Ni
model catalyst for PEM fuel cells: an ab-initio based multiscale modeling ap´
proach, PhD thesis, Ecole
Normale Sup´erieure, Lyon, (2011) http://www.theses.
fr/2011ENSL0694.
36. R. Ferreira de Morais, A. A. Franco, P. Sautet, and D. Loffreda, ACS Catal., 5, 1068
(2015).
37. K. Griffiths and D. Bonnett, Surf. Sci., 177, 169 (1986).
38. D. R. Monroe and R. P. Merrill, J. Catal., 65, 461 (1980).
39. https://www.mecs.u-picardie.fr/doku.php.
40. S. Meng, E. G. Wang, and S. Gao, Phys. Rev. B, 69, 195404 (2004).
41. N. M. Markovic, T. J. Schmidt, V. Stamenkovic, and P. N. Ross, Fuel Cells, 1, 105
(2001).
42. V. Viswanathan et al., J. Phys. Chem. C, 116, 4698 (2012).
43. K. H. Xue, E. McTurk, L. Johnson, P. G. Bruce, and A. A. Franco, J. Electrochem.
Soc., 162(4), A614 (2015).
44. K. H. Xue, T. K. Nguyen, and A. A. Franco, Journal of the Electrochemical Society,
161(8), E3028 (2014).
Downloaded on 2015-04-01 to IP 31.33.112.81 address. Redistribution subject to ECS terms of use (see ecsdl.org/site/terms_use) unless CC License in place (see abstract).