Multi-physics simulations of lithiation-induced stress in Li{}_{\rm 1+x}Ti{}_{2}O{}_{4} electrode particles

Multi-physics simulations of lithiation-induced stress in LiTiO electrode particles

Tonghu Jiang Department of Materials Science and Engineering, Johns Hopkins University, Baltimore, MD 21218 USA    Shiva Rudraraju Department of Mechanical Engineering, University of Michigan, Ann Arbor, MI 48109 USA    Anindya Roy Department of Materials Science and Engineering, Johns Hopkins University, Baltimore, MD 21218 USA    Anton Van der Ven Materials Department, University of California, Santa Barbara, CA 93106 USA    Krishna Garikipati Department of Mechanical Engineering, University of Michigan, Ann Arbor, MI 48109 USA Department of Mathematics, University of Michigan, Ann Arbor, MI 48109 USA    Michael L. Falk Department of Materials Science and Engineering, Johns Hopkins University, Baltimore, MD 21218 USA Department of Mechanical Engineering and Department of Physics and Astronomy, Johns Hopkins University, Baltimore, MD 21218 USA

Cubic spinel LiTiO is a promising electrode material as it exhibits a high lithium diffusivity and undergoes minimal changes in lattice parameters during lithiation and delithiation, thereby ensuring favorable cycleability. The present work is a multi-physics and multi-scale study of LiTiO that combines first principles computations of thermodynamic and kinetic properties with continuum scale modeling of lithiation-delithiation kinetics. Density functional theory calculations and statistical mechanics methods are used to calculate lattice parameters, elastic coefficients, thermodynamic potentials, migration barriers and Li diffusion coefficients. These quantities then inform a phase field framework to model the coupled chemo-mechanical evolution of electrode particles. Several case studies accounting for either homogeneous or heterogeneous nucleation are considered to explore the temporal evolution of maximum principle stress values, which serve to indicate stress localization and the potential for crack initiation, during lithiation and delithiation.

I Introduction

Scientists and engineers are exploring a vast materials space for cheaper, higher energy-density, and reliable battery materials, in response to the growing demand of portable electronics and electric vehicles.Hautier et al. (2011) Understanding the suitability of the cathode and anode materials, both computationally and experimentally, has been central to battery research. Li-ion batteries are the most used rechargeable batteries today, and a considerable fraction of battery research is focused on improving the characteristics of Li-ion batteries and finding better candidates for the electrodes.

In this work, we are interested in LiTiO for its possible use as an electrode material. LiTiO is known to form different polymorphs upon lithium insertion, and the common phases include rutile, anatase, spinel, and brookite.Cava et al. (1984); Zachau-Christiansen et al. (1988) The electric potential for lithium insertion into these structures is relatively low (1.5V ), enhancing the suitability of LiTiO as an anode material. Cubic spinel LiTiO and LiTiO have been shown in experiments to exhibit high lithium diffusivity, and a negligible change of lattice parameter during (dis)charging.Krtil and Fattakhova (2001); Colbow et al. (1989); Wang et al. (1999) A large change in stress during the charging process can result, via the elastic coefficients, from significant changes in lattice parameters. The high stress can initiate fracture, limiting cycleability. Low variations in lattice parameters can thus promote better cycleability.

This work presents a multi-physics and multi-scale case-study in which we investigated LiTiO computationally, with the aim of predicting, from first principles, the stresses that will arise within the material during repeated charging and discharging of electrode particles. In particular, we are interested in spinel LiTiO, which belongs to the space group . In spinel LiTiO, 16 and 32 sites are occupied by Ti and O, respectively. We emphasize that the goal of the present work is to demonstrate how the various computational techniques at different length scales could be applied in tandem, and we used a high temperature of 800 K to magnify the kinetic effects. While specific quantitative predictions regarding electrode failure are not possible without continued development of these methods, the work presented here demonstrates that interesting insights into the chemo-mechanical processes that can limit cycleability can nonetheless be obtained.

During charge and discharge, LiTiO undergoes first order phase transformations. For x0, LiTiO is a solid solution phase, in which the tetrahedral 8 sites are occupied by lithium ions; after lithium insertion beyond x = 0, a second phase forms, in which lithium ions reside in octahedral 16 sites. Between x = 0 and 1, lithium intercalation in LiTiO involves a two-phase coexistence process. Figure 1 shows the structures of the ideal and phases in (a) and (b). The thermodynamic and electronic properties of LiTiO have been studied by experiments and first-principles calculations.Krtil and Fattakhova (2001); Wagemaker et al. (2005); Bhattacharya and Van der Ven (2010)

Figure 1: (a) and (b): The ideal structure of phase LiTiO and phase LiTiO, respectively. The lithium ions reside in the tetrahedral sites and the octahedral sites formed by oxygen anions in respective LiTiO and LiTiO. (c) Connection between lithium sites. The tetrahedral sites are represented by larger atoms than octahedral sites.

Lithium diffusion mechanisms and phase boundary mobility strongly influence overall battery performance. Important factors that affect macroscopic lithium insertion and removal during a topotactic two-phase reaction within a crystalline electrode are interfacial free energies and coherency strain energy. Lithium ion diffusion within the interface region must sustain the propagation of the phase boundary in a direction that will lower the overall free energy. Lithium ion diffusion in most solid solution phases of electrode materials is anisotropic due to the crystal structure. However, because of its cubic symmetry, lithium diffuses isotropically in the solid solution phase of the spinel LiTiO . This is in contrast to two-dimensional diffusion in layered LiCoOVan der Ven et al. (2001) and one-dimensional diffusion in LiFePO.Ouyang et al. (2004); Morgan et al. (2004) As for the coherency strain energy of two-phase coexistence, the change in the lattice parameter is observed to be less than 1% between the and phases of LiTiO, and hence may not induce significant deformation of the interface. It is therefore a reasonable approximation to calculate homogeneous free energies at equilibrium lattice parameters from first principles and combine these with the ab initio computed change in lattice parameters and elastic moduli in continuum scale simulations to study the coupling between mechanics and chemistry during intercalation and deintercalation.

