Size effects in nanostructured Li-ion battery cathode particles

Cathode materials for Li-ion batteries exhibit volume expansions on the order of 10% upon maximum lithium insertion. As a result internal stresses are produced and after continuous electrochemical cycling damage accumulates, which contributes to their failure. Battery developers resort to using smaller particle sizes in order to limit damage and somemodels have beendeveloped to capture the effect of particle size on damage. In this paper,wepresent a gradient elasticity framework,which couples the mechanical equilibrium equations with Li-ion diffusion and allows the Young’s modulus to be a function of Li-ion concentration. As the constitutive equation involves higher order gradient terms, the conventional finite element method is not suitable, while, the two-way coupling necessitates the need for higher order shape functions. In this study, we employ B-spline functions with the framework of the iso-geometric analysis for the spatial discretization. The effect of the internal characteristic length on the concentration evolution and the hydrostatic stresses is studied. It is observed that the stress amplitude is significantly affected by the internal length, however, using either a constant Young’s modulus or a concentration dependent one yields similar results.


Introduction
Li-ion batteries are the most promising energy storage devices for portable devices and electric vehicles. The stability of commercial batteries is at least up to 1000 cycles, however, there still remain open questions that need to be addressed in order to improve their life-time and capacity. One of these issues is the volume expansions that the electrodes experience during Li-ion insertion and subsequent contractions during de-insertion. Such volume changes results in internal stresses that lead to damage. This problem mainly concerns failure in anode materials such as Si and Sn as they experience over 300% expansions and fracture occurs from the first few electrochemical cycles. Cathode materials such as LiNi 1/3 Co 1/3 Mn 1/3 O 2 , LiMn 2 O 4 , LiFePO 4 for LIBs, undergo a significantly lower volume expansion of approximately 10 %, resulting in damage accumulation over tens of cycles and contributing to the capacity retention over long term cycling. In [1,2] it was experimentally shown that continuous electrochemical cycling resulted in damage/fracture on the surface of LiMn 2 O 4 single crystal particles. Battery developers are exploring the use of nanoparticles as damage mechanisms are less severe at such scales, however, to understand the interplay between the particle size and performance, it is necessary to also perform modeling. Initial studies on the electro-mechanical behavior of cathodes [3,4], accounted for the diffusion of Li within classical mechanics, which cannot account for the particle size. Linear elastic fracture mechanics were also used in [2,5] in order to develop design criteria for materials selections and configuration that increase mechanical stability of anodes during cycling. Furthermore, several aspects of stress diffusion coupling, both one way and two way have been studied [1,2,5], however, a comprehensive numerical investigation has not been carried out to understand the effect of the size of the electrode particles. Moveover, recent first principle calculations suggest that the Young's modulus of cathode materials is dependent on the Li-concentration, and therefore a concentration dependent modulus needs to be accounted for in order to accurately capture intercalation induced mechanical degradation. In this paper, we use a fully coupled chemo-mechanical framework to study the size effects and the effects of stress-diffusion interactions in cathode materials of LIBs. To do so the framework of gradient elasticity (GRADELA) as put forth by Aifantis [6][7][8] is employed, which is a non-local theory that enhances the classical linear elasticity with only one extra material parameter (the internal length scale, ℓo) that multiplies the Laplacian of the Hookean stress, which is then added to the classical linear elasticity expression [6]. The resulting equations are numerically solved using the Galerkin approach. The addition of an extra material parameter increases the continuity requirement of the functions used to approximate the unknown fields. This is accomplished by employing the non-uniform rational B-splines for the spatial discretization.
The paper is organized as follows: Section 2 presents the governing differential equations for the diffusion of Liions and the corresponding gradient elasticity equations. Section 2.1 presents the weak form for the coupled differential equations and discusses the spatial and temporal discretization. The influence of various parameters such as the particle size, diffusion induced stress, stress induced diffusion, internal length scale on the concentration build up and the stresses in the particle are systematically studied in Section 3, followed by concluding remarks in the last section.

