Generalization of Porous Electrode Theory for Noninteger

29 Dec 2017 - Conventional porous electrode theories are formulated in integer dimensional space; they employ the volume averaging method to simplify ...
3 downloads 0 Views 2MB Size
Article pubs.acs.org/JPCC

Cite This: J. Phys. Chem. C XXXX, XXX, XXX−XXX

Generalization of Porous Electrode Theory for Noninteger Dimensional Space Jun Huang* College of Chemistry and Chemical Engineering, Central South University, Changsha 410083, PR China ABSTRACT: Conventional porous electrode theories are formulated in integer dimensional space; they employ the volume averaging method to simplify consideration of structural complexities and introduce volume fraction and tortuosity as two effective structural descriptors. In this work, I generalize the porous electrode theory for noninteger dimensional space using fractional calculus, in which fractional dimensionalities are invoked as quintessential structural descriptors. The present theory furnishes a theoretical framework for describing mass transport and reactions in fractal electrodes, which is of practical implications for modeling of porous electrodes in electrochemical energy devices that can be irregular, fragmented, and fractal.

diffusion electrodes by Frumkin et al.3 and the study of a fully electrolyte-impregnated porous electrode in alkaline batteries by Coleman.4 Assuming a constant reactant concentration and a linear polarization, Euler and Nonnenmacher investigated the current density distribution in a porous electrode in 1960.5 In ref 5, Euler and Nonnenmacher proposed the idea of transmission line analysis, which was adopted in the wellknown De Levie model published in 1963.6 In 1962, Newman and Tobias amended the model of Euler and Nonnenmacher by extending the kinetic polarization from linear to nonlinear regime using the Butler−Volmer equation and by considering the effect of variations in reactant concentrations on the kinetics of interfacial reactions.7 In 1973, Newman published the first version of the famous treatise entitled “electrochemical systems”, in which the methodology of modeling porous electrodes is now termed the Newman method.8 Battery models based on the Newman method have been serving as the standard approach since 1990s.8−10 The Russian school on porous electrode theory, represented by Chizmadzhev and Chirkov, has not been well-recognized yet. Their work, summarized in ref 11, focuses on the interplay of oxygen diffusion through porous networks and electrocatalytic reactions at the interface. Recent remarkable progresses are recapitulated below. Lai and Ciucci12 and the Bazant group13 reformulated the porous electrode theory in which the electrochemical potential is taken as the driving force of nonequilibrium thermodynamics in porous electrodes. Several groups have introduced the homogenization theory, which replaces the volume averaging method employed in the Newman method, to upscale the

1. INTRODUCTION In 1894, the Nobel laureate F. W. Ostwald stunningly asserted that electrochemistry is the only solution to the biggest challenge of his time, namely, securing affordable and clean energy.1 More than one century later, today, we are still looking for paths toward sustainable energy.2 Parallel and complementary electrochemical paths include capacitors, batteries, and fuel cells.2 A commonality they all share is the use of porous electrodes, as schematically portrayed in Figure 1. Being porous means that active materials have varied

Figure 1. Schematic illustration of a porous electrode.

distances away from the bulk metal and solution in the electrode. A question then arises: How efficiently are active materials utilized at different locations of the porous electrode? Theory and simulation have been playing a vital role in addressing this question. Theoretical interests in understanding the interplay between mass transport through the porous network and electrochemical processes at the electrode−electrolyte interface date back to 1940s, as exemplified by the study of porous gas © XXXX American Chemical Society

Received: October 9, 2017 Revised: December 10, 2017

A

DOI: 10.1021/acs.jpcc.7b09986 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

characteristic can be observed on scales only in the range of R0 ≪ R ≪ Rm, where R0 is the microscopic scale of constituent particles at the lower end of the spectrum of scales and Rm is the macroscopic scale of the porous electrode at the upper end of the spectrum of scales. The quintessential parameter of a fractal is the fractional dimensionality. Hausdorff dimension provides a rigorous definition of the dimension of an object, i.e., a line, a surface, or a space.42 From the engineering perspective, the dimension of an object is widely calculated by using the box-counting method.43 The basic idea is to cover the investigated object using a set of spherical balls of which the diameter is no larger than ϵ. Therefore, given a specific ϵ, one needs at least Nϵ balls to cover the object. The box-counting dimension of the investigated object is defined as

description of ion transport and reactions at the microscale to the macroscale.14−19 Conventional theories usually assume electroneutrality in the solution phase. By accounting for electrostatic interactions, the Eikerling group20−22 and the Bazant group14,23 incorporated surface charging effects into the porous electrode theory. Specifically, Eikerling et al. developed a theoretical framework for surface charging effects in electrochemical systems; the model self-consistently solves for the surface charge density as a function of electrode potential and solution properties.22 At the microscopic scale (material scale), phase-field models have been employed by the Ceder group24 and the Bazant group13,25 to describe phase transformations and microstructural evolutions in battery materials. Diffusion-induced stresses have been incorporated into the porous electrode theory, as exemplified by recent works by Ji and Guo26 and Wu and Lu.27 Making use of the volume averaging method, previous models captured essential features of a porous electrode without going into exact multiscale geometrical details, and they describe mass transport and reactions in integer dimensional space (Euclidean geometry) using ordinary calculus. This has become the standard approach since the 1970s.8 However, these models may be insufficient to describe porous electrodes that are irregular, fragmented, and fractal, viz., in forms that Euclidean geometry leaves aside. Fractal geometry, coined by Mandelbrot,28 treats objects with a noninteger dimension. The presented work aims at describing mass transport and reactions in fractal media using fractional calculus and thereby generalizing the porous electrode theory for noninteger dimensional space. Treating electrochemical processes in fractal media has been popular since the 1980s.29−39 The interests were partially aroused by experimental observation of the frequency dispersion phenomena, that is, the impedance obeys the relation Z ∝ (jω)−n with n < 1 for double-layer charging and n ≠ 0.5 for diffusion process. As a result, a lot of attention has been paid to relate the exponential factor n with the fractional dimensionality of a rough surface.30,32,34−36,38,39 Moreover, fractional calculus has been employed to describe diffusion processes in fractal media, so-called anomalous diffusion.33,37 The basic idea is to generalize the continuity equation by incorporating the fractional dimensionality into time derivatives. Interested readers are directed to book chapters written by Pyun and co-workers.40,41 However, a theoretical treatment of coupled diffusion and migration in multicomponent, concentrated electrolytes, as well as coupled faradic and nonfaradic reactions at the electrode−electrolyte interface in fractal porous electrodes is missing, which is the gap the presented paper intends to fill. The rest of this paper is organized as follows: First, the definition and determination methods of the fractional dimensionality and a selective introduction of fractional calculus are given in section 2, laying basis for deriving the generalized porous electrode theory in section 3. In the penultimate section, two rudimentary examples, one in the time-domain and the other in the frequency-domain, are given to demystify the theory, and the usefulness and future development of the presented approach are discussed.