We have used an array of computational tools to study the thermodynamics and kinetics of LiTiO. Density functional theory computations were used to obtain formation energies of LiTiO with different lithium configurations and migration barriers for lithium ion hops. The cluster expansion method was employed to parameterize the formation energy and migration barrier dependence on Li-vacancy disorder within LiTiO. Cluster expansions and the Metropolis Monte Carlo (MMC) method were used to determine thermodynamic potentials, such as the composition-dependent Gibbs free energy and interfacial free energies. Mobility and diffusion in the solid solution phase were estimated using a combination of first-principles and kinetic Monte Carlo (KMC) computations.

Informed by quantum-mechanics- and statistical-mechanics computations, we then performed continuum scale computations for Li-ion diffusion and the mechanical stress induced by phase transformations in a phase field framework. The continuum scale computations used the following inputs from first principles computations: the free energies of the two stable phases, phase kinetics, interfacial energy and changes in the lattice parameters and elastic moduli. These computations are helpful in understanding the charge-discharge (lithiation-delithiation) kinetics and mechanical deformation that can induce mechanical failure in the anode particles. We studied the effects of different charge-discharge cycles on phase boundary movement, charge localization and peak stress evolution.

In the next section we elaborate the computational methods used in this work. Sec. III contains the main results of this work – the multi-physics approach applied to the study Li-ion kinetics in LiTiO.

Ii Methods

In this section we describe the computational methods in this multi-physics approach, from atomistic first principles calculations to the continuum scale formulation. The thermodynamics and kinetics of the and phases of LiTiO have been previously studied by Bhattacharya et al.Bhattacharya and Van der Ven (2010) In this work, we use first principles computations of the thermodynamics, kinetics and elastic constants to inform continuum scale computations. This includes the composition dependence of lattice parameters and elastic coefficients. The continuum scale computations use a coupled phase field and finite strain mechanics multi-physics formulation. The use of the finite strain (nonlinear elasticity) formulation distinguishes this work, as phase field computations are traditionally limited to the infinitesimal strain approximation and thus only use linear elasticity formulations. In Sections II.1II.3 we discuss the first-principles methods, followed by the thermodynamics and the Li-ion kinetics calculated for the LiTiO system. In Section II.4, we elaborate on the continuum scale formulation.

ii.1 First principles calculations

We constructed a first-principles cluster expansion to describe the energy of LiTiO as a function of Li-vacancy order/disorder. To accomplish this, we borrowed formation energies of 48 Li-vacancy orderings in small supercells of the primitive LiTiO unit cell and 4 migration barriers from the work of Bhattacharya et al.,Bhattacharya and Van der Ven (2010) and extended this dataset by adding more energy information sampled in a larger configuration space. We used the VASPKresse and Furthmüller (1996a, b) plane wave code to relax structures with different lithium configurations and calculate their total energies. The parameters used in the calculations are referred to in Ref. 8.

The lattice parameters of LiTiO can be extracted from the relaxed structures. Because the cell shape does not necessarily remain cubic after relaxation, we calculated the volumes of the relaxed structures and then converted them to lattice parameters assuming the supercells were cubic. The elastic constants are key parameters for studying the intercalation-induced mechanical response of LiTiO. Because off-stoichiometric Li concentrations in LiTiOusually have arrangements that break the cubic symmetry of ideal LiTiO or LiTiO, it becomes impractical to deform the structures and extract elastic constants from the energy-strain relationship. We calculated elastic constants from the stress-strain relationship using finite-differences, obtained using VASP.Le Page and Saxe (2002) To calculate the composition dependence of elastic constants, we followed the method by Liu et al.Liu et al. (2005) to symmetrize the elastic constants by transforming the original elastic constants by the 48 symmetry operations belonging to the cubic point group.

The migration barriers of lithium were calculated with the nudged elastic band method implemented in VASP. We identified 18 symmetrically distinct Li hops between pairs of stable tetrahedral and octahedral sites over a range of Li concentrations in addition to the four lithium hops considered by Bhattacharya et al. Bhattacharya and Van der Ven (2010). The migration barriers were calculated in a conventional spinel cell.

ii.2 Cluster expansion

The cluster expansion method by Sanchez et al.Sanchez et al. (1984); Laks et al. (1992); Fontaine (1994) is a mathematical tool to describe the configuration dependence of any property of a multi-component crystal and has contributed greatly to the development of modern alloy theory in oxides.Ceder et al. (2000) It has been successfully applied in the study of various crystalline materials, particularly when analyzing phase stability, order-disorder transitions and diffusion.Wolverton and Zunger (1998); Van der Ven et al. (2001) In this study, we are interested in cluster expanding formation energies and migration barriers. The 0 K formation energy of a particular Li-vacancy ordering (labeled ) within LiTiO is defined as:


where is the total energy of a Li-vacancy ordering and where and are total energies of fully lithiated spinel (all octahedral sites occupied by Li) and delithiated spinel, respectively. The occupation of each lithium site, , within LiTiOis represented by a site occupation variable that is +1 if occupied by a lithium atom and -1 if vacant. The Li-vacancy ordering within LiTiO can then be uniquely determined by specifying all the occupation variables . A cluster expansion expresses the formation energy as a sum of products of effective interactions and polynomials of occupation variables associated with clusters of sites according to Sanchez et al. (1984):


