Mesoscale Phase-Field Modeling of Charge Transport in

Article pubs.acs.org/JPCC

Mesoscale Phase-Field Modeling of Charge Transport in Nanocomposite Electrodes for Lithium-Ion Batteries Shenyang Hu,* Yulan Li, Kevin M. Rosso, and Maria L. Sushko Pacific Northwest National Laboratory, Richland, Washington 99352, United States

ABSTRACT: A phase-field model is developed to investigate the influence of microstructure, thermodynamic and kinetic properties, and charging conditions on charged particle transport in nanocomposite electrodes. Two sets of field variables are used to describe the microstructure. One is comprised of the order parameters describing size, orientation, and spatial distributions of nanoparticles, and the other is comprised of the concentrations of mobile species. A porous nanoparticle microstructure filled with electrolyte is taken as a model system to test the phase-field model. Inhomogeneous and anisotropic dielectric constants and mobilities of charged particles, and stresses associated with lattice deformation due to Li-ion insertion/ extraction, are considered in the model. Iteration methods are used to find the elastic and electric fields in an elastically and electrically inhomogeneous medium. The results demonstrate that the model is capable of predicting charge separation associated with the formation of a double layer at the electrochemical interface between solid and electrolyte, and the effect of microstructure, inhomogeneous and anisotropic thermodynamic and kinetic properties, charge rates, and stresses on voltage versus current density and capacity during charging and discharging.

1. INTRODUCTION Nanocomposites of carbon and metal oxide nanoparticles comprise a new basis for high performance electrodes in Li ion batteries, essentially because of their substantial advantages in terms of higher electrode−electrolyte contact area and shorter diffusion distance for both electronic and Li ion transport.1−8 Use of small particles appears advantageous also because volume changes that accompany Li+ insertion/extraction that would normally lead to stress degradation of bulk materials are somewhat alleviated in nanoparticles, particularly when pore space is available to accommodate volume changes. Nanostructured TiO2 has attracted significant attention as a low voltage and high conductivity anode material for Li ion batteries. For TiO2, eight polymorphs are known: rutile, anatase, brookite, TiO2 B, TiO2 R, TiO2 H, TiO2 II, and TiO2 II.9 Anatase is generally considered to be the most electroactive Li-insertion host, while Li insertion into rutile is usually reported to be negligible.10,11 However, recent experiments show that nanometer-sized rutile has surprising advantages as compared to micrometer-sized rutile.7 Not only is the rate performance remarkable, but also a high capacity has been achieved at room temperature, in contrast to previous works on Li insertion into micrometer-sized rutile. Nanometer-sized © 2012 American Chemical Society

rutile is able to reversibly accommodate Li up to Li0.5TiO2 (168 mAh g−1) at 1−3 V versus Li+/Li with excellent capacity retention and high rate capability. Furthermore, experiments show that (1) the mobility of Li ions along the c-axis in rutile is much larger than that in the ab-plane; (2) irreversible phase transition from cubic TiO2 to rocksalt-type hexagonal LiTiO2 or amorphous-like structures may take place due to the volume expansion in the ab-plane during Li insertion into the nanometer-sized rutile; and (3) stresses associated with volume expansion could impact charged particle transport kinetics and structure integration. Although many mathematical models have been developed to describe the charged particle transport in fuel cells and batteries,12−19 microstructure, inhomogeneous, and anisotropic structural, thermodynamic, and kinetic properties and stresses are often ignored. These aspects have so far been challenging to incorporate because they are intrinsically multiscale phenomena that can span from the atomic scale to the macroscopic scale. Therefore, previous models have Received: July 9, 2012 Revised: November 30, 2012 Published: December 13, 2012 28

dx.doi.org/10.1021/jp3068014 | J. Phys. Chem. C 2013, 117, 28−40

The Journal of Physical Chemistry C

Article

Figure 1. (a) TEM image of mesoporous TiO2 particle composite,42 and (b) microstructure in a representative volume used in simulations.

nanosized grains. The phase-field model presented here thereby does not rely on parameter guessing. It takes into account the microstructure, anisotropic, thermodynamic, and kinetic properties, and elastic properties of nanostructured materials. We use TiO2 nanocomposites as a model system to demonstrate the capabilities of our approach for quantitative predictions of the electrochemical performance of the material as a function of its microstructure. Nevertheless, the qualitative conclusions presented here are more general and are applicable to similar systems.

limitations in predicting the outcome of such effects on nanocomposite electrode performance. Mesoscale phase-field methods have been successfully applied to predicting the evolution of materials properties in many processes including, for example, diffusion-controlled phase transitions in alloys,20−22 ferroelectric/ferromagnetic phase transitions,23,24 and microstructure evolution in irradiated materials.25,26 Garcia et al. proposed a theoretical framework to derive thermodynamically consistent equilibrium equations, kinetic driving forces, and governing equations for microstructure evolution in electrically and magnetically active materials.27 Han et al. first applied the phase-field model to intercalation process in a Li-ion battery with a regular solution model for LiFePO4.28 Various electrochemical processes such as electrodeposition and electrochemical impedance spectroscopy,29−32 and phase separation in LiFePO4 nanoparticles33,34 have also been simulated using phase-field models. Very recently, intercalation dynamics,35,36 the effect of elastic energy,37 electrode porosity, and phase transformations38 on intercalation dynamics have been incorporated consistently in phase-field models. However, in most cases, not all of the parameters of the models are known from experiment, and certain estimates have to be made, which reduces the predictive power of the approach. One of the solutions is to invoke atomistic simulations to make more informed predictions of thermodynamic properties of the material. Amajor challenge, however, is that thermodynamic properties calculated for bulk materials, that is, using periodic models, cannot be directly applied for the description of nanoparticles and nanocomposites. To overcome this problem, here we use a multiscale model, in which nanosizing effects are fully accounted for. The model uses a multiphysics approach, in which instead of formal consecutive upscaling, we introduce novel types of collective long-range interactions along with short-range effects of finer scale models to describe coupled Li ion and electron transport at the nanoscale.39−41 The collective long-range electrostatic and excluded volume interactions are introduced on the mesoscale (10−300 nm) via classical Density Functional Theory (cDFT) coupled with a Poisson−Nernst−Planck (PNP) formalism for dynamic effects. The mesoscopic free energy, which includes contributions from short-range activation dynamics of ions and electrons, derived from the atomistic models, is then used in a larger scale (micrometers) phase-field model to simulate charge transport in a network of

