The Impact of Compositionally Induced Residual Stress on Electrochemical Shock in Battery Electrode Particles

Large inhomogeneous volume changes can induce severe mechanical stresses in Li-ion battery electrode particles under high charging/discharging rates, thus causing failure and capacity fading. Varying the material composition of the particle in a speciﬁed way can impose residual stresses that reduce the resultant stress levels. Here, we model a LiMn 2-x Ni x O 4 spherical cathode particle with a semicircular surface crack using ﬁnite elements method on the basis of linear elastic fracture mechanics. We also present the effect of residual stresses on the resultant diffusion induced stresses and fracture behavior of the particle

Transition-metal oxides, such as LiCoO 2 and LiMn 2 O 4 , are currently the most widely used cathode materials for lithium-ion batteries because of their high specific capacity and low cost. 1,2 A large amount of research has been devoted to improving the capacity and lifetime of these materials. For example, reducing particle sizes can enhance battery performance. 3,4 Other approaches to improve the calendar life and capacity include novel structures [5][6][7] or composites, 8 and the use of various compositions. [9][10][11][12][13][14][15] Developing a graded material by varying the composition through the structure can also improve the performance. 16,17 In this paper, we consider the latter option, based on the concept of utilizing a graded material to improve the performance of an electrode particle on a mechanical basis.
Capacity fading mechanisms in lithium-ion batteries are not fully understood, largely because they appear to be influenced by a large number of different phenomena. 18 Lithium ion insertion/extraction induces nonuniform volume changes during the cycling process, resulting in mechanical stresses in the particle. Resolving the mechanical behavior of the electrode particle can help understanding these phenomena. Utilizing linear elastic fracture mechanics (LEFM), Aifantis and Dempsey obtained stability diagrams for various configurations of nanocomposite electrodes to predict the extent of stable crack growth of the radial cracks at the matrix-particle interface. 19 Christensen and Newman modeled the deformation of spherical LiMn 2 O 4 cathode particles with Eulerian strain description and predicted the onset of fracture using tensile stress criterion. 20 Zhang et al. carried out finite difference simulations of spherical particles and finite element simulations of ellipsoidal particles to study intercalation induced stresses. 21 Using LiMn 2 O 4 as the cathode material, they reported that it is desirable to synthesize electrode particles with smaller sizes and larger aspect ratios to reduce intercalation induced stresses. A more recent study used the extended finite element method based fracture analysis to investigate the effects of current density, particle size and aspect ratio on fracture propagation. 22 Cheng and Verbrugge established a criterion for crack propagation within a spherical electrode using the stored total elastic energy. 23 Zhao et al. constructed a fracture map for a LiCoO 2 particle using dimensionless parameters for diffusion kinetics and fracture mechanics. 24 Woodford et al. derived a fracture mechanics failure criterion for the electrode particles and produced an electrochemical shock map that shows regimes of failure depending on charge rate, particle size, and the material's fracture toughness. 25 They demonstrated galvanostatic charging of a spherical LiMn 2 O 4 particle and indicated that small surface flaws can grow unstably under certain conditions. The current study investigates the influence of residual stresses, which are generated by concentration gradient of the components, on fracture behavior of spherical electrode particles. LiMn 2-x Ni x O 4 is selected as the cathode particle material and the concentration distri-bution of Ni in the particle is represented by a given anticipated sample function. This analysis focuses on stresses induced by variations in the Ni/Mn ratio, where complexities due to phase transformations are intentionally neglected in this initial treatment. Single phase cubic spinels have been synthesized for compositions up to x = 0.5. 9,15 It is well documented that lithiation of these phases can result in significant two phase behavior with significant miscibility gaps, however, recent work also shows that disordered cation arrangements can lead to single phase solid solutions during lithiation/delithiation cycles. 26 This type of single phase material is used as the basis for our model. This type of compositional tempering approach has been used in a variety of other glass and ceramics applications, to strengthen the materials by generating residual compressive stresses at the surface and suppressing crack-growth. For example, ion exchange is a compositional tempering process used to strengthen glasses, where small alkali ions near the surface are replaced by larger ones that are introduced from a molten salt solution. 27 This approach can be quite effective when fracture initiates at surface flaws. However, this type of compressive stress at the surface also induces a compensating tensile region within the body. Thus, if the crack tip is within this tensile region or if the body has an internal flaw, the tempering can also promote fracture. 28,29 This distinction also applies to the approach that we analyze in this paper, in that inducing compositional stresses near the particle surface should mitigate fracture that is initiated at surface flaws. The compensating tensile stresses inside of these particles could promote fracture if there are critical flaws here instead.
The coupled field problem of the particle is solved by using a two dimensional (2D) finite elements method (FEM) code developed in MATLAB. During the lithium intercalation/deintercalation (discharge/charge) process, the volume change for the cathode materials stays in a range (∼10%) that the assumptions of linear elasticity can be made without loss of generality. However the stress-diffusion coupling results in a nonlinear equation system, which is handled by the Newton Raphson scheme. The fracture behavior of the particle is analyzed with the LEFM theory as in the method proposed by Woodford et al. 25