The expansion coefficients (, , etc.) are refered to as effective cluster interactions (ECI) and can be determined from first-principles using a variety of inversion methods that minimize a cross validation score van de Walle and Ceder (2002). Models (i.e. selection of non-zero ECI in a truncated cluster expansion) can be generated with genetic algorithms Hart et al. (2005), by deploying Bayesian approaches using prior knowledge,Mueller and Ceder (2009) or by implementing compressive sensing approaches.Nelson et al. (2013) We applied standard methodologies in our work, and the reader is referred to the cited references for further details.

ii.3 Thermodynamics and Kinetics

First-principles calculations and cluster expansions together provide us with the composition dependence of lattice parameters, formation energies, and elastic coefficients at zero-temperature. The application of Monte Carlo simulations to the cluster expansion of the formation energy enables the calculation of the free energies of bulk phases and interfaces at finite temperature, as described below. Li-ion kinetics is also determined at finite temperature using kinetic Monte Carlo methods that rely on cluster expansions to describe the configuration dependence of migration barriers and the energies of the end states of Li hops.

ii.3.1 Homogeneous free energy

We employed thermodynamic integration, a commonly used free energy calculation technique, to calculate the grand thermodynamic potential ,van de Walle and Ceder (2002); Binder et al. (2011); Frenkel et al. (1997) which is related to the Gibbs free energy according to . At constant chemical potential the grand thermodynamic potential can be calculated using:


At constant temperature, the grand thermodynamic potential can be calculated with:


ii.3.2 Interfacial free energy

Interfacial energy and surface energy play important roles when the size of the electrode crystallite approaches the nanoscale. Following Binder,Binder et al. (2011) we calculated interfacial free energies, , by thermodynamic integration from to of the expression:


The is the excess energy due to the existence of an interface, defined by the following equation:


Here is the total energy of the system of two coexisting phases, while and are the total energies of pure and phases of the same size as in the studied system, and is the fraction of phase in the system.

ii.3.3 Lithium ion mobility calculation

Mobility in a homogeneous crystalline solid is a quantity that describes collective transport of mobile atoms or ions. We first calculated the self-diffusion coefficient using kinetic Monte Carlo simulations to approximate the Kubo-Green expressionRichards (1977)


In the above equation, is the displacement of the -th Li ion after time , is the number of Li ions, and is the dimension of the crystal. The self-diffusion coefficient measures the random walk process of the geometric center of mass of the mobile lithium ions. The lithium ion mobility can be obtained by converting through:


where is the composition of the homogeneous system.

ii.4 The phase field model

Phase field models are commonly employed to study the evolution of spatially continuous order parameters by gradient flow kinetics. They have been widely adopted for modeling phase transformations described by the evolution of conserved (Cahn-Hilliard model Cahn and Hilliard (1958)) or non-conserved order parameters (Allen-Cahn modelAllen and Cahn (1972)). In this work, the existence of stable and phases at Lithium compositions of  (LiTiO) and  (LiTiO), respectively, and a two-phase co-existence at intermediate values of the composition allows for a classical double-well representation of the free energy as illustrated in Figure 3.

In modeling the diffusive kinetics and mechanical response of electrode particles, phase field treatments have traditionally taken into account the interfacial energy, elastic strain energy and anisotropies of the crystalline electrode material. Past studies used this approach to study the formation and growth of the solid-electrolyte interface layer,Deng et al. (2013) revealing a diffusion-limited process. Similarly, Kao et al. Kao et al. (2010) found via experiment and phase field modeling that the phase transition pathway depends on the overpotential in the lithium ion phosphate battery. Further, Lithium intercalation into nanometer-sized electrode particles,Burch and Bazant (2009) charge transport in TiO,Hu et al. (2013) and large mechanical stress due to phase segregation in LiMnOHuttin and Kamlah (2012) have also been modeled using the phase field method.

In this work, we consider a phase-field model driven by anisotropic interfacial and elastic strain energies obtained from first principles computations. Traditionally, most phase field treatments coupled with mechanics were limited to the infinitesimal-strain assumption, which does not satisfy frame invariance of the elastic strain energy function. In the presence of large strains, this induces spuriously high stresses due to unaccounted rotations and renders the model susceptible to predicting failure by the wrong mechanisms. However, our treatment is framed within nonlinear kinematics (finite strains), which nullifies the effect of rotations on the elastic strain energy with mathematical exactness. As stated in the introduction, for spinel , the diffusion of lithium is close to isotropic due to high symmetry of the lattice. We therefore consider an isotropic mobility for lithium-ion diffusion in the stable phases and the interface.

In this model, we consider the following free energy:


where is the composition, is the homogeneous free energy density and is the Green-Lagrange strain tensor. The latter is a frame invariant measure of strain, whose components are given by . Also in component form is the deformation gradient tensor and is the displacement vector.

The first two terms in the integrand of Equation (9) represent the homogenous free energy density and the orientation dependent gradient energy density of the interface, respectively. The form of the homogeneous free energy considered is shown in Figure 3 and the orientation dependence of the interfacial energy density is given by Figure 4. Here, is related to the interfacial energy asCahn and Hilliard (1958)


where is the difference between homogeneous free energy density and the common tangent line between compositions at and (Figure 3). The exact form of homogeneous free energy density, , is not very influential in the resulting phase microstructure or the dynamics by which it evolves.Vaithyanathan et al. (2002) The magnitude of is usually adjusted in order to reproduce the experimental or calculated interfacial free energy so long as the resulting interface width in the continuum scale computation is reasonable.

The third term in the integrand of Equation (9) is the elastic strain energy density, which is a function of the Green-Lagrange strain and the composition. Assuming a St. Venant-Kirchhoff type material model, the specific form of the elastic strain energy density considered here is given by:


where is the fourth order elasticity tensor with major and minor symmetries and cubic anisotropy, and is the stress-free strain introduced to model changes in the lattice parameter as a function of the composition. The composition dependence of the elastic moduli and of the lattice parameter are obtained from first-principles computations described in the next section, and the composition dependence is depicted in Figure 2.