Theoretical formulation
In this study, we assume that the response of a particle that is being lithiated is governed by the fully coupled chemo-mechanical problem involving two fields: the Liion concentration and the elastic displacements. Consider a linear elastic body occupying an open domain B ⊂ R 2 and bounded by a surface with an unit outward normal, n. Let (u, c) : B → R 2 be the displacement and the concentration at a point x of the elastic body when it is subjected to external tractions t : t → R 2 and external flux J : J → R 2 . The boundary is assumed to be decomposed into: The boundary value problem for the coupled chemo-mechanical problem in the absence of body forces (in elasticity) and the source terms (in the diffusion equation) is: subject to the following boundary conditions: Here u is the displacement vector, is the Cauchy stress tensor, c represents the concentration and J represents the diffusive flux. The diffusive flux is given by the gradient of the chemical potential as: where M is the mobility of the lithium ions and the chemical potential in an ideal solid solution is expressed as: where o is is the chemical potential of a species under a known pressure, R is a gas constant, T is the absolute temperature, X is the molar fraction of lithium ions, is the partial molar volume of lithium ions and h is the hydrostatic stress. Using Equation (4), the flux can be written as: where D = MRT is the diffusion constant and h = 1 3 tr( ) is the hydrostatic stress. The influence of the concentration on the stress field can be modelled by modifying the stress field within the gradient elasticity [6] by thermal analogy (similarly as in the recent works of [9,10]) by: where C = f (E, ) is the material constitutive matrix, E, are the Young's modulus and the Poisson's ratio, respectively, is the infinitesimal strain tensor, c = 3 c is the strain due to concentration, c = c − c i is the change in the concentration and c i is the initial concentration, ℓo is the internal characteristic length.
In this study, we are interested in understanding the stress and the concentration distribution in a spherical particle. Hence, owing to spherical symmetry, the above equations become one-dimensional and are given by: where r and t are the stress in the radial and the tangential directions. The stress-strain relations within the GRADELA framework [6] are It has been observed experimentally that the Young's modulus electrode materials depends on the Li-ion concentration [11]. In this work, we assume a linear dependence: where Eo is the Young's modulus of the material, which is independent of the concentration, prior to Li-insertion, cm is the maximum Li concentration in the particle and k is a material constant. The flux and the concentration are related by: and the strain-displacement relations are given by:

Weak form, spatial and temporal discretiztion
The two-way coupled interaction between the mechanical degradation and the diffusion process is governed by Equation (2). In this section, we develop the corresponding weak form of the governing equations and the solution methodology.

Weak form
The weak form is obtained by invoking the Galerkin orthogonality, wherein, the error due to the assumed test spaces for the displacement and the concentration are perpendicular to the trial spaces within the domain. Mathematically, the weak formulation is: find u ∈ U and c ∈ C such that ∀v ∈ V : where, U and C are the displacement and the concentration trial spaces, respectively and V is the test space: where the space W( ) includes linear displacement fields and the concentration field. The domain is partitioned into elements h , and on using shape functions a that span at least the linear space, we substitute vector-valued trial and test functions (12) and apply a standard Galerkin procedure to obtain the following mechanical and the chemical internal force vector, P int = {P u , P c } T : B u is the strain-displacement matrix, computed using the derivatives of the shape functions. The external nodal forces arising from the boundary tractions and the flux are considered in the global residual vector. For temporal discretization, we employ the backward Euler time integration scheme. The damping matrix is given by: All other components of the damping matrix are zero. The Galerkin procedure leads to a coupled non-linear system of equations (see Equation (14)). Next, we describe the linearization procedure based on the Taylor's series expansion of the residual vector. Let d = {u, c} T denote the vector of unknown displacements and concentration. Let d i+1 = d i + d be the solution at the (i+1) th iteration that depends on the solution at the i th step and a small increment d. The residue at the (i + 1) th step is: Using Taylor's expansion and upon eliminating the higher order terms, we get the following linearized form for the above equation as: where, T is the tangent stiffness matrix, given by: In the present study, we consider only Dirichlet boundary conditions and hence the derivative of the force term is zero. The tangent stiffness matrix is given by: The resulting nonlinear equations are solved using the Newton-Raphson method.