Governing Equations
Stress equilibrium.-The static equilibrium equation for the body in terms of stress σ and body force b fields is For an isotropic linear elastic body, stress-strain relation regarding the volumetric expansion induced by diffusion can be written as where c is the concentration change of the diffusion species from a reference value, and is the partial molar volume of the solute. Here, c r and r are the predefined concentration and corresponding partial molar volume of the component, which generates the residual Journal of The Electrochemical Society, 162 (7) A1282-A1288 (2015) A1283 stresses. The strain field ε in terms of displacement field u, and the elastic stiffness tensor C is defined as where E is Young's modulus and ν is Poisson's ratio.
Diffusion equation.-Conservation of species gives the following diffusion equation The species flux J can be written as 21,30 where M is the mobility of species, and μ is the chemical potential. The modified chemical potential in an ideal solid solution can be defined as where μ 0 is a constant, R is gas constant, T is absolute temperature, X is the molar fraction of lithium ion, and σ h is the hydrostatic stress (σ h = σ ii /3). The last term of this equation can be defined as a stress contribution to the potential: The materials of interest here are not expected to exhibit ideal solution behavior, and thus Eq. 6 is an oversimplification. However this approach has been used by a number of other researchers, 21,25 largely because accurate descriptions of the excess free energy of mixing are not available for these materials.
Assuming temperature is uniform, and substituting Eq. (6-8) into Eq. 5, species flux is obtained as where D = M RT is the diffusion coefficient. Substituting Eq. 9 into 4 the diffusion equation will be

Solution Process
The nonlinear stress-diffusion coupling equation system, which consists of the stress equilibrium Eq. 1 and diffusion Eq. 10, is solved by using FEM with a Newton Raphson scheme. Details about application of the FEM can be found in standard references. 31,32 Discretization in time is accomplished with the Crank-Nicolson method. The stress potential Eq. 7 is also satisfied independently. In the formulations, note that the independent variables are the nodal displacements u, concentrations c and the stress potentials μ s .
Defining a kinematically admissible virtual displacement field δv, the equation of stress equilibrium is forced by the following virtual work equation where t denotes the traction field.
The displacement, virtual displacement and concentration fields can be interpolated by using same shape functions [12] Substituting Eqs. 2, 3 and 12 into Eq. 11 and applying the integration by parts gives the residual In matrix form, Eq. 13 is K aibk The stiffness matrices and the force vectors are given in Appendix A.
Similarly, the FEM discretization can be applied to diffusion equation 10. FEM equations can be obtained by using the Galerkin method.
Approximation functions for the FE discretization of the concentration field, corresponding weight functions and the stress potential field are given as [15] Applying integration by parts to Eq. 14 and after some manipulations Using the Crank-Nicolson method, the diffusion equation is discretized in time as follows where t, n and θ denote time-step, time-step number and implicitness parameter. Here, θ = 1, and this constitutes a fully implicit scheme. For the stress potential equation 7, FEM discretization can be obtained using trial functions [18] and the residual vector will be For each time step, the linear equation system in terms of increments obtained from the consistent tangents of residual equations will be satisfied in the Newton Raphson scheme (see Appendix A). A more general and detailed finite element model for electrode materials can be found in Ref. 33.