Using a variational approach, we can derive the governing equations for the dynamics of non-equilibrium chemistry and for mechanical equilibrium. The local chemical potential is obtained by evaluating the variational derivative of the free energy , and is given by


The conservation of mass governs the non-equilibrium chemistry, and is given by


where is the constant mobility. This partial differential equation is complemented by boundary conditions:


where the first of the above boundary conditions is a consequence of assuming equilibrium at boundaries, while the second represents a boundary influx, . Taking variations with respect to the displacement, we obtain the governing equation of quasi-static mechanical equilibrium (conservation of linear momentum), on the basis that elastic equilibrium is established much more rapidly than chemical equilibrium:

Following standard variational arguments this condition leads, via Euler-Lagrange equations, to the following partial differential equation:


where, are components of the stress tensor. The partial differential equation in (15) is subject to boundary conditions:


where and are the Dirichlet and Neumann boundaries, respectively. In this formulation, the stress tensor inherits a composition dependence from the strain energy density .

Equations (13) and (15) are the governing partial differential equations for non-equilibrium chemistry and for mechanical equilibrium, respectively. We use the Finite Element Method (FEM) to solve these partial differential equations in a weak (integral) formulation. The details of the weak formulation, finite dimensional discretization and the solution schemes adopted are beyond the scope of this paper, but the methods considered are standard in the FEM literature. The numerical framework to solve the system of equations was implemented in an in-house C++ FEM code built on top of the deal.II library,Bangerth et al. (2007) and uses the Sacado library of the Trilinos project Heroux et al. (2005) for algorithmic differentiation to generate the Jacobian matrix. Some of the details of the weak formulation, finite dimensional discretization and numerical framework can be found in related publications by the authors.Rudraraju et al. (2014, 2016)

Iii Results and Discussion

A significant fraction of the atomistic calculations needed to connect to the continuum scale calculations have been reported before by Bhattacharya et al. in Ref. 8, and we restrict our report to the results from the atomistic calculations that are new, or essential.

iii.1 Results of first-principles calculations

iii.1.1 Lattice parameters and elastic constants

Figure 2 (top panel) shows lattice parameters as a function of number of lithium ions in a conventional spinel cell. In the common phase field models the lattice parameter is assumed to be linearly dependent on composition following Vegard’s law. In Figure 2 we see that the lattice parameter contracts after the phase transition from to , while in the phase the lattice parameter increases with lithium ion density. Although the trend of lattice parameter is clear, the relative difference between LiTiO and LiTiO is small, i.e, 1% for this phase transition.

Figure 2: (Top panel) lattice parameter of LiTiO calculated by DFT. Each red dot represents a lattice parameter of a unique structure. The dots in circle (in blue) represents the lattice parameter of the structure with the lowest energy. (Bottom panel) Red squares, green triangles and purple circles represent elastic coefficients C, C, and C respectively. The solid lines through the group of points are third-order-polynomial fit to the calculated elastic coefficients.

Figure 2 (bottom panel) shows the lithium composition dependence of the elastic coefficients C, C and C. The individual points correspond to calculated values, whereas the lines drawn through them are third order polynomial fits. As can be seen in Figure 2, C increases monotonically as the lithium composition increases. C increases rapidly in the solution phase from x1 to x0 and slowly in the two phase region. The overall increase of C from x=1 to x=1 could be as large as 50 GPa or 30% of its value at x1. C, C and C generally depend on lithium composition in a nonmonotonic way. The assumption made in many continuum scale modelsHu and Chen (2001) that elastic constants depend linearly on composition may not be satisfied particularly for C and C.

iii.2 Cluster expansions

iii.2.1 Cluster expansion of LiTiO 

We are interested in building a cluster expansion capable of predicting energies in single phases and within interfaces. The lithium configurations at the interface region correspond to high energy states; thus they are not well described by cluster expansions that have been optimized to more accurately predict low energy configurations. Seko et al. proposed a procedure that improves the cluster expansion model by adding structure samples in a systematic way such that the correlation between the structures is reduced.Seko et al. (2009) In this work, we simply add more training configurations to the cluster expansion fit that have random lithium arrangements in a conventional cell. Thirty-six structures corresponding to the two end states of 18 lithium transitions in the NEB calculations were added in the cluster expansion fit. This improved the cluster expansion prediction of energy differences between the end states of Li hops, and therefore also led to a more accurate prediction of migration barriers, as discussed in Section II.2. No more than 50 clusters were included in the cluster expansion fit to avoid overfitting, and a total of 92 formation energies of different configurations were used in the cluster expansion fit. The final cluster expansion includes 11 pair clusters, 9 triple clusters, and 3 quadruple clusters. The error and cross validation score of the cluster expansion is 18 meV and 27 meV per primitive cell.

To verify that this cluster expansion yields the correct phase behavior at 0 K, a series of Metropolis Monte Carlo calculations were performed that started at high temperature and that were slowly cooled down to 0 K to search for the lowest energy states. The lowest energy states obtained from Metropolis Monte Carlo calculation were then confirmed with a genetic algorithm search. The results obtained from these calculations agree with the previous calculations reported by Bhattacharya in Ref. 8, and we refer the reader to Figure 2 therein and the associated text.

iii.3 Thermodynamics and Kinetics

iii.3.1 Homogeneous free energy

Figure 3 shows the calculated Gibbs free energies of the and phases with Monte Carlo data obtained in a 444 supercell at T=800 K, denoted by blue dots. We can see that the phase is stable in only a very narrow lithium composition interval around 1. The homogeneous free energy was fit with a double-well function (the solid black curve in Figure 3) for use in the phase field model.

Figure 3: Homogenous free energy representation obtained by fitting a double-well polynomial to the free energy values obtained from grand canonical Monte Carlo calculations. The blue dots are the free energy values obtained from grand canonical Monte Carlo calculations.

iii.3.2 Interfacial free energy