Spatial and temporal discretization
As noted earlier, the presence of the coupling between the diffusion and the stress and the gradient elasticity framework requires a high-order continuous function to discretize the unknown fields, i.e., the displacement ur and the concentration c. The GRADELA framework introduces the Laplacian of the strain components, which requires a C 1 continuous shape functions. One approach to solve this problem is to employ the commuting property of the operators as discussed in [12].
In this study, we employ the non-uniform rational B-spline (NURBS) basis functions within the framework of isogeometric analysis to represent the unknown fields. This is because the continuity of the basis functions can be adapted to suit the geometric modeling and during the discretization. Thus, the single set of basis functions are used to represent both the displacement field and the concentration within the spherical particle. We give here only a brief introduction to NURBS. More details on their use in FEM are given in [13,14]. Within this framework, the unknown field, (u, c) is approximated by: where u h and c h are the nodal variables and R( ) are the NURBS basis functions. The key ingredients in the construction of NURBS basis functions are: the knot vector (a non decreasing sequence of parameter values, i ≤ i+1 , i = 0, 1, · · · , m − 1), the control points, P i , the degree of the curve p and the weight associated to a control point, w. A p th degree NURBS basis function is defined as follows: where w i are the weights for the i th B-spline basis function N i,p ( ). The i th B-spline basis function of degree p, denoted by N i,p is defined as [15]: The B-spline basis functions has the following properties: (i) non-negativity, (ii) partition of unity, ∑︀ interpolatory at the end points. As the same function is also used to represent the geometry, the exact representation of the geometry is preserved. It should be noted that the continuity of the spline functions can be tailored to the needs of the problem. Moreover, the spline function has limited support. When employed to approximate the FE solution space, the resulting stiffness matrix has similar properties to the stiffness matrix computed by employing Lagrange shape functions. Figure 1 shows the third order NURBS for an open knot vector = {0, 0, 0, 0, 1/3, 1/3, 1/3, 1/2, 2/3, 1, 1, 1, 1}.