Numerical Example
Constant diffusion flux at the surface of a LiMn 2-x Ni x O 4 spherical particle.-Fast charging in cathode particles can create high Li concentration gradients, leading to diffusion induced tensile stresses and fracture. This problem, which is described as an electrochemical shock, 25 will be considered as an example in this work.
The spherical particle under constant surface flux is modeled using 4 node quadrilateral finite elements. The FEM mesh and the scheme of the problem are shown in Fig. 1. Considering the axisymmetry of the problem, only the first quadrant is discretized by 500 elements with 541 nodes.
The lithium concentration flux given in Eq. 9 can be defined in terms of current density i n as J = i n F where F is the Faraday's constant. The relationship between the current density and charge rate (C-Rate) for a spherical particle is αρr max 3 The boundary and initial conditions for the axisymmetric problem can be written as displacements in x and y directions; u x | x=0 = 0, u y | y=0 = 0 The dimensionless variables in the formulations arē The material properties for the spherical particle are listed in Table I.
The end of the galvanostatic charging process (deintercalation) is chosen as a critical state that the Li concentration reaches zero at the surface of the particle. In this case, the high charging rate results in a steep concentration gradient, which induces a tensile stress field at the outer region of the particle. Since the cathode material is brittle, this stress field can reduce the fracture strength by opening the surface flaws, causing unstable growth of surface cracks and leading to the failure of the particle. Thus, imposing residual compressive stress field at the outer region can enhance the fracture resistance of the particle by mitigating the tensile stress field. Substitution of Mn with a suitable transition metal through the radius in a gradual manner can induce residual compressive stresses at the outer region of the particle. Replacing Mn with Ni in LiMn 2 O 4 causes the unit cell volume to shrink, 9 and thus can be used as a mechanism to create the compressive zone. Hence, in our case Ni is selected as a substitutional element and its atomic ratio is allowed to vary with the radius according to a specified function. Assuming that the Ni substitution does not affect the material properties considerably, we treat the material properties as constants.
The partial molar volume r depending on the Ni concentration is calculated by using the unit cell volumes for LiMn 2  Gradual decrease of Ni concentration from center to the surface of the particle will create shrinkage in the middle and the outer region will be compressed to sustain the balance. Here, we propose a plausible approximate function for the concentration distribution of Ni as a function of dimensionless radial coordinate (see Fig. 2). We can thus represent the distribution function in our FEM code easily (see the Appendix B for the relationship with diffusion profile in a sphere).
In the following analysis, the value m = 10 is chosen as an initial example for numerical calculations. The concentration distribution  function corresponds to Ni concentrations that are maximized at the center and zero at the surface. To facilitate comparisons between these results and Ref. 25, the radius of the particle (r max ) is taken as 21 μm.
Compressive residual tangential stresses (σ θ ) of the particle along the radius for the two different maximum Ni concentrations (x max ) can be seen in Fig. 3. For a particle under galvanostatic charging flux ofJ = 0.92 or at C-Rate 5 (h −1 ), Li concentration and tangential stress profiles at the end of charge are obtained (see Fig. 4). In Fig. 4, the results for x max = 0 correspond to LiMn 2 O 4 particle, which does not have any residual stress effect, and can be compared to those of Ref. 25. The effect of the residual stresses on the concentration field is insignificant for this case. However, the residual stresses reduce the tensile tangential stress at the outer region considerably.
The failure criterion for the brittle particle can be determined by utilizing the LEFM in terms of stress intensity factor (SIF). It is assumed that, the crack will grow if the SIF for the crack exceeds a critical value (K IC ). This critical value, typically referred to as the fracture toughness, is a material property as a measure of resistance to extension of a crack. Here, we only consider the mode-I fracture toughness, which dictates failure for simple tensile loading. The SIF for a semicircular surface flaw on the particle can be obtained by using the corresponding tangential stress profiles along the radius   and employing the following weight function approach of Mattheck et al. 34 (see details in Ref. 25).
where K re f is the reference stress intensity factor, u r is the crack opening displacement field and a is the crack length. Stress intensity factors are calculated as a function of crack length for each different composition of the particle (see Fig. 5). As seen in the figure, the residual compressive tangential stresses at the outer region reduce the stress intensity factor particularly for small initial flaws. Increasing the maximum concentration of Ni flattens the SIF curves and reduces the maximum values of SIF considerably. Note that, the negative stress intensity factors imply crack closure, which is not specifically described in the method employed here. Hence, the results for these cases should be interpreted accordingly.
To the best of our knowledge, the fracture toughness of LiMn 2 O 4 is not known. Thus, 1 MPa m 1/2 is used here as an acceptable estimate. For an ideal brittle material, the fracture toughness is independent of the crack extension. Hence, if the crack growth driving force is an increasing function of the crack length, the crack will grow unstably. If it is a decreasing function, the crack will grow stably. 35 This stability condition for crack growth can be expressed in terms of SIF as This condition implies that the maximum point of the SIF curve determines the transition from unstable to stable region. Fig. 5 shows that increasing x max lowers the peak value of the SIF curve, while the peak position shifts to higher flaw sizes. The effect of the shape of the Ni concentration profile (determined by m) on the maximum SIF   Using Fig. 6a, we can also find the m value which corresponds to the lowest maximum SIF for a prescribed x max . Fig. 6b enables us to see how the peak position of the SIF curve shifts, for different x max values. These can help tailor the Ni concentration distribution based on fracture behavior. Furthermore, Fig. 7 shows the change in the SIF curves as m varies for a fixed maximum Ni concentration (x max = 0.5).
The maximum SIF values, which occur at the end of the galvanostatic charging for different C-Rates and particle sizes, are calculated. These values are shown as contours of specific fracture toughness values in an electrochemical shock map (see Fig. 8). This map defines the onset of fracture by using the K IC values as criteria. The left sides of the curves represent lower C-Rates, which result in SIF values lower than K IC, thus defining the safe zone. On the contrary, the right sides define the unsafe zone at which fracture occurs. The map in Fig. 8 is plotted on logarithmic axes for particle sizes 5, 6 . . . 10, 15 . . . 50 μm and C-Rates 1, 2 . . . 5, 10, 40 and 100. Contour curves for the specified K IC levels of 0.5, 1 and 2 M Pa √ m are calculated with the MATLAB 4 griddata method. Since we used a relatively small number of data points, the contour curves are not perfectly smooth. However, the shift of the safe zone to higher C-Rates is clearly seen in Fig. 8. According to this map, the critical C-Rate for a 20 μm particle having K IC level of 1 M Pa √ m can be increased from 1.5 to 6 h −1 by choosing x max = 0.3 and m = 4.
As noted in the Introduction, the compositional tempering strategy proposed here is only effective if fracture initiates at surface flaws. This is known to occur in many ceramic materials, however, if critical flaws are located inside of the particles then the tensile stresses here should promote particle fracture. If this is the case then the compositional tempering could be reversed to improve fracture resistance (i.e., by varying composition to create tensile stress near the surface and compressive stress inside of the particles).