We constructed 844 and 1055 Monte Carlo supercells to calculate finite temperature interfacial free energies (111 is a conventional cell). The periodic screw boundary condition was imposed on the system in the Monte Carlo computation. The implementation details of this boundary condition can be found inJiang (2013). A total of 15 nonequivalent orientations of the formula , where M, N = 0,1,2,4, were selected in the 844 supercell, and similarly 21 orientations were selected for the 1055 supercell. The Monte Carlo computations started out at 1200 K (below the order-disorder temperature) with sharp interfaces and were run for 15000 sweeps. The temperature was lowered by 50 K after every other sweep until the temperature reached 0 K. The final systems at 0 K are lower in energy than the initial systems with sharp interfaces. However, they are not guaranteed to be the lowest energy states because the final energy may be cooling-rate dependent.

Since the continuum scale phase field computations are performed on a 2D domain, the interfacial energy anisotropy within the plane (001) is considered, and the long axis is directed along [110]. The resulting polar plot of the interfacial energy at 800 K as calculated in the 1055 system is shown in Figure 4.

Figure 4: Polar plot of 2D interfacial free energy (eV/Å). The 0 deg and 45 deg correspond to [100] and [110] direction, respectively.

iii.3.3 Lithium ion kinetics

High lithium ion mobility, which in part determines the battery charge/discharge rate, is one of the major criteria for selecting promising lithium ion battery electrodes.Sebastian and Gopalakrishnan (2003); Van der Ven et al. (2013) As most intercalation electrodes experience phase transformations during battery charge/discharge (lithiation/delithiation), lithium ion diffusion within the crystallites is determined not only by its mobility in the single phases but also across interfaces in the two-phase region. Despite its importance, only an “apparent diffusion coefficient” is usually measured for the two-phase region by fitting experimental data to equations derived from Fick’s law of diffusion and its reliability is doubtful.Zhu and Wang (2010) Numerical results with phase field models show that the error induced is large when the lithium composition reaches the two phase region, especially within the spinodal region.Han et al. (2004) Attempts have been made to develop models that explicitly treat phase boundary movement, e.g.,the moving boundary model,Funabiki et al. (1999) and the mixed control phase transformation model.Zhu and Wang (2010)

Lithium ion transport in the electrodes may be either limited by diffusion in the stable phases or controlled by the mobility of the interface. It has been reported that the rate of lithium insertion into graphite during two-phase reactions is determined by diffusion in the stable phases.Funabiki et al. (1999) In LiTiO, in contrast, the measured diffusion coefficient in the stable phases is much higher than that in the two-phase region, which suggests that phase boundary mobility is rate limiting during Li insertion or removal.Li et al. (2012) Analysis of lithium diffusion in vanadium oxide films also suggests that the kinetics of the two phase process is controlled by the phase boundary mobility.Bae and Pyun (1996) For more complex cases, asymmetry of lithium ion mobility in stable phases leads to asymmetrical charge/discharge processLi et al. (2012) and pinning of phase boundary movement.Choi et al. (1998) Despite this dependency, most continuum scale models assume a diffusion-limited process to simplify the problem.Hong et al. (2006)

In the present work we systematically studied the mobility of the interface, and found that the mobility in the interface region is about three orders of magnitude higher than the phase, indicating that the mobility is in fact diffusion limited.Jiang (2013) For the continuum scale computations in the solid-solution phase, we consider a constant mobility corresponding to the self diffusion coefficient value of 10 cm/s, and as explained earlier, the mobility is assumed to be orientation independent.

We calculated the lithium ion diffusion coefficient in the solution phase ranging from 1 to  0, and related it to the global mobility via Eq. 8. The kinetic Monte Carlo simulations of diffusion in the phase using an energy landscape searching algorithm indicated that hops occurred exclusively between nearest neighbor tetrahedral sites; therefore only this type of transition was allowed in subsequent kinetic Monte Carlo simulations and the energy landscape searching algorithm was no longer used. Because the lithium diffusion coefficient in the solution phase of LiTiO is isotropic, a simple periodic boundary condition for studying diffusion of lithium in (100) was used, with a 666 unit cell at T=800 K. Straightforward KMC calculations with 500 to 1000 sweeps (in each KMC calculation, at each Li composition) were performed, and the corresponding diffusion coefficients were then calculated with the help of Eq. 7. Figure 5 shows the calculated mobility in the solid solution phase at T=800 K. The trend of the curves is the same as the diffusion coefficient calculated by Bhattacharya et al.Bhattacharya and Van der Ven (2010) The diffusion coefficient reported there was smaller as it was calculated at room temperature, in contrast to 800 K in the present work.

Figure 5: Self-diffusion coefficient as a function of composition in the phase.

iii.4 Continuum scale simulations

The continuum scale computations are primarily aimed at (a) understanding the evolution of lithium ion composition and diffusion kinetics during charging/discharging cycles, and studying the effects of heterogenous nucleation sites or defects on the overall diffusion kinetics, and (b) studying the mechanical response and stress variations during the diffusion processes due to the dependence of the elastic moduli and lattice parameters on the lithium ion composition. These two studies have been constructed so as to provide a qualitative picture of the efficiency of lithiation/delithiation processes in terms of charge distribution and localization, and also to understand the effect of spatial and temporal stress fluctuations on the material degradation that can affect the cycleability of these electrode particles. First, we describe the problem setup (geometry, initial conditions and boundary conditions), and then present the results of the continuum scale studies.

iii.4.1 Problem setup