2. DESCRIPTION OF THE PHASE-FIELD MODEL Porous nanocomposites have complicated structures. Figure 1a shows a porous carbon-based TiO2 nanocomposite consisting of TiO2 particles with different sizes and orientations, carbon, and electrolyte.42 To simplify the modeling problem, some assumptions are made. First, we assume that electrons can quickly reach and leave all of the interfaces between TiO2 particles and electrolyte through the conductive carbon coating, while Li ions can quickly reach and leave the interface through the electrolyte. With this assumption, the microstructure of the carbon coating itself can be ignored. Figure 1b schematically shows the representative volume for a network of TiO2 nanoparticles chosen for our model. Nanoparticles have different sizes and orientations, determined by the orientation of the c-axis of the rutile TiO2 crystal structure. Second, the electrode is considered as a phase (denoted as α phase) of primarily TiO2, electrons, and Li ions. To facilitate the discretization of the governing equations, three sublattices are introduced in the electrode phase: the firstsublattice is the host lattice occupied by TiO2 and vacancies; the second is occupied by Li ion and vacancies, and the third is occupied by electron and vacancies. Third, the electrolyte is considered a phase (denoted as β phase) of primarily three components: Li+, anion A−, and solvent M. We use two sets of phase-field variables to describe the microstructure. One is the set of order parameters ξi (i = 1,2,....m), where m is the total number of particles. The order parameter ξi is 1 inside the ith particle, while it is 0 outside the ith particle, and varies from 1 to 0 across the interface among different particles and electrolyte. All of the thermodynamic and kinetic properties in particles and in the electrolyte are described by the order parameters. The other set of variables is the set of concentrations Xi, i = 1−5. X1, X2, and X3 are the mole fractions of electrons, Li+, and A− per mole TiO2 lattice, respectively. X4 and X5 are the mole fractions of 29



Article

Figure 2. Free energies of electrode TiO2 and electrolyte LiFP6/PC phases versus electron and Li ion mole fractions for different applied electric potentials from classical Density Functional Theory (cDFT) calculations.

TiO2 and solvent M per mole TiO2 lattice, respectively. The molar volume of TiO2 is Ω, which is a function of Li ion concentration if we assume that the solubility of A− in TiO2 and the molar volume of electrons are both zero. Therefore, the concentrations of different species can be calculated by Ci = Xi/ Ω, i = 1−5. Finally, we assume that Li ion insertion and extraction causes a small volume change. With this assumption, the deformation can be described by the linear elasticity theory. In the framework of the phase-field model, the total free energy of an isothermal system of charged particlesis described in terms of phase-field variables ξi and Xj as ⎡

E=

n

5

j=1

i−1

fV (ξj , Xi , T ) =

⎣

+

∑ j=1

⎤ 1 1 2 |∇ξj| + ρφ + σijεij ⎥ dV ⎥⎦ 2 2 2

F Ω

⎧ ∂Xi δE ⎫ ⎬, = ∇·⎨Mi(r) ·∇ ∂t δXi ⎭ ⎩

κj2

(1)

∂ξj ∂t

i=1

(3)

= −L

δE , δξj

i = 1, 2, ...5 (4)

j = 1, 2, ...n (5)

where Mi(r) is the mobility tensor of the ith species, L is the interface mobility or phase-field ξi, and (δF/δYi) is the variation of total free energy with respect to the phase-field Yi. Because charged particle diffusion is a slow process, mechanical equilibrium is assumed to be reached simultaneously: ∂σij

n

∑ ziXi

j=1

where fα(Xi,T) and fβ(Xi,T) are free energy densities of electrode and electrolyte phases, respectively. h(ξi) = ξ3i (6ξ2i − 15ξi + 10) is an interpolating function, which is 0 inside TiO2 particles or electrode, and 1 outside TiO2 particles or in electrolyte. g(ξi) = ξ2i (1 − ξi)2 is a double-well function. w0 is the double-well height. The gradient energy coefficient κj and double well height w0 are determined by the interface thickness and interfacial energy. By minimizing the total free energy of the system, we arrive at the following governing equations, which describe charge transport and microstructure evolution:

where f V is the bulk free energy per unit volume, and T is the absolute temperature. The second term in eq 1 is the potential well at the interface between electrolyte and TiO2 particles. The coefficients wi, i = 1,2,...5, are determined by the energy barriers for changed particles to pass through the interface. The third term is the gradient energy. κj are the gradient energy coefficients for the phase field. The fourth term is the electrostatic energy, where ϕ is the electrostatic potential, and ρ=

n

∑ {[1 − h(ξj)]fα (Xi , T )

+ h(ξj)fβ (Xi , T ) + w0g (ξj)}

∫V ⎢⎢fV (ξ1, ...ξm; X1, ...Xn; T ) + ∑ (∑ wXi i)g(ξj) m

1 Ω

∂rj (2)

= 0,

i , j = 1, 2, 3 (6)

In addition, for the electric potential, Poisson’s equation must be satisfied everywhere:

is the charge density, zi is the valence of species j, and F is Faraday’s constant. The last term is the deformation energy density associated with lattice parameter changes due to Li insertion and extraction during charging and discharging. σij and εij are the stress and the strain tensors, respectively. The bulk free energy density can be constructed as

∇·[ε(r)∇ϕ] + ρ = 0

(7)