Conclusions
During the fast charging/discharging process of an electrode particle, the tangential stresses at the surface can reach high levels due to the large inhomogeneous volume changes. These stresses can induce growth of the pre-existing surface flaws which can reduce the performance of the particles considerably. Varying the composition of the material can generate residual compressive stresses at the surface region and mitigate the crack opening effect of the diffusion induced stresses.
In this paper, we modeled a spherical LiMn 2-x Ni x O 4 cathode particle, which has a specified Ni concentration distribution along the radius. With this graded material concept, we showed that the SIF value of the surface cracks can be reduced, thus preventing crack growth. Hence, we conclude that the mechanical reliability of the particle under high charging rates can be improved by using chemical changes to induce residual compressive tangential stresses.

Acknowledgments
We gratefully acknowledge support from the National Science Foundation, award DMR-1410946 (Ceramics program). The authors would also like to acknowledge valuable input and comments from Prof. Allan F. Bower of Brown University, and thank two anonymous reviewers for their helpful comments and suggestions.

Appendix A: Finite Element Formulation
The stiffness matrices and the force vectors for the equations are given as The linearized system of equations is obtained by the application of Newton Raphson scheme to the residual vectors. where

Appendix B: Diffusion Profile in a Sphere
The well-known solution for one dimensional diffusion problem in a sphere 36 is given as where a is the radius of the sphere and c 0 is the surface concentration. The Ni concentration distribution can be related to this diffusion solution by taking c max = 1, c 0 = 0, a = 1. Hence, the following dimensionless equation can be obtained: We can find the corresponding diffusion profile that depends on Dt for each m value using least squares approximation (see Table BI).