We consider a circular domain of radius as the geometry of the electrode particle, and all the computations presented here are restricted to two dimensions (Figure 6). Mechanically, the electrode particles are only subjected to a minimal set of Dirichlet boundary conditions on the displacement field to prevent rigid body motion. The composition field is driven by radially symmetric flux (Neumann) boundary conditions that induce an inward/outward flux of lithium ions. The temporal variation and cycling of the boundary flux (inward/outward) are depicted in Figure 7. These conditions on the displacement field and the Lithium ion flux remove any null space from this two-field (composition and displacement) multi-physics formulation and ensure the numerical stability of the simulations. Finally, since the computations are initial boundary value problems (IBVPs), one also needs to specify initial conditions on the composition field. We use the initial conditions to specify either a system in which homogenous nucleation dominates (no heterogeneous nucleation sites) or a system in which nucleation is dominated by the presence of nucleation sites, modeled here as circular regions of radius (Figure 6). The motivation for including nucleation sites is driven by the physical understanding that these electrode particles are not spatially uniform and may have material defects or other heterogeneities, which enhance nucleation. It is expected that heterogeneous nucleation dominates over homogeneous nucleation in all except defect-free systems.

Recall that since the diffusion kinetics is very slow compared to elastic wave speeds, we assume the material to always be in elastic equilibrium and neglect inertial effects. Thus we model the mechanics problem as quasi-static without inertial terms. The initial conditions for the mechanics problem are zero displacement field with no pre-stress.

Figure 6: Problem geometry (circular domain of radius ) and composition initial conditions corresponding to homogenous nucleation (no nucleation sites) and heterogenous nucleation (2, 3 and 5 nucleation sites). Shown are the locations of the nucleation sites, each modeled as a circular region of radius with a composition value . The legend indicates the Lithium ion composition.
Figure 7: The temporal variation of the Lithium ion flux applied at the boundary to model charging (lithiation) and discharging (delithiation) cycles. The flux magnitude considered here models an initial charging (), followed by a fast discharging () and charging () and finally slow discharging ().

iii.4.2 Lithium diffusion kinetics

In this study, we consider four different IBVPs, each with different initial conditions used to model the homogenous/heterogenous nucleation. The first set of computations were carried out assuming a homogenous domain with no nucleation sites, and the subsequent computations considered multiple spatially distributed nucleation sites. Here, we only present results for the problems with two, three, and five nucleation sites whose sizes and spatial distributions are shown in Figure 6. As previously stated, the lithiation/delithiation process is modelled by the application of an axisymmetric inward/outward flux on the surface of the circular domain. The time-varying flux profile considered in these studies is shown in Figure 7, where a positive value indicates influx and a negative value indicates outflux.

The temporal evolution of lithium composition during the lithiation/delithiation cycles shows patterns of phase evolution within the electrode particle. Figure 8 shows the time evolution of lithium ion composition in the presence of zero, two and five nucleation sites. With no nucleation sites there is no breaking of the axisymmetry of the composition field. However, even in this case, the rapid lithiation/delithiation leaves behind regions of depletion () and accumulation (). This suggests that to avoid such phase localization within electrode particles, the lithiation/delithiation rates must be carefully controlled. These patterns of lithium localization are more prominent as well as complex in computations with heterogeneous nucleation sites, which break the symmetry. In these cases, two important observations are that (i) lithium depletion/accumulation occurs primarily around nucleation sites, and that (ii) heterogeneous nucleation leads to complex spatial patterns of lithium distribution as seen in the two nuclei and five nuclei case at time , , and .

Figure 8: Time snapshots of Lithium ion composition during the lithiation/delithiation cycles. Shown are the composition fields at time , , and for the problems with no nuclei, two nuclei and five nuclei. The blue regions with composition close to are in the phase and the red regions with composition close to are in the phase. Intermediate values of the composition correspond to the two-phase regions.

iii.4.3 Mechanical response and stress evolution

An important factor in the performance degradation of electrode materials is the mechanical failure of electrode particles through fractureKlinsmann et al. (2016); Miehe et al. (2016). While we do not explicitly model material degradation by fracture/damage mechanisms in the current work, we can still obtain a useful qualitative picture of material degradation by tracking the stress distribution and evolution of peak stress values during the lithiation/delithiation cycles. The peak stress values are important, since fracture and voiding in brittle and ductile materials, respectively, are initiated at points with critical values of the relevant stress measures such as the maximum principal stress and the related Mode-I/Mode-II opening tractions for fracture, and the maximum hydrostatic stress for voiding. In this work we identify the temporal evolution of the maximum principle stress values from which the other measures can be obtained. The spatial distribution of stress and stress localization observed in the presence of zero, two and five nucleation sites are shown in Figure 9.

The computations indicate that the peak stress values occur at the surfaces of the nucleation sites and at the sites of lithium localization. This spatial distribution of peak stresses is explained by the fact that the interfaces are phase boundaries characterized by discontinuously changing lattice parameters and elastic moduli, which induce large stresses. Further, the spatial distribution of peak stress values is strongly dependent on the heterogeneity induced by the nucleation sites. As seen in Figure 9, the peak stress profiles are significantly different for the IBVPs with two and five heterogeneous nucleation sites when compared to the case with no heterogeneous nucleation sites.

Using a maximum principle stress criterion for fracture, one can expect to see crack initiation and evolution from material points in the domain where the peak maximum principle stress value exceeds a critical value. In the absence of data about the critical value of the maximum principle stress for crack initiation, we examine the evolution of peak stress values in the electrode particle and the effect of lithiation/delithiation cycles by plotting the time evolution of the peak maximum principle stress magnitude (Figure 10). These simulations indicate that the peak stress values are higher in the case of heterogenous nucleation and that the peak values increase with the number of nucleation sites. This observation, coupled with the spatial distribution of stress shown in Figure 9 indicate that there is an enhanced possibility of fracture and voiding in the vicinity of nucleation sites.