where ε(r) is the dielectric tensor. Therefore, charged particle transport and microstructure evolution can be obtained by solving the governing eqs 4 and 5 together with the equilibrium conditions 6 and 7. 30



Article

3. THERMODYNAMIC AND KINETIC PROPERTIES OF THE SYSTEM 3.1. Electrochemical Potentials. If we view electrons and Li ions as solutes in two sublattices in the host TiO2 lattice, the free energy for dilute solution can be simply written as

fβ (Xi , T ) = f β0 (Xi , T , ϕ) +

1 ρ(ϕ1 − ϕ) 2

(10)

where ϕ is the applied voltage, and ϕ0 and ϕ1 are voltages of electrode TiO2 and electrolyte versus Li+/Li, respectively. At equilibrium, the electrochemical potential for any species should be equal in the two phases, that is:

2

1 fα (X1 , X 2 , T ) = [∑ Xiμi0 + RT (Xi ln Xi Ω i=1

∂fα (Xi , ϕ , T ) ∂Xi

1 ρϕ (8) 2 where Xi are the mole fraction of electrons and Li ions, respectively, R is the molar gas constant, T is the temperature, and μi are the chemical potentials of electrons and Li ions, respectively. However, it is known that the concentration X2 of Li ion may become as high as 0.5 in TiO2 during charging. In such a high concentration regime, the chemical potential of the species becomes concentration dependent, and the expression for the ideal solution free energy is no longer valid. Therefore, we used classical Density Functional Theory (cDFT)39,40 to calculate the free energies of TiO2 and LiFP6 in propylene carbonate (PC) in terms of given electric potential and the concentrations of Li ions and electrons (Figure 2). It is clear that the chemical potential strongly depends on concentrations. The equilibrium concentrations in TiO2 and electrolyte phases can be calculated by the common tangent of two free energies. We can find that the equilibrium concentrations of electron and Li ion in electrolyte are almost independent of applied electric potential, while the equilibrium concentrations of Li ions and electrons in TiO2 increase as applied electric potential increases. Figure 3 shows typical discharge−charge capacity + (1 − Xi) ln(1 − Xi)] +

= X i = X iα *

∂fα (Xi , ϕ , T ) ∂Xi

X i = X iβ *

(11)

On the basis of equilibrium concentrations from the cDFT calculations and experiments, the functions f0α(Xi,T,ϕ) and f0β(Xi,T,ϕ) can be determined. Figure 4 shows the phase diagram calculated from the fitted free energy functions, which were used in the following simulations.

Figure 4. Potential−composition phase diagram used in the simulations, illustrating the bulk equilibrium between electrode TiO2 and electrolyte phases. The lines denote different values of potentials (ϕ0 − ϕ). ϕ0 = 3 V is the voltage of TiO2 versus Li+/Li. ϕ is the applied voltage. The electrolyte phase has the equilibrium mole = 0, Xα,eq = Xα,eq = 0.1, Xα,eq = 0, and Xα,eq = 1.0, which fraction: Xα,eq 1 2 3 4 5 are independent of applied voltage. The electrode phase has the = 0, Xβ,eq = 1.0, and Xβ,eq = 0, while equilibrium mole fraction: Xβ,eq 3 4 5 β,eq β,eq X1 = X2 depends on the applied voltage ϕ.

3.2. Elastic Energy. As can be deduced from Figure 1, the porous nanocomposite is an elastically inhomogeneous system. The electrode TiO2 phase and electrolyte LiFP6/PC phase have different elastic constants. Although TiO2 rutile crystals have anisotropic elastic properties,43,44 we assume that the TiO2 particles are elastically isotropic. C44 = 124 GPa, C11/C44 = 3, and C12/C44 = 1 are used.44 We also assumed that the electrolyte does not sustain stress, and thereby the corresponding elastic constants are zero. Li ion insertion or extraction will cause the lattice constant change. If we assume that the variation of stress-free lattice parameter of TiO2, a, with the given concentrations of Li ion obeys Vegard’s law, the local stress-free strain caused by the defect inhomogeneity is given by

Figure 3. Discharge−charge capacity profile of mesoporous crystalline TiO2 at various rates (1C − C/10) between voltage limits of 1 and 3 V (from ref 1).

profiles in terms of cycling rates in a mesoporous TiO2 electrode. It can be seen that the capacity profile converges when the rate decreases. We assume that the capacity profile is controlled by thermodynamic properties when the rate is low. In other words, the two phases electrode and electrolyte at the interface reach dynamic equilibrium when the rate is less than a critical value. For instance, it is 10/C for the case shown in Figure 3. Therefore, the equilibrium concentrations of Li ions in the electrode for a given voltage can be calculated from the converged capacity profile. We propose general free energy formulas of electrode and electrolyte phases as 1 fα (Xi , T ) = f α0 (Xi , T , ϕ) + ρ(ϕ0 − ϕ) (9) 2

εij* = ε 0(X 2 − X 2α ,eq )δij

(12)

where ε = (1/a) da/dX2 is the expansion coefficient of the lattice parameter due to the insertion of Li ion, and δij is the Kronecker-delta function. For a multiparticle system, the total stress-free strain tensor is 0

31



Article

height w0) and interface properties (the interfacial energy γ and interface thickness 2λ0) can be derived, giving

n

εijtot * =

∑ ε0(X 2 − X 2α ,eq)δijh(ξm) m=1

The elastic energy density U

U def =

def

(13)

is calculated by:

1 λijklεijelεklel 2

γ= (14)

(15)

(16)

D loc

where ε(r) is dielectric tensor. Because the dielectric constants in TiO2 and in the electrolyte are different, the system is dielectrically inhomogeneous. In addition, for single-crystal TiO2 particles, the dielectric constants are anisotropic; for instance, εc = 175ε0, εa = εb = 50ε0 in TiO2 rutile crystals.46 The subscripts a,b,c denote the crystallographic directions. Also, ε0 is the permittivity of free space. The dielectric constant of propylene carbonate solvent is 65ε0. If we consider single-crystal TiO2 particles, in the global coordinate system the dielectric tensor can be calculated as ε(r) = AεlocAT