Results
In this section, the effect of the stress-induced diffusion and the diffusion induced stress in a spherical particle is numerically studied using the isogeometric analysis framework. The effect of the dependence of the Young's modulus on the concentration and the internal characteristic length is also studied. The internal length, modelled within the GRADELA framework, requires boundary conditions for the higher stress order terms, which in this study are assumed to be zero on the boundary. The domain is assumed to be subjected only to diffusion boundary conditions and only 'essential' boundary conditions for elasticity are enforced to yield a non-singular system of equations. For the diffusion type boundary condition, only a potentiostatic boundary condition is considered. For the numerical study, both the domain and the unknown field variables are discretized with NURBS basis functions with order p = 3, this allows us to account for the higher order terms due to the GRADELA framework. However, one could choose higher order continuous functions, i.e. p > 3 using NURBS basis functions. Based on successive refinement, a total of 100 NURBS quadratic elements were found to yield accurate solutions. The NURBS basis functions are rational functions, hence conventional numerical quadrature rules are less efficient. We adopt a nearly optimal quadrature rule proposed in [16] to integrate the terms in the stiffness and mass matrices.
The material properties for the Mn 2 O 4 Cathode material are [17]: Young's modulus, Eo = 10 GPa, partial molar volume = 3.497 ×10 −6 m 3 mol −1 , diffusion coefficient D = 7.08×10 −15 m 2 /s and the maximum concentration cm = 2.29×10 4 mol m −3 . The Young's modulus is assumed to be a function of concentration with k = 1 GPa [18] (c.f. Equation (9)) and c i is the initial Li-ion concentration, which without loosing generality is assumed to be zero, corresponding to the beginning of lithiation. Before proceeding with the analysis, the governing equations are normalized such that the nonlinear iteration scheme does not face any convergence issues. The possible convergence issues could be attributed to the different orders of magnitude of the physical parameters of the chemical and the mechanical properties. The governing equations are normalized as [19]: where is the characteristic length scale of the domain considered for the numerical simulation, r, u and t are the spatial coordinate, displacement and time with '*' representing dimensionless quantities. The material parame-ters are scaled by /RT such that the stress field becomes dimensionless. Due to symmetry, the particle is modelled as being a one-dimensional problem with the following boundary conditions: where cm is the maximum concentration and a is the radius of the particle. The boundary at r = a is maintained at cm for all time steps and the simulation is stopped when the particle is fully lithiated, i.e. c(r, t) = cm. Two different particle sizes are considered of radius a = 20 m and 50 m, therefore a full electrochemical cycle is not considered but only the first maximum lithiation. Since small strain deformation is considered, the volume changes are not significatn and the radius (r = a) is not affected by the Li-ion concentration. Figure 2 shows the evolution of the normalized concentration, normalized radial displacement, normalized radial stress and normalized hydrostatic stress at the center of the particle (r * = 0.5). The influence of the concentration dependent Young's modulus and internal length scale is examined. To study the effect of the internal length scale, two different values for ℓo are considered (= 0, 0.4), with ℓo = 0, reducing the equations to those of classical elasticity. Furthermore, for comparison purposes, the case where the Young's modulus is a constant and not dependent of on the Li-ion concentration is also plotted. It is observed that, using either a constant or concentration dependent Young's modulus yields similar results. The gradient elasticity solutions (i.e. when ℓo is different from zero) yield higher stress and concentration values at the particle center.
Next, the influence of the particle size, the internal length and the effect of the concentration dependent Young's modulus on the evolution of the concentration and the hydrostatic stress at the center of the particle are shown in Figure 3 for two different particle sizes of radius a = 20 m and 50 m. For comparison, the results are shown only until 500 min. This is because the larger particle takes more time to reach the maximum concentration. From the parameters examined, it is seen that the value of the internal length does not significantly affect the concentration profile at the center of the particles. Similarly, considering a concentration dependent elastic modulus gives a similar concentration and hydostatic stress profiles, as when a constant modulus was assumed during lithiation. The most significant effect is that of the particle size on the hydrostatic stress, which is also significantly affected by the value of the internal length, but up to 100 mins of lithiation. It is seen that for the smaller particles this stress  drops to zero after about 340 min, while the larger ones experience significant stresses even beyond 500 min. Table 1 shows the effect of the internal length and the concentration dependent Young's modulus on the maximum hydrostatic stress in the particle with a radius of a = 20 m. It can be seen that the internal length, ℓo has a signif-icant influence on the maximum hydrostatic stress in the particle whilst the concentration dependent Young's modulus has minimum impact. This is because although the Young's modulus is a function of concentration, the rate of change of Young's modulus dE/dc is small. Hence, we do not see an appreciable influence of the concentration dependent Young's modulus on the maximum hydrostatic stress. It should be noted here that a difference in response between the 20 and 50 micron particles is observed even in the case of ℓo = 0. This is because, the larger the particle, the slower is the rate at which the lithium ion concentration and the hydrostatic stress increases inside the particle. Hence, the size effects are obtained due to the coupling with concentration.

Conclusions
The present article employed gradient elasticity coupled with stress assisted diffusion within the isogeometric analysis in order to capture the effect that the particle size has on the maximum value of the hydrostatic stress at the center of the spherical particle during lithiation. Comparing with the case of classical elasticity, it is seen that the concentration profiles are similar either with or without gradient considerations, however, the hydrostatic stress was significantly higher in the case of gradient elasticity. As experimental observations have shown that the elastic modulus is dependent on the Li-ion concentration, the case of a concentration dependent elastic modulus was considered, however, similar concentration and hydrostatic stress profiles were obtained as when a constant modulus value was assumed. This can justify the use of constant moduli in the implementation of more complex frameworks, that account for damage and fracture within the active material.
tion for supporting this work through the CMMI grant (CMI-1762602).