Further, we see a convergent behavior with respect to number of nucleation sites, clearly shown by the similar peak stress profiles for the cases of two, three and five heterogeneous nucleation sites. This indicates that in the limit of numerous randomly located heterogeneous nucleation sites, one could obtain a characteristic peak stress distribution profile for a given electrode particle geometry and lithiation/delithiation cycles. Some details of our observations, such as the stress profile and its time evolution, may be influenced strongly by the specific lithiation/delithiation cycles chosen in these simulations, the circular geometry, and other IBVP-specific choices. However, the broad qualitative trends of the peak stress profiles and their time evolution can be very important to the design of lithiation/delithiation cycles and the understanding of the electrode particle performance and possible degradation.

Figure 9: Time snapshots of maximum principal stress distribution during the lithiation/delithiation cycles. Shown are the stress fields at time , , and for the problems with no nuclei, two nuclei and five nuclei. The unit of stress in the contour plots is GPa.
Figure 10: Variations in the peak values of the maximum principal stress during the lithiation/delithiation cycles in the presence of 0, 2, 3 and 5 nucleation sites. The flux variation indicating lithiation/delithiation cycles is superposed for reference.

Iv Conclusion

The atomistic calculations on the LiTiO system reported here add to the calculations reported in Ref. 8, and supply the necessary quantities used in the continuum scale studies. Here we have reported the elastic coefficients C, C and C of LiTiO , as well as the Gibbs free energy and the interfacial energy between the and the phases of LiTiO at 800 K. The calculated mobility indicates that transport in LiTiO is diffusion-limited. This same approach could be implemented, and perhaps automated, for a wide array of electrode materials in order to provide the groundwork for analogous simulation studies. This may encourage wider application of these methods to study battery electrode materials.

For the continuum scale studies, we considered a phase field framework for modeling the chemo-mechanical problem of phase transformations coupled with mechanics. The numerical framework is based on a coupled phase field and finite strain mechanics formulation. The latter choice distinguishes this work as traditionally phase field computations are limited by the infinitesimal strain approximation.

The continuum scale simulations provide insights to the kinetics (spatial distribution of Lithium ion composition) and mechanics (stress profiles and peak stress value evolution) and how these are affected by the lithiation/delithiation processes. Specifically, the peak stress profile is observed to be associated with the presence and the density of nucleation sites and its time evolution is associated with the rates of the imposed lithiation/delithiation cycles. As discussed, these stress profiles and peak stress values are important for understanding fracture and voiding.

Thus, continuum scale studies, informed by first principles, MMC and KMC studies, have a potential role to play in optimizing lithiation/delithiation cycles and predicting the effects of nucleation sites on the overall performance and mechanical degradation of electrode particles. First principles parameterizations of continuum simulations allow us to better understand battery performance via the simulation of lithiation/delithiation dynamics and studying the resulting Lithium ion diffusion kinetics and mechanical response, and thus potentially contribute to the design optimization of Lithium ion batteries.