εloc

κj w0

(21)

⎡ Da 0 0 ⎤ ⎢ ⎥ = ⎢ 0 Db 0 ⎥ ⎢ ⎥ ⎣ 0 0 Dc ⎦

(22)

Also, in the global coordinate system, the diffusivity can be calculated using the rotation matrix 19 as

D(r) = AD locAT

(23)

In eq 4, Mi(r) is the mobility of the ith species. It is related to the diffusivity Di(r) by Mi(r) = (1 − Xi)XiDi /RT

(24)

Figure 5 shows a plot of measured dQ/dV versus potential of lithiated/delithiated mesoporous crystalline TiO2. It is

(17)

where εloc is the dielectric tensor in the local coordinate system. It is

⎡ εa 0 0 ⎤ ⎢ ⎥ = ⎢ 0 εb 0 ⎥ ⎢ ⎥ ⎣ 0 0 εc ⎦

(20)

where α is a dimensionless coefficient equal to ∼2.2. In an electrochemical system, negative and positive particles may separate at the interface and form a double layer. The coefficient w0 gives the height of the energy barriers at the interface, which will affect the interfacial energy as well as the double layer structure. In the present work, they are chosen to reproduce a double layer structure obtained using a onedimensional PNP/cDFT model,33 which simulates ion and electron distributions at the interface from first principles. 3.5. Anisotropic Mobility of Species. Experimental and theoretical studies showed that the Li ion diffusion coefficient along the c-axis is approximately Dc = 10−10 m2/s, while it is only about Da = Db = Dab = 10−19 m2/s in the ab-plane.1 Therefore, the local diffusivity matrix is