ln Nϵ (1) ln ϵ Mass dimension is another common definition, which describes how the mass of an object scales with its size, that is44 D b = − lim

ϵ→ 0

ij R yz m M Dm = M 0jjj zzz j R0 z (2) k { where M0 is the mass of an elementary particle. Fractional calculus is a well-established subject of integrals and derivatives of arbitrary order. Interested readers are referred to pertinent textbooks.45,46 In what follows, we briefly introduce key definitions and equations that will be used in subsequent model development. This section is included for the purpose of transparency and clarity of the present theory. Definitions of Riemann−Liouville (RL) fractional integrals and derivatives are given below. Let f(x) be a continuous function on the real axis. Let Ω = [a, b] (−∞ < a < b < ∞) be a finite interval on the real axis  . T h e l e f t - s i d e d R L f r a c t i o n a l i n t e g r a l of o r d e r α ∈  (ℜ(α) > 0) is defined as D

(Iaα+f )(x) ≔

1 Γ(α)

∫a

x

f (t ) d t (x − t )1 − α

(3)

and the right-sided RL fractional integral of order α is defined as (Ibα−f )(x) ≔

1 Γ(α)

∫x

b

f (t ) d t (t − x)1 − α

(4)

where Γ(x) is the gamma function, Γ(x + 1) = x Γ(x), and Γ(1/2) = π . Fractional integrals are used to describe integral over objects with a fractional dimensionality. The single-variable integral over a fractal in 1 (1D Cartesian space) with a dimensionality of α is given by45−47

∫ f (x ) d α x = =

2. FRACTIONAL DIMENSIONALITY AND FRACTIONAL CALCULUS Porous electrodes can be pictured as a dense object consisting of pores or holes with a fractal structure. The fractal

π α /2 Γ(α) α (I f )(x) Γ(α /2) π α /2 Γ(α /2)

∫ f (x) x α−1 dx

(5)

α

We could rewrite eq 5 as ∫ f(x) d x = ∫ f(x) dos1(α, x) dx with dos1(α, x) = πα/2xα−1/Γ(α/2), as recommended by Tarasov.47 Tarasov discussed different expressions of dos1(α, x) in ref 47. Then, we have, dαx = dos1(α, x) dx. As a result, B

DOI: 10.1021/acs.jpcc.7b09986 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C dos1(α, x) serves as the numerical factor introduced to transform the element in α-dimensional space, dαx, to the element in 1D space, dx. Equation 5 indicates that an integral over a fractal in noninteger dimensional space can be described using a fractional integral. Considering an arbitrary fractal W with a fractional dimensionality of αi (i = 1, ..., n) in each Cartesian coordinate in n (n = 1, 2, and 3), we can express the integration of a function f(x) over the fractal W as

∫ f (x) dα x= ∫ f (x) dosn(αW , x) dnx W

n

(6)



where x = ∑i = 1 x e1, dnx = ∏i n= 1 dxi, αW = ∑i n= 1 αi, and the density of state, dosn(αW, x), is written as n

dosn(αW , x) =

∏ i=1

π αi /2 αi − 1 xi Γ(αi /2)

(7)

Tarasov interpreted dosn(α, x) as how permitted states of elementary particles constituting the fractal W are packed in the space n.45,46 This interpretation can be elucidated by considering the mass of a spherical fractal, as the mass is proportional to the number of elementary particles packed in the fractal W. In terms of a spherically symmetric function, we have45,47

∫ f (r) d

αW

r=

∫ f (r) dos3(αW , r) dr

Figure 2. Schematic illustration of a fractal with a fractional dimensionality of α and a 2D section of the fractal. The electrode− electrolyte interface at which interfacial reactions occur is also fractal and has a fractional dimensionality of β. In the 2D section, part of the area is occupied by the electrolyte phase. The electrolyte phase part in 2 is fractal and has a fractional dimensionality of γ.

μi = μi0 + kBT ln(ai) + zieϕ (8)

(9)

where μi0 is the standard chemical potential, kB is the Boltzmann constant, T is the temperature, zi is the charge number, e is the elementary electron charge, ϕ is the potential, and ai is the absolute chemical activity ai = γici (10)

α W /2 α W −1

with dos 3 (α W , r) = 2π r /Γ(α W /2). Note that dimensionless variables should be used in eq 8 to avoid imbalance in dimension. We define a dimensionless radius, r̃ = r/R0 with R0 being the radius of microscopic constituent particles, and a dimensionless mass, M̃ = M/M0 with M0 being the mass of microscopic constituent particles. The dimensionless mass of a spherical fractal W with a fractional dimensionality of αW and a dimensionless radius of R̃ W is R̃ α expressed as, M̃ = 2παW/2/Γ(αW/2) ∫ 0 W rαW−1 dr̃ = παW/2R̃ WWΓ3 ̃ ̃ (αW/2 + 1). It is readily seen that M = 4πRW/3 when αW. = 3 As a result, dosn(α, x) can be interpreted as the “density of state” or “number density of elementary particles” of the fractal in the space n .

with ci as the concentration and γi as the activity coefficient. The transport of species i is driven by the gradient of μi. Assuming that the electrochemical system is close to equilibrium, the flux of species i is written as Ni = −Mici∇μi

(11)

where Mi is the mobility. Substituting eq 9 into eq 11 gives

3. DESCRIPTION OF TRANSPORT AND REACTIONS We consider an isotropic porous electrode S in an ndimensional space n (n = 1, 2, and 3). The porous electrode is fully filled with a binary z/z electrolyte, as schematically shown in Figure 2. The electrolyte phase E is fractal and has a dimensionality of α ≤ n. Therefore, the dimensionless mass of electrolyte contained in the porous electrode is given by M̃ e ∼ (R̃ e)α, with R̃ e being the dimensionless radius. On the basis of the definition in eq 1, the dimensionality of the solid phase is equal to n.48 The dimensionality of the electrode−electrolyte interface in space n − 1 is denoted as β, which is generally larger than (n − 1). β is a key parameter related to interfacial reactions in the porous electrode. Given a section of the porous electrode in n − 1, the part occupied by electrolyte phase has a dimensionality of γ < (n − 1), as shown in Figure 2. We begin with describing the ion transport in the electrolyte phase. The electrochemical potential of species i in the electrolyte phase is expressed as

Ni = −Dichem∇ci −

zieDici ∇ϕ kBT

(12)

where

ij ∂ ln γi yz zz Dichem = Dijjj1 + j ∂ ln ci zz{ k

(13)

is the chemical diffusivity considering the concentrationdependent activity coefficient, and Di = MikBT is the dilutesolution limit diffusivity. 3.1. Fractional Conservation Law. Let us consider an arbitrary region W in the electrolyte phase E. The boundary ∂W is composed of an electrode−electrolyte interface part, ∂Eeei, where interfacial reactions take place and an electrolyte phase part, ∂Eeps, where ion fluxes flow, as shown in Figure 2. As E and ∂E are assumed to be homogeneous, the dimensionality of W is equal to that of E, α, and the dimensionality of ∂Eeei is equal to that of ∂E, β. The dimensionality of ∂Eeps is given by γ. C

DOI: 10.1021/acs.jpcc.7b09986 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

mass conservation law is obtained when α = n, γ = β = n − 1 in integer dimensional space, namely, ∂ci/∂t = −∇·Ni − aRi. For a spherically symmetrical W, eq 17 can be represented as

The total amount of ions stored in W is given by the integral, ∫ Wci dαx, where dαx denotes α-dimensional volume differential. As a result, the change rate of ions stored in W is written as d(∫ Wci dαx)/dt, which is transformed to an integral in regular space n , ∫ Wd(ci dosn(α, x))/dt dnx, according to eq 6. Chemical reactions in the bulk electrolyte phase are excluded. Therefore, the change of ions stored in W has two sources only. One source is the ion flux, Ni, driven by concentration gradient (diffusion) and potential gradient (migration) through the ∂Eeps, which is given by, ∫ ∂W Ni dγx = ∫ ∂EepsNi dγx, with dγx being γ-dimensional area differential, as Ni = 0 at ∂W except ∂Eeps, ∂W\∂Eeps. Note that dγx may be replaced by dSγ in the literature.45 According to fractional integral, we have, ∫ ∂EepsNi dγx = ∫ ∂EepsNi(dos(n−1)(γ, x) d(n−1)x) As we have transformed the fractional integral over ∂Eeps in γdimensional space to a regular integral in (n − 1)-dimensional space, we can further employ the Gauss’s law, well-established in integer dimensional space, to transform a integral over the boundary in the (n − 1)-dimensional space to a integral over the object in the n−dimensional space. Gauss’s law leads to

∫∂E

∇·(Ni dos(n − 1)(γ , r)) dos(n − 1)(β , r) ∂ci =− −a Ri ∂t dosn(α , r) dosn(α , r) (18)

In the rest of this paper, equations expressed in the Cartesian coordinate can be transformed to its counterpart in the polar coordinate in the similar manner. 3.2. Mass Transport. Applying eq 12 as well as the relation c+ = c− = ce for cations and anions leads to N+ = −D+chem∇ce −

N − = −D−chem∇ce −

(14)

Je = z+F N+ + z −F N −

The other source is electrochemical reactions (including both charge transfer reactions and double layer charging), Ri, at ∂Eeei, which is expressed as, ∫ ∂WRi dβx = ∫ ∂EeeiRi dβx, with dβx or dSβ being the β-dimensional area differential. Ri is positive when ions are consumed in electrochemical reactions. We assume that W is small enough that the electrochemical reaction rate is homogeneous; therefore

∫∂E

R i d βx =

eei

∫∂E

ij z 2FeD c z 2FeD−ce yzz zz∇ϕ Je = −jjjj + + e + − kBT z{ k kBT − (z+D+chem + z −D−chem)F ∇ce

∫W a dos(n−1)(β , x)R i dnx

(15)

σe =



z+2FeD+ce z 2FeD−ce + − kBT kBT

D±chem = z+D+chem + z −D−chem Je = −σe∇ϕ − D±chemF ∇ce

∫W a dos(n−1)(β , x)R i dnx

Dchem ±

−a

dos(n − 1)(β , x) dosn(α , x)

z+Dchem +

(26)

z−Dchem −

Dchem +

Dchem −

when = + = 0, namely, = as z+ = −z− for a binary z/z electrolyte. ∇ce can be expressed as a function of N+ and N− by eliminating ∇ϕ via algebraic calculations of eqs 19and 20

(16)

Equation 16 is valid for an arbitrary region W. Therefore, we can drop the integral operation and obtain

dosn(α , x)

(25)

and the specific relation Je = −σe∇ϕ

∂ci =− ∂t

(24)

Then, we have the general relation

∫W ∇·(Ni dos(n−1)(γ , x)) dnx

∇·(Ni dos(n − 1)(γ , x))

(23)

and a differential diffusivity between cations and anions of the electrolyte as

d (ci dosn(α , x)) dnx dt

=−

(22)

by substituting eqs 19and 20 into eq 21. Let us define the observed conductivity of the electrolytic solution as

R i dos(n − 1)(β , x) d(n − 1)x

where a is termed as the volumetric electrochemical surface area, defined as, a = ∫ ∂W d(n−1)x/∫ W dnx and a = 3/r for a solid spherical ball with a radius of r. On the basis of the conservation law, the increasing rate of concentration of ions stored in W is equal to the sum of ion flux entering W through ∂Eeps and production rate of ions in interfacial reactions at ∂Eeei. Therefore, we write

∫W

(21)

which is expanded to

eei

=

(20)

for anions, respectively. Chemical diffusion coefficients, are used in the diffusion term, while dilute-solution limit diffusion coefficients, Di, are used in the migration term. The current density in the electrolyte phase is carried by both cations and anions, which is defined as

x)

∫W ∇(Ni dos(n−1)(γ , x)) dnx

z −eD−ce ∇ϕ kBT

Dchem , i

eps

=

(19)

for cations and

(n − 1)

Ni(dos(n − 1)(γ , x) d

z+eD+ce ∇ϕ kBT

∇ce =

Ri

z −D−N+ − z+D+N − z+D+D−chem − z −D−D+chem

(27)

which is transformed to

(17)

∇ce =

which is the fractional conservation law for ion concentrations in fractal media. It is readily seen from eq 17 that the ordinary D

−t −N+ − t+N − Damb

(28) DOI: 10.1021/acs.jpcc.7b09986 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

of which the reaction rate Rct, positive-defined for the reduction reaction, is described using the Butler−Volmer equation:

given the following definitions: ambipolar diffusivity Damb: Damb =

z+D+D−chem − z −D−D+chem z+D+ − z −D−

surf ji C z+ i αFη zy zz R ct = k0(C Mb z+)(1 − α)(CMb)α jjjj Mb expjjj− k RT { k C M z+ i (1 − α)Fη yzyzz zzzz − expjjjj zz RT k {{

(29)

transference numbers: t± =

±z±D± z+D+ − z −D−

(30)

Combining eqs 21 and 28 leads to N+ = −Damb∇ce −

where k0 is the standard rate constant, CbMz+ and CbM are the concentrations of metal cations in bulk solution and bulk metal, respectively, Csurf Mz+ is the concentration of metal cations at the reacting surface, α is the transfer coefficient, and η is the overpotential:

t+ Je z −F

(31)

Substituting eq 31 into eq 17 for cations gives ∇·(dos(n − 1)(γ , x) Damb∇ce) ∂ce = ∂t dosn(α , x) +

t+ Je

(

∇· dos(n − 1)(γ , x) z dosn(α , x)

−F

η = ϕM − ϕesurf − Eeq

) − a dos

(n − 1)(β ,

x)

dosn(α , x)

Equation 32 describes ion transport in a multicomponent, concentrated electrolyte with coupled diffusion and migration in an arbitrary dimensional space. Moreover, an equation for charge transport is needed to close the theory, as given below. 3.3. Charge Transport. According to the charge conservation law, the accumulation/consumption rate of charge carried by ion flux through ∂Eeps in the fractal W must be compensated by the accumulation/consumption rate of charge involved in electrochemical reactions at the electrode−electrolyte interface, ∂Eeei. Charge involved in the electrochemical reactions is calculated as a fractional integral over ∂W, given by

R dl =

Je dγ x=

eps

∫∂E

∇·(dos(n − 1)(γ , x) Damb∇ce) ∂ce = ∂t dosn(α , x)

n

i 1 ∇·jjjdos(n − 1)(γ , x) z + + + k dosn(α , x)

(33)

(

dos(n − 1)(β , x)

= −az+FR+

Je

F

(41)

(42)

and

Jy i ∂ce = ∇·(Damb∇ce) + ∇·jjjj(1 − t+) e zzzz ∂t F{ k

(35)

3.4. Electrochemical Reactions. The term R+ in eqs 32 and 35 should be expressed as a function of the potential difference across the EEI and the ionic concentration at the EEI. In what follows, we consider the following electroplating reaction as a typical charge transfer reaction M z+ + z+e → M

) yzzz{

ijij 1 t yz J yz ∂ce = ∇·(Damb∇ce) + ∇·jjjjjjj + + zzz e zzzz j ∂t z − z{ F { kk z+

(34)

Combined, the charge conservation in W is described as ∇·(Je dos(n − 1)(γ , x))

t+ z−

Regarding a regular porous electrode in integer dimensional space, i.e., α = n, β = γ = n − 1, eq 41 is simplified to the wellknown equation

Je (dos(n − 1)(γ , x) d(n − 1)x)

∫W ∇·(Je dos(n−1)(γ , x)) dnx

(40)

3.5. Summary and Nondimensionalization of Derived Equations. Substituting eq 35 into eq 32 leads to

eps

=

(39)

R+ = R ct + R dl

In deriving eq 33, we have employed the assumption that electrochemical reactions are uniform over ∂Eeei. The integration of ion flux over ∂Eeps in the electrolytic solution is given by

∫∂E

Cdl ∂(ϕ − ϕM) F ∂t

where Cdl is the interfacial capacitance. R+ is given by the sum of eqs 37 and 39:

∫∂W z+FR+ dβx = ∫∂W z+FR+ dos(n−1)(β , x) d(n−1)x ∫W dos(n−1)(β , x) z+FR+ d x

(38)

with ϕsurf being the electrolyte phase potential at the reacting e surface, and Eeq the equilibrium potential of the reaction surf corresponding to CbMz+ and CbM. Note that Csurf are Mz+ and ϕe determined by eqs 32 and 35. The double-layer charging rate is expressed as

R+ (32)

=a

(37)

(43)

for the special case z+ = −z− = 1. Equation 42 shares the same form of the controlling equation describing mass conservation in the Newman method.9,10,12,13 Substituting eq 25 into eq 41 gives

(36) E

DOI: 10.1021/acs.jpcc.7b09986 J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

i i 1 ∇·jjjdos(n − 1)(γ , x) jjjDamb − z + + ∂ce k k = ∂t dosn(α , x)

(

i 1 ∇·jjjdos(n − 1)(γ , x) z + + − k dosn(α , x)

(

t+ z−

)

t+ z−

Article

y yzz chem z zz∇cez ± {

)D

{

y σe ∇ϕzzz F

approach, two average structural descriptors are introduced; they are the porosity ϵ, viz., the volume fraction of electrolyte phase in porous electrodes, and the tortuosity τ. Transport parameters are multiplied by a factor of ϵ/τ to reflect the effects of structure on species transport. As is well-known,49,50 the Newman method has the following disadvantages: (i) The relation between ϵ and τ is ambiguous. (ii) It is difficult to determine intrinsic transport properties due to the obstacle of ambiguities in ϵ and τ. The recently developed homogenization method improves over the volume averaging method in two aspects.14−19 First, effective parameters, including ϵ and τ, can be determined from the microstructure and from solving relevant equations. Second, the accuracy of the volume averaging method is up to the scale δ (the ratio between micro- and macroscales), while that of the homogenization method can be up to the scale δ2. In consequence, the volume averaging method is limited to cases of δ ≪ 1, while the homogenization method can release this limitation. However, the homogenization method has two disadvantages. First, it is usually limited to periodic distribution of microscopic particles, although it can be extended to general cases in principle. Second, the equation set derived from the homogenization method is more complicated than that in the Newman method, thus limiting its application in engineers. The present theory introduces the fractional integral approach to tackle structural complexities, where fractional dimensionalities of the porous electrode are employed as structural descriptors. The structural information is reflected in position-dependent dos terms. Intrinsic transport properties are used in the present theory. Moreover, the equation set derived from the present theory has a simpler form compared with that derived from the homogenization method. In its present form, this theory employs integer time derivatives. The presented theory can be extended in future by employing fractional time derivatives, such as in theories of fractional diffusion (also termed as anomalous diffusion).37 We note that the presented framework can be employed to treat heat transfer in fractal media with necessary modifications. Interested readers are referred to works from Tarasov51 and Strichartz.52

{ (44)

for mass transport. Substituting eq 25 into eq 35 gives ∇·(dos(n − 1)(γ , x) (σe∇ϕ + D±chemF ∇ce)) dos(n − 1)(β , x)

= az+FR+

(45)

for charge conservation. Regarding a regular porous electrode of integer dimension, i.e., α = n, β = γ = n − 1, eq 45 returns to the well-known equation ∇·(σe∇ϕ + D±chemF ∇ce) = az+FR+

in the Newman method, ∇·(σe∇ϕ +

9,10,12,13

D±chemF ∇ce)

(46)

and

= aFR+

(47)

for the special case z+ = −z− = 1 with Dchem = Dchem − Dchem ± + − . For the sake of generality, simplification and dimension balance, we introduce dimensionless variables: cẽ = ϕ̃ =

D±chem ce x , x ̃ = , a ̃ = al , D±̃ chem = , Damb cref l

lR+ tD σeRT ϕF , t ̃ = amb ,ξ= , R+̃ = , 2 2 RT Dambcref l Dambcref F (48)

where cref is the reference concentration, l is the characteristic length of the porous electrode, and ξ is a dimensionless parameter reflecting the contribution of migration relative to diffusion in ionic transport. By this point, eqs 45 and 46 can be represented in dimensionless forms: i i 1 ∇·jjjdos(n − 1)(γ , x̃) jjj1 − z + + ∂cẽ k = k ̃ ∂t dosn(α , x̃)

(

i 1 ∇·jjjdos(n − 1)(γ , x̃) z + + k −ξ dosn(α , x̃)

(

t+ z−

t+ z−

y yzz chem z zz∇cẽ z ±

)D̃

)∇ϕyzzz{̃

{

{

4. SIMPLE EXAMPLES Applications of the present theory to modeling of practical battery electrodes are ongoing and will be reported in future. To demystify the theory, we consider simple examples to elucidate main features of the present theory. Before going into the model simulation, we first discuss the three key structural parameters, α, β, and γ. Generally, we have 3 ≥ α > 2, 3 > β ≥ 2, and 2 ≥ γ > 1. α can be determined from the size−mass relation according to eq 2. The range of 2 < α ≤ 3 is expected because the electrolyte phase constitutes a fraction of the volume of the electrode. A simple example to estimate β is given as below. First, let us cut a 2D section of the porous electrode, as shown in Figure 3. Then, we can calculate the fractional dimensionality of the electrode−electrolyte interface in the 2D section, denoted as αline ≥ 1. Thereby, we have β ≥ αline + 1 ≥ 2. Regarding γ, we have 1 < γ ≤ 2 because the electrolyte phase constitutes a fraction of the 2D section. Figure 3 exhibits a regular case of β = 2, γ = 2 corresponding to a regular electrode−electrolyte interface and a representative case of β > 2, γ < 2 corresponding to an irregular electrode−electrolyte interface. Investigation on the intrinsic relations between α, β and γ is out the scope of this

(49)

for mass conservation and ∇·(dos(n − 1)(γ , x̃)(ξ∇ϕ ̃ + D±̃ chem ∇cẽ )) dos(n − 1)(β , x̃)

= az̃ +R+̃

(50)

for charge conservation, respectively. 3.6. Discussion and Comparison. Equations 49 and 50 constitute the dimensionless controlling equations for coupled diffusion and migration as well as coupled faradic and nonfaradic reactions in a porous electrode with an arbitrary dimensionality which is fully impregnated with a multicomponent, concentrated electrolyte. In the Newman method and its derivatives, the volume averaging method is invoked to upscale microscopic equations to the macroscale and simplify the consideration of structural complexities.9,10,12,13 In this F

DOI: 10.1021/acs.jpcc.7b09986 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C The boundary conditions to close eq 51 are cẽ = 1

(52)

at x̃ = 1 and ∂c ẽ /∂x ̃ = j ̃

Figure 3. Two cases of the electrode−electrolyte interface: (A) regular case and (B) irregular case.

study. As a special circumstance, we tentatively assume that α = 3 − δD, β = 2 + δD, γ = 2 − δD, with δD defined as the dimensionality deficiency of the electrolyte phase in a porous electrode. From the definition, we know, 0 ≤ δD < 1. Other parameters used in model simulation are given in Table 1. Table 1. Model Parameters for the Simple Case variable

meaning

value

l a z+ Damb Cdl σe n

thickness volumetric electrochemical surface area charge number of cations ambipolar diffusivity double layer capacitance ionic conductivity dimension of the space

100 μm 106 m−1 1 10−5 m2s−1 20 μFcm−2 1 S cm−1 3

As a rudimentary and classical case, we consider the description of ionic transport across a porous medium toward a planar electrode at x̃ = 0, where the metal deposition reaction takes place, as schematically illustrated in Figure 4A. Note that we neglect reactions in the porous electrodes, namely, R̃ + = 0 in eqs 49 and 50. In other words, we are treating a space-only fractional diffusion model. Consequently, the controlling equation is simplified to ∇·(dos(n − 1)(γ , x)̃ ∇cẽ ) ∂cẽ = ̃ ∂t dosn(α , x)̃

(53)

at x̃ = 0, with j ̃ = jl/(z+FDambcref) being the dimensionless current density applied at x̃ = 0. Figure 4B shows the simulated dimensionless ionic concentration at x̃ = 0, c̃e(x̃ = 0), as a function of the dimensionless time, t.̃ The dashed and dotted lines describe ionic transport in integer dimensional space with a series of the dimensionless current density j ̃ in the range of 0.05 ≤ j ̃ ≤ 0.15 (every 0.01), corresponding to a series of Damb. The dark thick line is derived from a noninteger dimensional space with α = 2.8, γ = 1.8, and j ̃ = 0.1. By comparing the dark thick line (the fractal case) and the dashed and dotted line of j ̃ = 0.1 (the regular case), we find that c̃e(x̃ = 0) decreases slower for the fractal case in the early stage (t ̃ < 0.2). On the contrary, a decrease in c̃e(x̃ = 0) is greater for the fractal case in subsequent stages, indicating that ionic transport is more impeded in fractal electrodes. In conventional theories, Damb is multiplied by a factor of ϵ/τ to effectively reflect the effect of structure on species transport. As a result, we let the dimensionless current density j ̃ vary in the range of 0.05 ≤ j ̃ ≤ 0.15 (every 0.01), accounting for variations in Damb in an equivalent manner. It is found that the concentration profile of the fractal case is paradigmatically different, which cannot be described using conventional theories formulated in integer dimensional space, regardless of adjusting the diffusion coefficient. Next, we consider the impedance response of the same case in Figure 4A in the frequency domain. In so doing, we first transform eq 51 to the frequency domain using Laplace transformation: jωcẽp =

∇·(dos(n − 1)(γ , x)̃ ∇cẽp) dosn(α , x)̃

(54)

where c̃pe is the perturbation in c̃e due to stimulation of a sinusoidal current density with a frequency of ω at x̃ = 0. The boundary conditions to close eq 54 are

(51)

in 1D (x-coordinate).

Figure 4. (A) Schematic illustration of ionic transport across a fractal medium toward a planar electrode where the metal deposition reaction takes place. (B) Dimensionless ionic concentration at x̃ = 0, c̃e(x̃ = 0), as a function of the dimensionless time, t.̃ The dashed and dotted lines describe ionic transport in integer dimensional space with a series of the dimensionless current density j ̃ corresponding to a series of Damb. The dark thick line is derived from diffusion in noninteger dimensional space with α = 2.8, γ = 1.8. Note that we actually treat a 2D XY section of the fractal porous electrode in numerical simulation, as shown in panel (A). In consequence, we used dosn(α, x̃) = dos2(α − 1, x̃) and dos(n−1)(γ, x̃) = dos1(γ − 1, x̃) in the simulation. G

DOI: 10.1021/acs.jpcc.7b09986 J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C cẽp = 0



(55)

*E-mail: [email protected].

(56)

ORCID

Jun Huang: 0000-0002-1668-5361

at x̃ = 0, meaning that an unit dimensionless current density is applied at x̃ = 0. The dimensionless diffusion impedance is defined as Z̃ = cẽp(x ̃ = 0)

AUTHOR INFORMATION

Corresponding Author

at x̃ = 1, saying that the concentration is unaltered here, and ∂cẽp /∂x ̃ = 1

Article

Notes

The author declares no competing financial interest.



(57)

ACKNOWLEDGMENTS This work is partially supported by the starting fund for new faculty members at Central South University. Professor Jianbo Zhang at Tsinghua University, Professor Michael Eikerling at Simon Fraser University and Dr. Thomas Kadyk at Technische Universität Braunschweig are acknowledged for helpful discussions and critical reading.

which is normalized to RTl/(z+F)2Dambcref. Figure 5 shows the dimensionless diffusion impedance derived from eq 54 for different values of δD. δD = 0



REFERENCES

(1) Ostwald, W. Die Wissenschaftliche Elektrochemie Der Gegenwart Und Die Technische Der Zukunft. Z. Elektrochem. Angew. Phys. Chem. 1894, 1, 81−84. (2) Chu, S.; Cui, Y.; Liu, N. The Path Towards Sustainable Energy. Nat. Mater. 2016, 16, 16−22. (3) Frumkin, A. Distribution of the Corrosion Process Along the Tube Length. Zh. Fiz. Khim. 1949, 23, 1477−1482. (4) Coleman, J. J. Dry Cell Dynamics: The Bobbin. Trans. Electrochem. Soc. 1946, 90, 545−583. (5) Euler, J.; Nonnenmacher, W. Stromverteilung in Porösen Elektroden. Electrochim. Acta 1960, 2, 268−286. (6) de Levie, R. On Porous Electrodes in Electrolyte Solutions. Electrochim. Acta 1963, 8, 751−780. (7) Newman, J. S.; Tobias, C. W. Theoretical Analysis of Current Distribution in Porous Electrodes. J. Electrochem. Soc. 1962, 109, 1183−1191. (8) Newman, J.; Thomas-Alyea, K. E. Electrochemical Systems; John Wiley & Sons: Englewood Cliffs, NJ, 2012. (9) Doyle, M.; Fuller, T. F.; Newman, J. Modeling of Galvanostatic Charge and Discharge of the Lithium/Polymer/Insertion Cell. J. Electrochem. Soc. 1993, 140, 1526−1533. (10) Newman, J.; Tiedemann, W. Porous-Electrode Theory with Battery Applications. AIChE J. 1975, 21, 25−41. (11) Chizmadzhev, Y. A.; Markin, V.; Tarasevich, M.; Chirkov, Y. G. Macrokinetic Processes in Porous Media (Fuel Cells); Nanka Publisher: Moscow, 1971 (in Russian). (12) Lai, W.; Ciucci, F. Mathematical Modeling of Porous Battery ElectrodesRevisit of Newman’s Model. Electrochim. Acta 2011, 56, 4369−4377. (13) Ferguson, T. R.; Bazant, M. Z. Nonequilibrium Thermodynamics of Porous Electrodes. J. Electrochem. Soc. 2012, 159, A1967− A1985. (14) Schmuck, M.; Bazant, M. Homogenization of the Poisson– Nernst–Planck Equations for Ion Transport in Charged Porous Media. SIAM J. Appl. Math. 2015, 75, 1369−1401. (15) Allaire, G.; Brizzi, R.; Dufrêche, J.-F.; Mikelić, A.; Piatnitski, A. Role of Non-Ideality for the Ion Transport in Porous Media: Derivation of the Macroscopic Equations Using Upscaling. Phys. D 2014, 282, 39−60. (16) Allaire, G.; Brizzi, R.; Dufrêche, J.-F.; Mikelić, A.; Piatnitski, A. Ion Transport in Porous Media: Derivation of the Macroscopic Equations Using Upscaling and Properties of the Effective Coefficients. Computat. Geosci. 2013, 17, 479−495. (17) Ciucci, F.; Lai, W. Derivation of Micro/Macro Lithium Battery Models from Homogenization. Transp. Porous Media 2011, 88, 249− 270. (18) Allaire, G.; Mikelić, A.; Piatnitski, A. Homogenization of the Linearized Ionic Transport Equations in Rigid Periodic Porous Media. J. Math. Phys. 2010, 51, 123103.

Figure 5. Dimensionless impedance of bounded diffusion in a fractal with the dimensionality deficiency δD = 0, 0.1, and 0.2, respectively: (A) Nyquist plot; (B) Bode plot of phase angle.

corresponds to the regular bounded diffusion, in which case the impedance structure consists of a 45° line in the highfrequency range and a semicircle in the low-frequency range. Regarding bounded diffusion in noninteger dimensional space, where δD becomes larger than zero, the structure of diffusion impedance changes accordingly. The tilted line in the highfrequency range deviates from a 45° line and becomes steeper when δD becomes larger, which is seen clearer in the Bode plot in Figure 5B. Moreover, the size of the low-frequency semicircle grows with increasing δD, indicating that ionic transport is more impeded in fractal electrodes, which is consistent with time-domain analysis in Figure 4B.

5. CONCLUSION Classical porous electrode theories are formulated in integer dimensional space. In this work, the porous electrode theory is generalized for noninteger dimensional space using fractional calculus, furnishing a theory for describing mass transport and reactions in irregular, fragmented and fractal porous electrodes that the classical theory leaves aside. The present theory invokes geometrical dimensionalities as quintessential parameters of the porous electrode structure, while classical theories use the volume fraction and the tortuosity between which the relation is usually ambiguous. The generalized theory is reduced to the classical theory when the geometrical dimensionalities of a porous electrode are integers. A simple case demystifies the present theory. It is found that ionic transport in fractal electrodes described by the present theory is paradigmatically different from that described using conventional theories regardless of adjusting transport properties. H

DOI: 10.1021/acs.jpcc.7b09986 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C (19) Foster, J. M.; Gully, A.; Liu, H.; Krachkovskiy, S.; Wu, Y.; Schougaard, S. B.; Jiang, M.; Goward, G.; Botton, G. A.; Protas, B. Homogenization Study of the Effects of Cycling on the Electronic Conductivity of Commercial Lithium-Ion Battery Cathodes. J. Phys. Chem. C 2015, 119, 12199−12208. (20) Huang, J.; Zhang, J.; Eikerling, M. Theory of Electrostatic Phenomena in Water-Filled Pt Nanopores. Faraday Discuss. 2016, 193, 427−446. (21) Chan, K.; Eikerling, M. A Pore-Scale Model of Oxygen Reduction in Ionomer-Free Catalyst Layers of Pefcs. J. Electrochem. Soc. 2011, 158, B18−B28. (22) Huang, J.; Malek, A.; Zhang, J.; Eikerling, M. H. NonMonotonic Surface Charging Behavior of Platinum: A Paradigm Change. J. Phys. Chem. C 2016, 120, 13587−13595. (23) Storey, B. D.; Bazant, M. Z. Effects of Electrostatic Correlations on Electrokinetic Phenomena. Phys. Rev. E 2012, 86, 056303. (24) Han, B. C.; Van der Ven, A.; Morgan, D.; Ceder, G. Electrochemical Modeling of Intercalation Processes with Phase Field Models. Electrochim. Acta 2004, 49, 4691−4699. (25) Bai, P.; Cogswell, D. A.; Bazant, M. Z. Suppression of Phase Separation in Lifepo4 Nanoparticles During Battery Discharge. Nano Lett. 2011, 11, 4890−4896. (26) Ji, L.; Guo, Z. Analytical Modeling and Simulation of Porous Electrodes: Li-Ion Distribution and Diffusion-Induced Stress. Acta Mech. Sin. 2017, DOI: 10.1007/s10409-017-0704-5. (27) Wu, B.; Lu, W. A Battery Model That Fully Couples Mechanics and Electrochemistry at Both Particle and Electrode Levels by Incorporation of Particle Interaction. J. Power Sources 2017, 360, 360−372. (28) Mandelbrot, B. B.; Pignoni, R. The Fractal Geometry of Nature; WH Freeman: New York, 1983; Vol. 173. (29) Liu, S. H.; Kaplan, T.; Gray, L. J. Ac Response of Fractal Interfaces. Solid State Ionics 1986, 18-19, 65−71. (30) Kaplan, T.; Gray, L. J.; Liu, S. H. Self-Affine Fractal Model for a Metal-Electrolyte Interface. Phys. Rev. B: Condens. Matter Mater. Phys. 1987, 35, 5379−5381. (31) Sapoval, B.; Gutfraind, R.; Meakin, P.; Keddam, M.; Takenouti, H. Equivalent-Circuit, Scaling, Random-Walk Simulation, and an Experimental Study of Self-Similar Fractal Electrodes and Interfaces. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1993, 48, 3333−3344. (32) Keddam, M.; Takenouti, H. Impedance of Fractal Interfaces: New Data on the Von Koch Model. Electrochim. Acta 1988, 33, 445− 448. (33) Nigmatullin, R. R. The Realization of the Generalized Transfer Equation in a Medium with Fractal Geometry. Phys. Status Solidi B 1986, 133, 425−430. (34) Sapoval, B.; Chazalviel, J. N.; Peyriere, J. Effective Impedance of a Non-Blocking Fractal Electrode. Solid State Ionics 1988, 28-30, 1441−1444. (35) Sapoval, B. Fractal Electrodes and Constant Phase Angle Response: Exact Examples and Counter Examples. Solid State Ionics 1987, 23, 253−259. (36) Nyikos, L.; Pajkossy, T. Fractal Dimension and Fractional Power Frequency-Dependent Impedance of Blocking Electrodes. Electrochim. Acta 1985, 30, 1533−1540. (37) Bisquert, J.; Compte, A. Theory of the Electrochemical Impedance of Anomalous Diffusion. J. Electroanal. Chem. 2001, 499, 112−120. (38) Liu, S. H. Fractal Model for the Ac Response of a Rough Interface. Phys. Rev. Lett. 1985, 55, 529−532. (39) Sapoval, B.; Chazalviel, J. N.; Peyrière, J. Electrical Response of Fractal and Porous Interfaces. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 38, 5867−5887. (40) Pyun, S.-I.; Shin, H.-C.; Lee, J.-W.; Go, J.-Y. Lithium Transport through Electrode with Irregular/Partially Inactive Interfaces. In Electrochemistry of Insertion Materials for Hydrogen and Lithium; Pyun, S.-I., Shin, H.-C., Lee, J.-W., Go, J.-Y., Eds.; Springer: Berlin, Heidelberg, 2012; pp 213−237.

(41) Go, J.-Y.; Pyun, S.-I. Fractal Approach to Rough Surfaces and Interfaces in Electrochemistry. In Modern Aspects of Electrochemistry; Vayenas, C. G., White, R. E., Gamboa-Adelco, M. E., Eds.; Springer: Boston, MA, 2006; pp 167−229. (42) Grassberger, P. On the Hausdorff Dimension of Fractal Attractors. J. Stat. Phys. 1981, 26, 173−179. (43) Liebovitch, L. S.; Toth, T. A Fast Algorithm to Determine Fractal Dimensions by Box Counting. Phys. Lett. A 1989, 141, 386− 390. (44) Perfect, E.; Rasiah, V.; Kay, B. D. Fractal Dimensions of Soil Aggregate-Size Distributions Calculated by Number and Mass. Soil Sci. Soc. Am. J. 1992, 56, 1407−1409. (45) Tarasov, V. E. Fractional Dynamics: Applications of Fractional Calculus to Dynamics of Particles, Fields and Media; Springer Science & Business Media: Berlin, 2011. (46) Miller, K. S.; Ross, B. An Introduction to the Fractional Calculus and Fractional Differential Equations; Wiley: New York, 1993. (47) Tarasov, V. E. Anisotropic Fractal Media by Vector Calculus in Non-Integer Dimensional Space. J. Math. Phys. 2014, 55, 083510. (48) The electrode phase S is the subtraction of the electrolyte phase E from the space n (n = 1, 2, and 3). One needs at least Nϵ(n) ∼ ϵ−n balls of radius ϵ to cover n (n = 1, 2, and 3), and Nϵ(E) ∼ ϵ−α balls of radius ϵ to cover E. As a result, Nϵ(S) ∼ (ϵ−n − ϵ−α) for S. Nϵ(S) ∼ ϵ−n when ϵ → 0 and α ≤ n. On the basis of the box-counting definition, the dimension of the electrode phase is equal to n. (49) Landesfeind, J.; Hattendorff, J.; Ehrl, A.; Wall, W. A.; Gasteiger, H. A. Tortuosity Determination of Battery Electrodes and Separators by Impedance Spectroscopy. J. Electrochem. Soc. 2016, 163, A1373− A1387. (50) Tjaden, B.; Brett, D. J.; Shearing, P. R. Tortuosity in Electrochemical Devices: A Review of Calculation Approaches. Int. Mater. Rev. 2018, 63, 47−67. (51) Tarasov, V. E. Heat Transfer in Fractal Materials. Int. J. Heat Mass Transfer 2016, 93, 427−430. (52) Strichartz, R. S. Analysis on Fractals. Not. Am. Math. Soc. 1999, 46, 1199−1208.

I

DOI: 10.1021/acs.jpcc.7b09986 J. Phys. Chem. C XXXX, XXX, XXX−XXX