This research was supported by the NSF CDI Type I grant: CHE1027729 “Meta-Codes for Computational Kinetics”. First principles calculations were performed on the Homewood High Performance Computing Cluster supported in part by the U.S. National Science Foundation, grant NSF-OCI-108849. The continuum scale computations were performed on a PRedictive Integrated Structural Materials Science (PRISMS) Center computing cluster at the University of Michigan. PRISMS is supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Materials Sciences and Engineering under Award #DE-SC0008637. AR was supported by NSF DUE1237992.


  • Hautier et al. (2011) G. Hautier, A. Jain, S. P. Ong, B. Kang, C. Moore, R. Doe, and G. Ceder, Chem. Mater. 23, 3495 (2011).
  • Cava et al. (1984) R. Cava, D. Murphy, S. Zahurak, A. Santoro, and R. Roth, J. Solid State Chem. 53, 64 (1984).
  • Zachau-Christiansen et al. (1988) B. Zachau-Christiansen, K. West, T. Jacobsen, and S. Atlung, Solid State Ionics 28, 1176 (1988).
  • Krtil and Fattakhova (2001) P. Krtil and D. Fattakhova, J. Electrochem. Soc. 148, A1045 (2001).
  • Colbow et al. (1989) K. Colbow, J. Dahn, and R. Haering, J. Power Sources 26, 397 (1989).
  • Wang et al. (1999) G. Wang, D. Bradhurst, S. Dou, and H. Liu, J. Power Sources 83, 156 (1999).
  • Wagemaker et al. (2005) M. Wagemaker, A. Van Der Ven, D. Morgan, G. Ceder, F. Mulder, and G. Kearley, Chem. Phys. 317, 130 (2005).
  • Bhattacharya and Van der Ven (2010) J. Bhattacharya and A. Van der Ven, Phys. Rev. B 81, 104304 (2010).
  • Van der Ven et al. (2001) A. Van der Ven, G. Ceder, M. Asta, and P. Tepesch, Phys. Rev. B 64, 184307 (2001).
  • Ouyang et al. (2004) C. Ouyang, S. Shi, Z. Wang, X. Huang, and L. Chen, Phys. Rev. B 69, 104303 (2004).
  • Morgan et al. (2004) D. Morgan, A. Van der Ven, and G. Ceder, Electrochem. Solid State Lett. 7, A30 (2004), URL
  • Kresse and Furthmüller (1996a) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996a).
  • Kresse and Furthmüller (1996b) G. Kresse and J. Furthmüller, Comput. Mater. Sci. 6, 15 (1996b).
  • Le Page and Saxe (2002) Y. Le Page and P. Saxe, Phys. Rev. B 65, 104104 (2002), URL
  • Liu et al. (2005) J. Z. Liu, A. Van de Walle, G. Ghosh, and M. Asta, Phys. Rev. B 72, 144109 (2005).
  • Sanchez et al. (1984) J. M. Sanchez, F. Ducastelle, and D. Gratias, Physica A 128, 334 (1984).
  • Laks et al. (1992) D. B. Laks, L. G. Ferreira, S. Froyen, and A. Zunger, Phys. Rev. B 46, 12587 (1992).
  • Fontaine (1994) D. D. Fontaine (Academic Press, 1994), vol. 47 of Solid State Physics, pp. 33 – 176, URL
  • Ceder et al. (2000) G. Ceder, A. Van der Ven, C. Marianetti, and D. Morgan, Modell. Simul. Mater. Sci. Eng. 8, 311 (2000).
  • Wolverton and Zunger (1998) C. Wolverton and A. Zunger, Phys. Rev. Lett. 81, 606 (1998).
  • van de Walle and Ceder (2002) A. van de Walle and G. Ceder, J. Phase Equilib. 23, 348 (2002).
  • Hart et al. (2005) G. L. W. Hart, V. Blum, M. J. Walorski, and A. Zunger, Nat. Mater. 4, 391 (2005).
  • Mueller and Ceder (2009) T. Mueller and G. Ceder, Phys. Rev. B 80, 024103 (2009).
  • Nelson et al. (2013) L. J. Nelson, G. L. Hart, F. Zhou, and V. Ozoliņš, Phys. Rev. B 87, 035125 (2013).
  • Binder et al. (2011) K. Binder, B. Block, S. K. Das, P. Virnau, and D. Winter, J. Stat. Phys. 144, 690 (2011).
  • Frenkel et al. (1997) D. Frenkel, B. Smit, and M. A. Ratner, Understanding molecular simulation: from algorithms to applications (Academic Press, Orlando, 1997).
  • Richards (1977) P. M. Richards, Phys. Rev. B 16, 1393 (1977), URL
  • Cahn and Hilliard (1958) J. W. Cahn and J. E. Hilliard, J. Chem. Phys. 28, 258 (1958).
  • Allen and Cahn (1972) S. M. Allen and J. W. Cahn, Acta Metall. 20, 423 (1972).
  • Deng et al. (2013) J. Deng, G. J. Wagner, and R. P. Muller, J. Electrochem. Soc. 160, A487 (2013).
  • Kao et al. (2010) Y.-H. Kao, M. Tang, N. Meethong, J. Bai, W. C. Carter, and Y.-M. Chiang, Chem. Mater. 22, 5845 (2010).
  • Burch and Bazant (2009) D. Burch and M. Z. Bazant, Nano Lett. 9, 3795 (2009).
  • Hu et al. (2013) S. Hu, Y. Li, K. M. Rosso, and M. L. Sushko, J. Phys. Chem. C (2013).
  • Huttin and Kamlah (2012) M. Huttin and M. Kamlah, Appl. Phys. Lett. 101, 133902 (2012).
  • Vaithyanathan et al. (2002) V. Vaithyanathan, C. Wolverton, and L. Chen, Phys. Rev. Lett. 88, 125503 (2002).
  • Bangerth et al. (2007) W. Bangerth, R. Hartmann, and G. Kanschat, ACM Trans. Math. Softw. 33, 24/1 (2007).
  • Heroux et al. (2005) M. A. Heroux, R. A. Bartlett, V. E. Howle, R. J. Hoekstra, J. J. Hu, T. G. Kolda, R. B. Lehoucq, K. R. Long, R. P. Pawlowski, E. T. Phipps, et al., ACM Trans. Math. Softw. 31, 397 (2005), ISSN 0098-3500.
  • Rudraraju et al. (2014) S. Rudraraju, A. Van der Ven, and K. Garikipati, Comput. Methods Appl. Mech. Eng. 278, 705 (2014), ISSN 0045-7825, URL
  • Rudraraju et al. (2016) S. Rudraraju, A. Van der Ven, and K. Garikipati, npj Comput. Mater. 2, 16012 (2016), URL
  • Hu and Chen (2001) S. Hu and L. Chen, Acta Mater. 49, 1879 (2001).
  • Seko et al. (2009) A. Seko, Y. Koyama, and I. Tanaka, Phys. Rev. B 80, 165122 (2009).
  • Jiang (2013) T. Jiang, A multi-physics study of lithium ion battery electrode materials, Ph. D. Thesis, Johns Hopkins University, MD (2013), URL
  • Sebastian and Gopalakrishnan (2003) L. Sebastian and J. Gopalakrishnan, J. Mater. Chem. 13, 433 (2003).
  • Van der Ven et al. (2013) A. Van der Ven, J. Bhattacharya, and A. A. Belak, Acc. Chem. Res. 46, 1216 (2013).
  • Zhu and Wang (2010) Y. Zhu and C. Wang, J. Phys. Chem. C 114, 2830 (2010).
  • Han et al. (2004) B. Han, A. Van der Ven, D. Morgan, and G. Ceder, Electrochim. Acta 49, 4691 (2004).
  • Funabiki et al. (1999) A. Funabiki, M. Inaba, T. Abe, and Z. Ogumi, J. Electrochem. Soc. 146, 2443 (1999).
  • Li et al. (2012) D. Li, P. He, H. Li, and H. Zhou, Phys. Chem. Chem. Phys. 14, 9086 (2012), URL
  • Bae and Pyun (1996) J.-S. Bae and S.-I. Pyun, Solid State Ionics 90, 251 (1996).
  • Choi et al. (1998) Y.-M. Choi, S.-I. Pyun, and J. M. Paulsen, Electrochim. Acta 44, 623 (1998).
  • Hong et al. (2006) J. Hong, C. Wang, and U. Kasavajjula, J. Power Sources 162, 1289 (2006).
  • Klinsmann et al. (2016) M. Klinsmann, D. Rosato, M. Kamlah, and R. M. McMeeking, J. Mech. Phys. Solids 92, 313 (2016), ISSN 0022-5096, URL
  • Miehe et al. (2016) C. Miehe, H. Dal, L.-M. Schänzel, and A. Raina, Int. J. Numer. Methods Eng. 106, 683 (2016), ISSN 1097-0207, URL
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description