where ε̅ij is the homogeneous strain characterizing the macroscopic shape and volume change, and δεij(r) is the heterogeneous strain of ∫ Vδεij(r) dV = 0. For elastic inhomogeneous solids, the elastic solution ((∂σij/∂rj = 0, j = 1,2,3) can be obtained using iteration methods.45 3.3. Electrostatic Energy. For a given charge distribution and applied voltage, the potential should satisfy Poisson’s equation: ∇·[ε(r)∇ϕ] + ρ = 0

3 2

2λ 0 = α 2

where the summation convention over the repeated indexes is used. λijkl = ∑nm = 1λαijkl(1 − h(ξm)) + λβijklh(ξm) is the elastic constant tensor of the system (its Voigt notation is Cij), and λαijkl and λβijkl are elastic constant tensors of electrode and electrolyte phases, respectively. εelij is the elastic strain, which is calculated by εijel = εij̅ + δεij(r) − εijtot *(r)

κj w0

(18)

A is a rotation matrix with components aij: ⎡ a11 a12 a13 ⎤ ⎢ ⎥ A = ⎢ a 21 a 22 a 23 ⎥ ⎢⎣ a31 a32 a33 ⎥⎦ ⎡ cos ψ cos ϕ cos θ − cos ψ cos ϕ cos θ sin ψ cos ϕ ⎤ ⎢ ⎥ − sin ϕ sin θ − sin ϕ cos θ ⎢ ⎥ ⎢ ⎥ = ⎢ cos ψ cos ϕ cos θ − cos ψ sin ϕ sin θ sin ψ sin ϕ ⎥ ⎢ ⎥ + cos ϕ sin θ + cos ϕ cos θ ⎢ ⎥ ⎢⎣ − sin ψ cos θ sin ψ sin θ cos ψ ⎥⎦

Figure 5. dQ/dV versus potential plot of lithiated/delithiated mesoporous crystalline TiO2 (from ref 1).

commonly observed in both bulk and nanostructured materials that the current has a peak at around 1.7 V during charging (increase in Li+ concentration in the anode) and around 1.9 V during discharging. From a thermodynamic viewpoint, the high current or high charged particle fluxes could be due to a high driving force and/or high mobility at specific voltages. The voltage dependence of the driving force is described by the

(19)

and AT is the transpose of A. 3.4. Interfacial Energy. Using equilibrium solutions of eqs 4 and 5, the relationships between phase-field model parameters (gradient energy coefficient κj and double well 32



Article

electrochemical potentials. The mobility may also depend on the voltage and Li concentration, especially when the Li concentration is high, which causes high pressure and reduce the mobility. Without a loss of generality, the effective diffusivity is assumed to depend on voltage, and experimental data from Figure 5 are used to introduce such dependence in Di (Figure 6). In the simulations, diffusivity of electrons in TiO2 is

Table 1. Model Parameters, and Thermodynamic and Kinetic Properties

Figure 6. Voltage dependence of diffusivity. All of the diffusivities Di of species in TiO2 are scaled by the factor D0 with D̅ i = D0Di.

5. RESULTS AND DISCUSSION 5.1. One-Dimensional Simulations. First, we evaluate the phase-field model by simulating the charging and discharging processes in one dimension with different rates dϕ*/dt*. The simulations start in a simulation cell 256 dr*1 × 1 dr*2 × 1 dr*3 with a TiO2 particle at the center of the cell. The particle has a radius R = 64 dr1*. The initial concentrations of different species are X1 = 0, X2 = X3 = 0.1, X4 = 0, and X5 = 1.0 in the electrolyte and X1 = X2 = 0.001, X3 = 0, X4 = 1, and X5 = 0 in the TiO2 particle. The initial voltage in the TiO2 particle is ϕ0 = 3 V, while ϕ1 = 1.7 V in electrolyte. During the simulation, the concentration of Li ion at the electrolyte side of the interface between the particle and electrolyte (order parameter: ξ = 0.09−0.99) is checked and updated to equilibrium concentration X2 = 0.1. The same amount of electrons is added or removed at the interface ξ = 0.09−0.99 to keep charge neutrality in the whole system. Simulations of the evolution of the average Li capacity in the TiO2 particle during discharging and charging with different charging rates, dϕ*/dt*, show good agreement with the experimental data presented in Figure 2 (Figure 7). As expected, for higher charge rate, Li capacity cannot reach saturation, while Li cannot be fully extracted in a higher discharge rate regime. In contrast, the discharging and charging processes become fully reversible with the decrease in charging rate. The charged particle concentrations are tracked during discharging and charging processes. Because of the symmetry of all initial conditions and the periodic boundary conditions with respect to the center of the simulation cell, it is sufficient to plot charged particle concentrations and voltage only in one-half of the simulation cell (Figures 8 and 9). Figure 8 shows the evolution of electron and Li ion concentrations, and voltage for the case of a high charging rate dϕ*/dt* = 0.0144. It is clear that the electron and Li ion concentrations near the interface increase and approach the equilibrium concentration Xα,eq = 1 Xα,eq = 0.4 (described by the free energy and phase diagram) in 2 the TiO2 phase. Inside the TiO2 particle, the concentrations slowly increase through correlated Li ion and electron hopping diffusion driven by the chemical free energy, evaluated at the

assumed to be the same as that for Li ion. In electrolyte, the diffusivity of electron, Li ion, and anion A− are also assumed to be equal to Dc.

4. NUMERICAL METHODS We use the following normalizations for our numerical calculations: t* = (C44Mc/l20)t = t/t0, t0 = l20/(C44Mc), ri* = ri/ l0, L* = Ll20/Mc, f * = f/C44, w0* = w0/C44, wi* = wi/C44, κj* = κ2j / (l20C44),ϕ* = ϕ/ϕ̅ , σ*ij = σij/C44 with l0 being a characteristic length and Dc being the diffusivity of Li ions along the c-axis of TiO2 rutile. Mc = DcΩ/(T9 ) is the mobility of Li ions along the c-axis in rutile. 9 is the ideal gas constant, C44 is the shear modulus of the considered material, ϕ̅ = 1 V, and T = 300 K. By solving eqs 4 and 5 together with mechanical equilibrium eq 6 and Poisson eq 7, we can obtain the temporal and spatial distribution of charged particle mole fraction Xi(r*,t*), electric potential ϕ*(r*,t*), and stresses σ*ij (r*,t*). We use the semiimplicit Fourier spectral method for the time-stepping and spatial discretization for evolution of eqs 4 and 5, and iteration methods47to solve the mechanical equilibrium eq 6 and Poisson eq 7 in this elastically and dielectrically inhomogeneous and anisotropic system. Table 1 lists the normalized model parameters and thermodynamic and kinetic properties. In our simulations, the applied voltage is used to control the charging and discharging processes.Voltage with a rate dϕ*/dt* is continuously applied to mimic charging (dϕ*/dt* > 0) and discharging (dϕ*/dt* < 0).Without spatial charge distribution, that is, ρ = 0, the voltage is ϕ1 in the electrolyte, while the voltage is ϕ*0 − ϕ* in TiO2. The formation of a double layer at the interface produces an additional voltage variation in the system. In addition, we assume that the diffusion of species in the electrolyte is quick so that the concentrations at the electrolyte side of the interface between TiO2 particles and electrolyte always reach equilibrium concentrations. The equilibrium condition at the interface is used to determine the fluxes of Li ion and electrons at the interface. 33

description

value

dt* l0 dr*1 = dr*2 = dr*3 κ*j w*0 w*1 = w*2 w*3 ,w*4 ,w*5 Ω l0 ϕ1 T Dc Da/Dc εc εa/εc εelectrolyte interfacial energy elastic constants ε0

0.0001 1.5 nm 1 0.0024 0.0012 0.5 0.0 3.76 × 10−5 m3/mol 3V 1.7 V 300 K 10−10 m2/s 1, 0.1, and 0.001 176ε0 1 and 0.37 50ε0 1.0 J/m2 C11/C44 = 3, C12/C44 = 1, C44 = 115 GPa 0.0, 0.05, and 0.08



Article

Li concentration does not drop to zero in the center of the nanoparticle after discharging. In the electrolyte phase, the electron and Li ion concentrations retain the equilibrium values Xα,eq = 0.0 and 1 = 0.1 as required by the phase diagram in Figure 2. The Xα,eq 2 right side panels show the voltage distribution evolution; it can be seen that the voltages inside the TiO2 particle and the electrolyte are almost constant. The variation of voltage near the interface is mainly due to charged particle separation, that is, double layer formation at the interface. Figure 9 shows the concentration evolution with a small charging rate dϕ*/dt* = 9.00 × 10−4. As compared to the results in Figure 8, the main difference is that the electron and Li ion concentration distribution inside TiO2 particles become much more uniform with the decrease in charging rate. This implies that the system reaches dynamic equilibrium. As a consequence, the discharging and charging processes are almost reversible (Figure 7). A weak segregation of anions A− at the electrolyte side of the interface is observed (Figure 8c). In contrast, electrons segregate at the TiO2 particle side of the interface (Figure 8b). The total fluxes of Li ion across the interface are plotted in Figure 10. dQ*/dt* is the amount of Li ion passing through the

Figure 7. Voltage versus capacity curves predicted by the simulations with different charge rates.

finer scale. As mentioned earlier, in the high charge rate regime, the discharging and charging processes are not reversible, and

Figure 8. Evolution of electron and Li ion concentrations and voltage (a) during charging and (b) during discharging with charge and discharge rates dϕ*/dt* = ±0.0144. The origin of the x-axis denotes the center of the TiO2 particle; 95−100 nm region denotes the interface, and 100−195 nm region is the electrolyte phase. 34



Article

Figure 9. Evolution of electron, Li ion, anion A−, and charge during charge with a small charge rate dϕ*/dt* = 9.00 × 10−4. The same notations as in Figure 8 are used.

voltage, which depends on the discharging and charging rates. It is interesting to find that the peak voltage increases with the decrease in charging rate, while it decreases with the decrease in discharging rate. 5.2. 2D Results. To examine the effects of electrode microstructure, mobility anisotropy, and elastic interaction on charge transport, we performed simulations in two and three dimensions. In this section, the phase-field model is used to simulate charge transport in a 2D representative volume as shown in Figure 1b. Cases were designed such that particles could have different sizes and orientations, anisotropic Li ion mobilities, and volume changes associated with Li ion insertion and extraction. The initial conditions, interface equilibrium assumption, and discharging and charging processes are the same as those described in section 5.1. 5.2.1. Effect of Diffusivity Anisotropy on Charged Particle Transport. We designed three different microstructures in the representative volume. In Case 1, three particles have either amorphous or polycrystalline structure, and the diffusivity in the particles is assumed to be isotropic (Da/Dc = 1). In Cases 2 and 3, the particles are assumed to have single-crystal structures and anisotropic diffusivities. In Case 2, the c-axes of the three TiO2 particles align in the same direction, and in Case 3, the c-

Figure 10. Current voltage curves for TiO2 particles in electrolyte for charging rates ranging from 0.0009 to 0.0144.

interface per unit time t* and is related to the current. It is clear that the model reproduces the peak of current at a specific 35



Article

All of these results can be consistently explained using a space charge model, which dictates charge accumulation at the boundaries of nanoparticles with diameters larger than the Debye screening length in the material.39,40 This space charge layer screens further Li ion and electron insertion and limits the capacity of larger nanoparticles. Thereby charge is preferentially accumulated in the smaller nanoparticles. The whole charging and discharging process is also simulated for Case 3 with a slower charging rate dϕ*/dt* = 0.009 (Figures 13 and 14). A slower charging rate promotes a larger

axes of the three particles are oriented at the angles of 120°, 45°, and 0° relative to the x-axis. Figure 11 shows the voltage versus capacity curves during charging with a charging rate dϕ*/dt* = 0.144 for three

Figure 11. Voltage versus capacity curves during charging with a rate dϕ*/dt* = 0.144 for three different structures. Case 1, isotropic diffusivity; Case 2, c-axis of three particles aligns in the same direction; and Case 3, three particles have different orientations as shown in the insets.

different microstructures. In these simulations, we assumed that the diffusivity of electron and Li ions in ab-planes is two magnitudes smaller than that along the c-axis, that is, Da/Dc = 0.01. Simulations show that reduced diffusivity in ab-plane slows the charge process. For example, in Case 2, diffusion between nanoparticles has to occur through the ab-planes, and the overall diffusivity through the network is smaller than that in Case 3. The electron distributions at points A, B, and C in the voltage versus capacity curves are presented in Figure 12. Li ion distributions are almost the same as those of electrons to retain charge neutrality. From the distributions, we can conclude that (1) the concentration near particle surfaces is much higher than that in their centers; (2) the anisotropy in diffusivities and in particle orientations causes a more inhomogeneous charged particle distribution; (3) the concentration near particle surfaces and along fast diffusion channels is much higher than that in other directions; and (4) there is higher electron (and Li ion) concentrations in smaller particles.

Figure 13. Voltage versus capacity curve during charging and discharging with a charging rate dϕ*/dt* = 0.009 for the microstructure of Case 3.

capacity of the network of nanoparticles: the average mole fraction of Li ion is about 0.3 for dϕ*/dt* = 0.009, while it is 0.19 for dϕ*/dt* = 0.144. With the charging rate dϕ*/dt* = 0.009, the charging and discharging processes are irreversible. There is about 0.1 mol fraction of Li left inside the particles when the voltage comes back to the original 3 V. Figure 14 plots the temporal evolution of electron and Li ion concentrations at points A−E during charging and discharging. Similar to 1D simulations, it can be seen that the model predicts almost identical distributions of electrons and Li ions, as dictated by the condition of detailed balance between Li ion and electron concentrations imposed at the interface for electroneutrality. The distributions of Li ion and electrons are

Figure 12. Electron distribution at points A, B, and C shown in Figure 11. Li ion distributions are almost the same as those of electrons for the three cases. The color bar denotes the value of concentrations. 36



Article

Figure 14. Snapshots of electron and Li ion distribution at points A−E in the curve shown in Figure 13.

stresses during charging and discharging. Simulations suggest that the pressure in the central region of the particle is always negative, or compressive, during both charging and discharging. The pressure near the surface changes from negative during charging to positive at the late stages (points D and E) of discharging cycle. The red color denotes the largest positive pressure, that is, tensile stresses, which may cause particle cracking. The largest compressive stresses are observed at the interface of two particles. The largest shear stress increases during charging and decreases during discharging. It is found that the largest shear stress develops at the edges of the interfaces between two particles. Therefore, the interface edges are potential sites for shear-induced cracking. The evolution of Li ion concentration, total charge z1X1 + z2X2 + z3X3, where zi is the valence of species j, and voltage is presented in Figure 17. The main observations are (1) the elastic interaction does not cause inhomogeneities in charged particle distributions; (2) inside the particle and in electrolyte, the total charge is almost zero; (3) near the interface, charge separation takes place leading to the formation of a double layer; and (4) the voltage inside the particles changes with dϕ*/dt* = ±0.009. However, because of isotropic mobilities of charged species leading to a uniform charge distribution inside the nanoparticles beyond the space-charge layer, the voltage distribution inside the nanoparticles is almost uniform. Diffusivity anisotropy (Da/Dc = 1, 0.1, and 0.001) and elastic interaction (lattice mismatch coefficient ε0 = 0.0, 0.05, and 0.08) are used as model parameters to investigate the effect of anisotropic diffusivity and elastic interaction on charge transport. Figure 18 presents the voltage versus capacity curves for a charge rate dϕ*/dt* = ±0.009. It follows from the figure that the elastic interaction slows Li ion insertion at the late stage of charging cycle when Li ion concentration becomes high and stress becomes large. It is expected that the capacity will decrease with increasing lattice mismatch ε0 = 0.0, 0.05, and 0.08. Introducing diffusivity anisotropy leads to a decrease in capacity due to slower Li ion diffusion in the ab-plane.

inhomogeneous along the surface of the particles. Li ions are quickly extracted from the surface layer of the particles during charging, while some residual concentration of Li ions remains at their centers. 5.2.2. Effect of Elastic Interaction on Charged Particle Transport. Li ion insertion/extraction changes the local lattice constant of TiO2 and causes stress generation. The deformation and stresses not only increase the total energy of the system and affect charged particle transport, but also affect the structural stability of the material. For example, cracking is often observed in large particles and at interfaces between particles and carbon.3 In this section, we examine the capability of our model to investigate the effect of elastic interaction on charged particle transport. A voltage versus capacity curve with elastic interaction under a charging rate dϕ*/dt* = 0.009 for Case 1 microstructure is presented in Figure 15. The curve is similar to that in Figure 11; however, the total inserted amount of Li ion is about 0.32, and the remaining Li concentration after discharging is about 0.035. Figure 16 shows the evolution of pressure and shear

6. CONCLUSIONS We have developed an electrochemical phase-field model for a nanocomposite electrode based on the following fundamental assumptions for thermodynamic and kinetic properties: (1) the equilibrium concentrations of Li ions and electrons in

Figure 15. Voltage versus capacity curve during discharging and charging with elastic interaction under a charging rate dϕ*/dt* = 0.009 for the Case 1 microstructure. 37



Article

Figure 16. Snapshots of pressure p/C44 = (σ11 * + σ22 * + σ33 * )/3 and shear stress σ23 * distribution at points A−E in the curve shown in Figure 15.

Figure 17. Snapshots of electron, Li ion, anion A−, and charge distributions at points A−E in the curve shown in Figure 15.

properties and the understanding of physics and mechanisms of different kinetic processes in the system of interest. In general, the input information can be obtained from atomistic simulations and experiments. However, because of the complexity of real systems and the limitations of atomistic simulations, it is challenging to calculate and/or measure some thermodynamic and kinetics properties such as double layer structure, interfacial reaction kinetics, concentration, and temperature dependences of diffusivity. One of the advantages of modeling is that the unknown thermodynamic and kinetic properties can be validated by parametrical simulations and comparison with measurable structure and properties. In addition, it should be pointed out that (1) some important material processes are not taken into account in the current model, for instance, the interfacial reaction kinetics, possible phase transitions, volume changes during Li insertion; and (2)

TiO2nanoparticles depend on voltage and can be obtained from measured voltage versus capacity curves with a slow charging rate; (2) the effective mobility of charged particles depends on voltage, which is related to the measured current versus voltage; and (3) the mobility of species in electrolyte is high so that equilibrium concentrations at the electrolyte side of the interface are always reached. The model takes into account porous microstructure, mobility anisotropy of diffusive species in the solid state, and deformation due to Li ion insertion and extraction. 1D and 2D simulations demonstrate that the model is able to capture the diffusion of charged particles, charging and discharging processes in complex nanocomposites, and predict the effect of microstructure, mobility anisotropy, and elastic interaction on charged particle transport. We would like to emphasize that quantitative phase-field simulations require accurate thermodynamic and kinetic 38



Article

(4) Shim, H. W.; Lee, D. K.; Cho, I. S.; Hong, K. S.; Kim, D. W. Nanotechnology 2010, 21, 255706. (5) Chen, J. S.; Liu, H.; Qiao, S. Z.; Lou, X. W. J. Mater. Chem. 2011, 21, 5687−5692. (6) Goodenough, J. B.; Kim, Y. Chem. Mater. 2010, 22, 587−603. (7) Hu, Y. S.; Kienle, L.; Guo, Y. G.; Maier, J. Adv. Mater. 2006, 18, 1421−1426. (8) Kerisit, S.; Rosso, K. M.; Yang, Z. G.; Liu, J. J. Phys. Chem. C 2009, 113, 20998−21007. (9) Kavan, L.; Gratzel, M.; Gilbert, S. E.; Klemenz, C.; Scheel, H. J. J. Am. Chem. Soc. 1996, 118, 6716−6723. (10) Van de Krol, R.; Goossens, A.; Schoonman, J. J. Phys. Chem. B 1999, 103, 7151−7159. (11) Sudant, G.; Baudrin, E.; Larcher, D.; Tarascon, J. M. J. Mater. Chem. 2005, 15, 1263−1269. (12) Subramanian, V. R.; Boovaragavan, V.; Diwakar, V. D. Electrochem. Solid-State Lett. 2007, 10, A255−A260. (13) Subramanian, V. R.; Diwakar, V. D.; Tapriyal, D. J. Electrochem. Soc. 2005, 152, A2002−A2008. (14) Golmon, S.; Maute, K.; Dunn, M. L. Comput. Struct. 2009, 87, 1567−1579. (15) Johan, M. R.; Arof, A. K. J. Power Sources 2007, 170, 490−494. (16) Subramanian, V. R.; Boovaragavan, V.; Ramadesigan, V.; Arabandi, M. J. Electrochem. Soc. 2009, 156, A260−A271. (17) Cai, L.; White, R. E. J. Power Sources 2011, 196, 5985−5989. (18) Gomadam, P. M.; Weidner, J. W.; Dougal, R. A.; White, R. E. J. Power Sources 2002, 110, 267−284. (19) Santhanagopalan, S.; Guo, Q. Z.; Ramadass, P.; White, R. E. J. Power Sources 2006, 156, 620−628. (20) Chen, L. Q.; Khachaturyan, A. G. Acta Metall. Mater. 1991, 39, 2533−2551. (21) Zhang, M. Y.; Wang, Y. X.; Chen, Z.; Dong, W. P.; Lai, Q. B.; Zhang, L. P. Rare Met. Mater. Eng. 2009, 38, 962−966. (22) Moriguchi, I.; Hidaka, R.; Yamada, H.; Kudo, T.; Murakami, H.; Nakashima, N. Adv. Mater. 2006, 18, 69−73. (23) Li, Y. L.; Hu, S. Y.; Liu, Z. K.; Chen, L. Q. Acta Mater. 2002, 50, 395−411. (24) Koyama, T. Sci. Technol. Adv. Mater. 2008, 9, 013006. (25) Hu, S. Y.; Henager, C. H.; Heinisch, H. L.; Stan, M.; Baskes, M. I.; Valone, S. M. J. Nucl. Mater. 2009, 392, 292−300. (26) Tonks, M.; Gaston, D.; Perrnann, C.; Millett, P.; Hansen, G.; Wolf, D. Nucl. Eng. Des. 2010, 240, 2877−2883. (27) Garcia, R. E.; Bishop, C. M.; Carter, W. C. Acta Mater. 2004, 52, 11−21. (28) Han, B. C.; Van der Ven, A.; Morgan, D.; Ceder, G. Electrochim. Acta 2004, 49, 4691−4699. (29) Guyer, J. E.; Boettinger, W. J.; Warren, J. A.; McFadden, G. B. Phys. Rev. E 2004, 69, 21603. (30) Guyer, J. E.; Boettinger, W. J.; Warren, J. A.; McFadden, G. B. Phys. Rev. E 2004, 69, 21604. (31) Okajima, Y.; Shibuta, Y.; Suzuki, T. Comput. Mater. Sci. 2010, 50, 118−124. (32) Gathright, W.; Jensen, M.; Lewis, D. J. Mater. Sci. 2012, 47, 1677−1683. (33) Bai, P.; Cogswell, D. A.; Bazant, M. Z. Nano Lett. 2011, 11, 4890−4896. (34) Tang, M.; Huang, H. Y.; Meethong, N.; Kao, Y. H.; Carter, W. C.; Chiang, Y. M. Chem. Mater. 2009, 21, 1557−1571. (35) Singh, G. K.; Ceder, G.; Bazant, M. Z. Electrochim. Acta 2008, 53, 7599−7613. (36) Bazant, M. Z., http://arxiv.org/abs/1208.15872012. (37) Cogswell, D. A.; Bazant, M. Z. ACS Nano 2012, 6, 2215−2225. (38) Ferguson, T. R.; Bazant, M. Z. J. Electrochem. Soc. 2012, 1204− 2934; http://arxiv.org/abs/. (39) Sushko, M. L.; Rosso, K. M.; Liu, J. J. Phys. Chem. Lett. 2010, 1, 1967−1972. (40) Sushko, M. L.; Rosso, K. M.; Liu, J. J. Phys. Chem. C 2010, 114, 20277−20283.

Figure 18. Effect of diffusivity anisotropy and elastic interaction on the charging process.

the phase-field model only investigates microstructure and kinetic property evolution in a very small representative volume to capture the interfacial properties in nanocomposites. Predicting the effect of microstructure of nanocomposites on the whole battery performance requires a multiscale simulation tool, which enables one to deliver the thermodynamic and kinetic knowledge from smaller length scale to larger length scale. In the present work, the simulations demonstrate the predictive capability of phase-field model for given thermodynamic and kinetics properties, but more physics has to be integrated into the model for real battery simulation applications. Despite all of these limitations of the current model, simulated quantities such as equilibrium concentrations, voltage versus capacity, and current versus voltage curves are in qualitative and quantitative agreement with experimental observations,1,42 suggesting that the model can be used to screen materials and electrode design options, and help predict the optimum properties for enhanced performance.

■

AUTHOR INFORMATION

Corresponding Author

*Tel.: (509) 371-6928. Fax: (509) 375-3033. E-mail: shenyang. [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work is supported by the Laboratory-Directed Research and Development Program at Pacific Northwest National Laboratory (PNNL). PNNL is a multiprogram national laboratory operated for DOE by Battelle under contract DEAC05-76RL01830.

■

REFERENCES

(1) Yang, Z. G.; Choi, D.; Kerisit, S.; Rosso, K. M.; Wang, D. H.; Zhang, J.; Graff, G.; Liu, J. J. Power Sources 2009, 192, 588−598. (2) Hu, Y. S.; Demir-Cakan, R.; Titirici, M. M.; Muller, J. O.; Schlogl, R.; Antonietti, M.; Maier, J. Angew. Chem., Int. Ed. 2008, 47, 1645− 1649. (3) Bridel, J. S.; Azais, T.; Morcrette, M.; Tarascon, J. M.; Larcher, D. Chem. Mater. 2010, 22, 1229−1241. 39



Article

(41) Sushko, M. L.; Rosso, K. M.; Zhang, J. G.; Liu, J. J. Phys. Chem. Lett. 2011, 2, 2352−2356. (42) Saravanan, K.; Ananthanarayanan, K.; Balaya, P. Energy Environ. Sci. 2010, 3, 939−948. (43) Koci, L.; Kim, D. Y.; de Almeida, J. S.; Mattesini, M.; Isaev, E.; Ahuja, R. J. Phys.: Condens. Matter 2008, 20, 345218. (44) Swamy, V.; Gale, J. D.; Dubrovinsky, L. S. J. Phys. Chem. Solids 2001, 62, 887−895. (45) Hu, S. Y.; Chen, L. Q. Acta Mater. 2001, 49, 1879−1890. (46) Parker, R. A. Phys. Rev. 1961, 124, 1719−1722. (47) Chen, L. Q.; Shen, J. Comput. Phys. Commun. 1998, 108, 147− 158.

40


Mesoscale Phase-Field Modeling of Charge Transport in

Recommend Documents