A Hybrid Phase Field Model for Fracture Induced by Lithium Diffusion in Electrode Particles of Li-ion Batteries
HHighlights: • A hybrid PF model is proposed to study fracture behavior of electrode during lithiation. • The PF model is applied to a Si NW electrode. • The concentration-dependent elastic properties is taken into account. • Hybrid model shows less tendency to crack growth than isotropic model.1 a r X i v : . [ c s . C E ] M a r Hybrid Phase Field Model for Fracture Induced by LithiumDiffusion in Electrode Particles of Li-ion Batteries
Masoud Ahmadi
Faculty of Mechanical Engineering, University of Guilan, P.O. Box 3756, Rasht, Iran
Abstract
Lithium-ion batteries (LIBs) of high energy density and light-weight design, have found wide applica-tions in electronic devices and systems. Degradation mechanisms that caused by lithiation is a mainchallenging problem for LIBs with high capacity electrodes like silicon (Si), which eventually can re-duce the lifetime of batteries. In this paper, a hybrid phase field model (PFM) is proposed to studythe fracture behavior of LIB electrodes. The model considers the coupling effects between lithium(Li) -ion diffusion process, stress evolution and crack propagation. Also, the dependency of Elasticproperties on the concentration magnitude of Li-ion is considered. A numerical implementation basedon a MATLAB finite element (FE) code is elaborated. Then, the proposed hybrid PF approach isapplied to a Nanowire (NW) Si electrode particle. It is shown that the hybrid model shows lesstendency to crack growth than the isotropic model.
Keywords:
Phase field model, Hybrid formulations, Lithium-ion batteries, Crack propagation,Finite element method, Electrode particles.
Nomenclature : Double dot product¯ σ Largest principal value of the effective stresses¯ f Integral term of Fβ A constant scalar D Damping matrix K Tangential stiffness matrix R Nodal residual vector
Email address: [email protected] (Masoud Ahmadi)
Preprint submitted to Journal of Power Sources Friday 20 th March, 2020
Displacement field σ Cauchy stress tensor/vector ε Infinitesimal Strain tensor/vector ε c Chemical strain tensor ε e Elastic strain tensor B I φ Cartesian derivative matrix B I c Cartesian derivative matrix B I u Usual strain matrix D
2D elasticity matrix D A constant vector D A constant vector Q Field variables vector · Dot product χ Relaxation constant of the fracture order¨ (cid:3)
Second material time derivative∆ Increment δ Variation operator˙ (cid:3)
First material time derivative η Threshold stiffness parameter for a fully broken region η φ Test function for φη c Test function for cη u Test function for u γ A constant scalar (cid:104) (cid:3) (cid:105)
Macaulay brackets 3n Natural logarithm D Fourth-order stiffness tensor µ Standard chemical potential ∇ Nabla operator ν Poisson’s ratioΩ Partial molar volume ∂ Partial differential operator φ Fracture order field σ p Hydrostatic stress (cid:3) T Transpose of a matrix J Ion flux vectord Total differential operator ξ u Elastic energy density function ξ + u Tension part of the elastic energy density function ξ − u Compression part of the elastic energy density function a Initial crack length c Li concentration field c Initial Li concentration c max Maximum Li concentration D A constant scalar E Young’s modulus F Total free energy f φ Local free energy density of φf c Local free energy density of c u Local free energy density of u g (cid:48) ( φ ) Derivative of energetic degradation function with respect to φg ( φ ) Energetic degradation function G cr Critical energy release rate h Thickness of NW electrode I Node index i Newton iteration index k φ Gradient energy coefficient of φk B Boltzmann constant k c Gradient energy coefficient of cl Regularization constant M Molecular mobility m Number of nodes of each element n Time step index N I Nodal values of shape functions N A Avogadro’s constant R Radius of NW electrode T Absolute temperature t Time V Volume I Second-order unit tensor 5 . Introduction
Due to high energy storage density, LIBs are widely used in different technologies such as portableelectronic devices and electric vehicles [1]. The working principle of LIB cells lies essentially in theelectrochemical potential driven redox reaction in the electrode active materials [2]. In LIBs, Li-ions transfer from the anode and diffuse through the electrolyte towards the cathode during chargeand when the battery is discharged, the respective electrodes change roles. Abundant efforts havebeen made to develop next-generation battery systems with novel electrode materials, however, mostcandidates with promising electrochemical potential have chemo-mechanical stability problems [3].Due to Graphite low volume changes (about 10%) under intercalation, Graphite is a common anodeamong the commercial batteries [4]. Silicon is a promising candidate material with a theoreticalcapacity of 4200 mA h/g, which is exceedingly higher than the theoretical capacity of Graphite with372 mA h/g. Nevertheless, Si shows a volume changes of about 300%, which causes high stress inelectrodes that eventually lead to fracture, fragmentation, or pulverization [5, 6, 2]. Cracking ofelectrodes under diffusion is one of the main reasons for the short life span of LIBs with high capacityelectrodes [7].For a better understanding of the fracture behavior of LIBs during lithiation, numbers of model-ing and computational simulations have been expanded, which can offer insights into the structuralreliability of electrodes in LIBs. Ryu et al. [8] used a fracture mechanics approach to study fracturebehavior of Si NWs during lithiation/delithiation process. They considered large deformation associ-ated with lithiation and effects of pressure gradients on the diffusion of Li in their study. Grantab andShenoy [9] proposed a method that accounts for the effects of pressure-gradients within the materialon the flux for studying lithiation-induced crack propagation in Si NWs. Their model consists ofFE simulations in which the pressure-gradients are computed numerically to capture the effect of thecrack-tip on the localized diffusion. Chen et al. [10] established an analytical model to study the stressevolution and crack propagations in spherical particle electrodes during phase transformation. Hu etal. [11] developed an analytical model to study the diffusion induced stress evolution and discussed thecrack growth by using stress intensity factor coupled with surface effects. They showed that smaller SiSpherical particles exhibit higher structural integrity. Wang et al. used Peridynamic theory to studythe lithiation induced stress and fracture in Si square thin film electrodes [12, 13], and Si sphericaland cylindrical nanoparticles [14]. Gwak et al. [15] investigated the stress evolution and dynamiccracking behavior of Si NWs using a chemo-mechanical model based on the large deformation theoryand the bilinear cohesive zone model.In sharp-interface modeling of phase-separating behavior, a discontinuity is imposed at the inter-face, and a subsequent strain mismatch will give rise to stress concentration at the interface region [16].6hase field modeling as one of the most prominent approaches for the simulation of microstructureevolution is introduced as an alternative to sharp-interface modeling. PFM is based on the thermo-dynamic description of non-equilibrium states in materials including interfaces and is formulated asa state variable in space and time, the evolution of which controls the pathway towards equilibrium[17]. Phase field fracture models were developed by regularizing the sharp crack topology by a diffusedamage band through the introduction of phase field parameter that discriminate the intact and bro-ken material [18, 19, 20, 21, 22, 23, 24]. Phase field fracture models have been employed in order toinvestigate the fracture behavior of the electrodes in LIBs during lithiation/delithiation process. Zuoand Zhao [7] developed a PFM coupling Li diffusion and stress evolution with crack propagation toinvestigate the behavior of Si thin film electrodes. Miehe et al. [25] studied chemo-mechanical inducedfracture in LIBs by PF fracture modeling. Zuo et al. [26] proposed a PFM coupling Li diffusion, finitedeformation, stress evolution with crack propagation to study spherical Si electrodes in LIBs. Guanet al. [27] studied the stress evolution and crack propagation in the SEI layer formed on the LiMn O electrode during the Li-ion diffusion by using a PFM.In the isotropic PF fracture models, cracking may arise in regions under compression, thus leadingto unphysical crack growth patterns. All the PF fracture models that have been discussed yet areisotropic models. In order to prevent cracking under compressive load states, a tension/compressionsplit of the strain energy is proposed in the so-called anisotropic models [28, 29]. In some PFsimulations of lithiation induced stress and fracture in electrode particles in the literature such as[30, 31, 32, 33, 34, 35, 36], the split of the strain energy into tensile and compressive parts is con-sidered. By using PFM of fracture coupled with anisotropic Cahn-Hilliard-type diffusion, Zhao et al.[30] investigated the electrochemical reaction in LIBs. A chemo-mechanical coupled computationalframework was formulated by Zhang et al. [31] to model diffusion induced large plastic deformationand PF fracture of Si electrodes. Klinsmann et al. developed a coupled model of Li diffusion, me-chanical stress and crack growth, that uses a PF method for the fracture to investigate crack growthin LiMn O spherical and cylindrical particles during Li extraction [32] and insertion[33] in the firsthalf cycle and in the 2nd half cycle [34]. Xu et al. [35] investigated the fracture behavior of cylinderelectrode particles by using a finite strain PF fracture model, in which phase segregation and elec-trochemical reaction on both particle surfaces and fracture surfaces have been taken into account.Nguyen et al. [36] developed a new formulation based on the PF method for modeling stress corrosioncracking induced by anodic dissolution. They coupled classical phase transition model for materialdissolution with the mechanical problem and applied their model to an aluminum alloy in a salinemedium.From the computational viewpoint, anisotropic models are significantly more expensive than the7sotropic ones, since the tension/compression split of the strain energy leads to nonlinear balanceof momentum equations [37]. To overcome this issue, the so-called hybrid PFM has been proposedin recent years [38, 39, 37, 40]. The hybrid formulation formally comprises features from both theisotropic and anisotropic models [37]. In this paper, the author proposes a FE-based hybrid PF modelto study the fracture behavior of electrode particle in LIBs during lithiation. The proposed modelwith coupling effects among Li diffusion, stress evolution and crack propagation, enables a significantreduction of computational cost in comparison with the available anisotropic models.
2. Phase field model
In this section, a PF model with hybrid formulation based on the total free energy of the systemis developed, which is characterized by elastic deformation, Li concentration and crack propagation.
The total free energy of the system is formulated as a functional of the parameters, u , c and φ [26] F = (cid:90) (cid:18) f u + f c + f φ + 12 k c ( ∇ c ) + 12 k φ ( ∇ φ ) (cid:19) d V , (1)where f u , f c and f φ are the local free energy densities of the displacement field u , Li concentrationfield c and fracture order field φ , and ∇ c and ∇ φ are the gradients of concentration and fracturerespectively. Also, k c and k φ represent the gradient energy coefficients.Considering the coupling between the elastic field and the fracture field, f u is described as follows f u = 12 g ( φ ) σ : ε = g ( φ ) ξ u , (2)where ξ u = σ : ε indicates the elastic energy density function, σ denotes the Cauchy stress tensorand ε is the infinitesimal strain tensor which is given by ε = 12 (cid:0) ∇ u + ( ∇ u ) T (cid:1) . (3)The energetic degradation function g ( φ ) couples the elastic field and fracture order parameter byconsidering the stiffness loss between an intact ( φ = 1) and a fully broken region ( φ = 0). The g ( φ )has to be monotonically decreasing and also satisfies g (0) = 0, g (1) = 1, and ∂g (0) /∂φ = 0. Herein,the fourth quartic polynomial function g ( φ ) = 4 φ − φ + η is chosen from the several availabledegradation functions in the literature [37, 41]. The parameter η is a residual stiffness for a fullybroken region and its value may not be too small because of stability issues [42]. The influence of theconcentration field on the elastic field is modeled in analogy to thermal effects [43, 44]. Therefore, thetotal infinitesimal strain in the material is decomposed into chemical and elastic parts as ε = ε e + ε c , (4)8here ε e is the elastic strain and ε c is the chemical strain which is assumed as a hydrostatic dilatationas ε c = (Ω c ) / I , where Ω is the partial molar volume and I is the second-order unit tensor. Theconstitutive equation for the mechanical stresses is given by σ = D : ε e = D : (cid:18) ε − c Ω3 I (cid:19) . (5)Therein, D is the fourth-order stiffness tensor.In general, f c is modeled by a regular solution model [45] or an ideal solution model [46]. Here,for f c it holds [46] ∂f c ∂c = µ + N A k B T c max ln c , (6)where the constant parameter µ , is the standard chemical potential, N A is Avogadro’s constant, k B is the Boltzmann constant, c max is the maximum Li concentration, and T is the absolute temperature.From the [29, 47], f φ is modeled as follows f φ = G cr l (1 − φ ) , (7)where, G cr is the critical energy release rate and l controls the width of the transition area betweenthe unbroken and the broken region. Also, the k φ parameter in Eq. (1), is assumed to be k φ = G cr l .In the present multi-field modeling, the following coupling effects are considered: • Li diffusion affects elastic field by the chemical strain ε c . • Fracture affects elastic field by the energetic degradation function g ( φ ).It should be noted that no direct coupling between Li concentration and fracture is considered in thismodel. Nevertheless, the coupling effects among Li diffusion, elastic field and crack propagation arecaptured via bridging interactions, i.e. Li concentration can cause elastic deformation which can leadto crack propagation. For each variation of (1) the Cahn-Hilliard equation and the Ginzburg-Landau equation can bederived through a variational method [48]. Due to the characteristic time of the elastic field being farless than the concentration field and fracture field, the evolution equation of displacement is assumedto be quasi-static. The governing equations of the system are derived by solving the Ginzburg-Landau equations for the evolution of locally non-conserved parameters u and φ , and the Cahn-Hilliard9quation for the evolution of locally conserved parameter c , [7] ∂ u ∂t = − δFδ u = ∇ ∂ ¯ f∂ ( ∇ u ) − ∂ ¯ f∂ u = 0 , (8) ∂φ∂t = − χ δFδφ = χ (cid:18) ∇ ∂ ¯ f∂ ( ∇ φ ) − ∂ ¯ f∂φ (cid:19) , (9) ∂c∂t = ∇ M ∇ (cid:18) δFδc (cid:19) , (10)where ¯ f is the integral term of the total free energy of the system in (1), and M and χ are molecularmobility and relaxation constant of the fracture order filed respectively.Equation (8) is the balance of linear momentum which can be written explicitly as ∇ · ( g ( φ ) σ ) = 0 . (11)According to (1) and (6), (10) finally yields the mass conservation equation which can be writtenexplicitly as follows ∂c∂t + ∇ · J = 0 , (12)where J is the ion flux vector given by J = − M k B T (cid:20) ∇ c − (cid:18) Ω N A K B T (cid:19) c ∇ ( g ( φ ) σ p ) (cid:21) , (13)with hydrostatic stress σ p , in turn defined by σ p = 13 σ : I . (14)So, the evolution of the Li concentration can be written in a more convenient form by inserting (13)into (12) as [7] ∂c∂t − M k B T ∇ c + M c k c N A c max ∇ c + M k c N A c max ∇∇ c ·∇ c + M c Ω N A ∇ ( g ( φ ) σ p )+ M Ω N A ∇ c ·∇ ( g ( φ ) σ p ) = 0 . (15)Following [49, 7], for convenience in the numerical simulations, the fourth order terms in (15) areneglected in numerical study in Sec. 3.By inserting Eq. (1) and (7) to Eq. (9), the crack propagation equation becomes1 χ ∂φ∂t + G cr l ∇ φ + G cr l (1 − φ ) − g (cid:48) ( φ ) ξ u = 0 . (16)The g (cid:48) ( φ ) is the derivative of energetic degradation function with respect to φ .Equations (11), (15) and (16) are the strong form of governing equations of the electrode forisotropic PF model. 10 .3. Hybrid formulations The isotropic models allow for cracking in regions under compression and interpenetration of thecrack faces, hence it yields physically unrealistic crack evolution patterns. To address this issue, theanisotropic models suggest a tension/compression split of elastic energy density function, i.e., f u = g ( φ ) ξ + u + ξ − u , (17)such that (16) turns into 1 χ ∂φ∂t + G cr l ∇ φ + G cr l (1 − φ ) − g (cid:48) ( φ ) ξ + u = 0 , (18)where the energetic degradation function g ( φ ) is only applied to the tension part of the elastic energydensity function.Although the anisotropic formulation prevent crack growth in compressive regions, it leads tononlinear balance of momentum equations, hence makes the numerical treatment more expensive [37].The general idea of hybrid (isotropic-anisotropic) formulation is to retain a linear momentum balanceequation and also prevent crack growth in regions under compressive load states. To this end, the notseparated form of f u in (2) is kept from isotropic model, while (16) is replaced by (18) from anisotropicmodel which contains tension part of the elastic energy density. Therefore, Equations (11), (15) and(18) are the complete set of governing equations for hybrid PF model. The only seeming disadvantageof the hybrid formulation is that it is variationally inconsistent [37], however, it does not violate thesecond law of thermodynamics [50].The tension part of the elastic energy density function, ξ + u is given by [50] ξ + u = 12 E (cid:104) ¯ σ (cid:105) , (19)where E is Young’s modulus, ¯ σ denotes the largest principal value of the effective stresses which isdescribed in Appendix A, and Macaulay brackets (cid:104) (cid:3) (cid:105) are defined as (cid:104) x (cid:105) = max { x, } . Utilizing the tests functions η c for c , η φ for φ , η u for u , and integrate over the domain, results inthe weak forms of (11), (15) and (18) as (cid:90) V ∇ η u : ( g ( φ ) σ ) d V = 0 , (20) (cid:90) V ˙ c η c d V + (cid:90) V ( M k B T ) ∇ c · ∇ η c d V − (cid:90) V (cid:18) M c Ω N A (cid:19) ∇ ( g ( φ ) σ p ) · ∇ η c d V = 0 , (21) (cid:90) V χ ˙ φ η φ d V + (cid:90) V G cr l ∇ φ · ∇ η φ d V − (cid:90) V G cr l (1 − φ ) η φ d V + (cid:90) V g (cid:48) ( φ ) ξ + u η φ d V = 0 , (22)11here ˙ c and ˙ φ are the material time derivative of the Li concentration field and fracture order parame-ter respectively. It should be noted that, the external nodal forces applied on the electrode boundariesare considered in the global residual vector by means of the FE code.
3. Finite element implementation
In the implementation of numerical method, matrix-vector notation is more convenient than tensornotation. Hence, we switch to the matrix-vector notation for the rest of this paper. In matrix-vectornotation, a second-order symmetric tensor is expressed using a vector, while a fourth-order symmetrictensor is expressed using a matrix.Since the numerical example considered in the present work can be accurately and effectively betreated by a two-dimensional plane strain model. So, the σ, ε and u from (5) and (3) are defined invector form as σ = [ σ xx , σ yy , σ xy ] T , ε = [ ε xx , ε yy , ε xy ] T , u = [ u x , u y ] T . (23) By taking nodal variables as u , c and φ , and utilizing the element shape functions, (20), (21) and(22) are discretized by isoparametric elements u = m (cid:88) I N I u I , ε = m (cid:88) I B I u u I , (24) c = m (cid:88) I N I c I , ∇ c = m (cid:88) I B Ic c I , (25) φ = m (cid:88) I N I φ I , ∇ φ = m (cid:88) I B Iφ φ I , (26)where m is the number of nodes of each element and N I is the nodal values of shape functions.The B I u is the usual strain matrix, and B Ic and B Iφ are the Cartesian derivative matrices, which areexpressed as B I u = ∂N I ∂x ∂N I ∂y∂N I ∂y ∂N I ∂x , B Ic = B Iφ = ∂N I ∂x∂N I ∂y . (27)12ithout losing the generality of formulation, test functions can be expressed as the variations ofdifferent fields, namely η u = m (cid:88) I N I δ u I , ∇ η u = m (cid:88) I B I u δ u I , (28) η c = m (cid:88) I N I δc I , ∇ η c = m (cid:88) I B Ic δc I , (29) η φ = m (cid:88) I N I δφ I , ∇ η φ = m (cid:88) I B Iφ δφ I . (30)So, the variational formulations of (20), (21) and (22) are achieved as (cid:90) V (cid:0) δ u I (cid:1) T (cid:0) B I u (cid:1) T g ( φ ) σ d V = 0 , (31) (cid:90) V N I ˙ c δc I d V + (cid:90) V ( M k B T ) ( B Ic ) T ∇ c δc I d V − (cid:90) V (cid:18) M Ω N A (cid:19) (cid:0) B Ic (cid:1) T ∇ ( g ( φ ) σ p ) c δc I d V = 0 , (32) (cid:90) V χ ˙ φ N I δφ I d V + (cid:90) V G cr l ( B Iφ ) T ∇ φ δφ I d V − (cid:90) V G cr l N I (1 − φ ) δφ I d V + (cid:90) V N I g (cid:48) ( φ ) ξ + u δφ I d V = 0 . (33)Stress vector and hydrostatic stress from (5) and (14) can be decomposed as σ = D ε + D c = D B I u u I + D N I c I , (34) σ p = 13 (cid:104) ( D ) T ε + D c (cid:105) = 13 (cid:104) ( D ) T B I u u I + D N I c I (cid:105) . (35)where D is the 3 × D , D and D for plane-stress assumption are expressed as D = − E Ω3(1 − ν ) [1 , , T , D = E − ν [1 , , T , D = − E Ω3(1 − ν ) , (36)and for plane-strain assumption are expressed as D = − E Ω3(1 − ν ) [1 , , T , D = E − ν [1 , , T , D = − E Ω1 − ν , (37)in which ν is Poisson’s ratio. So, ∇ ( g ( φ ) σ p ) can be calculated as ∇ ( g ( φ ) σ p ) = g ( φ ) ∇ σ p + ∇ g ( φ ) σ p = g ( φ ) (cid:18) (cid:104) ( D ) T ∇ ε + D ∇ c (cid:105)(cid:19) + (cid:0) φ (1 − φ ) ∇ φ (cid:1) σ p . (38)For convenience in differentiation, ( D ) T ∇ ε = D ∗ ∇ ε ∗ , which D ∗ for plane-stress and plane-strainrespectively is D ∗ = E − ν , D ∗ = E − ν , (39)13nd ∇ ε ∗ = u x,xx u y,yx u x,xy u y,yy = n (cid:88) I B I u ∗ u Ii , B I u ∗ = ∂ N I ∂x ∂ N I ∂y∂x∂ N I ∂x∂y ∂ N I ∂y . (40) The field variables vector is defined as Q = { u , c, φ } T , hence based on the quasi-static as-sumption for the displacement, the first and second material time derivative vectors are defined as˙ Q = (cid:110) , ˙ c, ˙ φ (cid:111) T and ¨ Q = (cid:110) , ¨ c, ¨ φ (cid:111) T . An implicit time integration method, a Newmark method withtime step size ∆ t n +1 = t n +1 − t n is used for temporal discretization, where the following approxima-tions for ˙ Q and ¨ Q at time t n +1 are considered [51, 52],˙ Q n +1 = γβ ∆ t n +1 (cid:0) Q n +1 − Q n (cid:1) + (cid:18) − γβ (cid:19) ˙ Q n + (cid:18) − γ β (cid:19) ∆ t n +1 ¨ Q n , (41)¨ Q n +1 = 1 β (∆ t n +1 ) (cid:0) Q n +1 − Q n (cid:1) − β ∆ t n +1 ˙ Q n − − β β ¨ Q n , (42)where γ and β are constants parameters. Considering both accuracy and stability, γ and β are set tobe γ = β = 0 .
5. For more details about possible choices of γ and β , see [52, 53]. The nonlinear set of equations of system, (31), (32) and (33), can be written in the following form R = C ( ˙ Q n +1 ) + P ( Q n +1 ) = 0 . (43)The incremental equation system, Newton method is applied to determine the unknown vector Q n +1 .This leads to the iterative scheme with iteration index i , (cid:20) γβ ∆ t n +1 D + K ( Q n +1 i ) (cid:21) ∆ Q n +1 i +1 = − R ( Q n +1 i ) , (44) Q n +1 i +1 = Q n +1 i + ∆ Q n +1 i +1 , (45)which will be calculated at each time step. The nodal residual vector, R I = (cid:110) R I u , R Ic , R Iφ (cid:111) T , can beexpressed as R I u = (cid:90) V (cid:0) B I u (cid:1) T g ( φ ) σ d V , (46a) R Ic = (cid:90) V ( M k B T ) ( B Ic ) T ∇ c d V − (cid:90) V (cid:18) M Ω N A (cid:19) (cid:0) B Ic (cid:1) T ∇ ( g ( φ ) σ p ) c d V , (46b) R Iφ = (cid:90) V G cr l ( B Iφ ) T ∇ φ d V − (cid:90) V N I (cid:18) G cr l (1 − φ ) − g (cid:48) ( φ ) ξ + u (cid:19) d V . (46c)14he tangential stiffness matrix K and damping matrix D can be derived as K IJ = ∂ P I ∂ Q J , D IJ = ∂ C I ∂ Q J , (47) K IJ = K IJ uu K IJ u c K IJ u φ K IJc u K IJcc K IJcφ K IJφ u K IJφc K IJφφ , D IJ = D IJcc
00 0 D IJφφ . (48)The nodal contributions to the K IJ matrix result from K IJ uu = ∂ P I u ∂ u J = (cid:90) V g ( φ ) (cid:0) B I u (cid:1) T D B J u d V , (49a) K IJ u c = ∂ P I u ∂c J = (cid:90) V g ( φ ) (cid:0) B I u (cid:1) T D N J d V , (49b) K IJ u φ = ∂ P I u ∂φ J = (cid:90) V g (cid:48) ( φ ) (cid:0) B I u (cid:1) T σ N J d V , (49c) K IJc u = ∂P Ic ∂ u J = (cid:90) V − (cid:18) M Ω N A (cid:19) ( B Ic ) T (cid:20) g ( φ ) (cid:18) D ∗ B J u ∗ (cid:19) + (12 φ (1 − φ ) ∇ φ ) (cid:18)
13 ( D ) T B J u (cid:19)(cid:21) c d V , (49d) K IJcc = ∂P Ic ∂c J = (cid:90) V ( M k B T ) ( B Ic ) T B Jc d V − (cid:90) V (cid:18) M Ω N A (cid:19) ( B Ic ) T ∇ ( g ( φ ) σ p ) N J d V − (cid:90) V (cid:18) M Ω N A (cid:19) ( B Ic ) T (cid:20) g ( φ ) (cid:18)
13 D B Jc (cid:19) + (12 φ (1 − φ ) ∇ φ ) (cid:18)
13 D N J (cid:19)(cid:21) c d V , (49e) K IJcφ = ∂P Ic ∂φ J = (cid:90) V (cid:18) M Ω N A (cid:19) ( B Ic ) T (cid:2) g (cid:48) ( φ ) ∇ σ p + σ p (12 φ (1 − φ ) B Jφ + 6 φ (4 + 6 φ ) ∇ φ N J ) (cid:3) c d V , (49f) K IJφ u = ∂P Iφ ∂ u J = (cid:90) V N I (cid:18) ∂ ( g (cid:48) ( φ ) ξ + u ) ∂ u (cid:19) d V , (49g) K IJφc = ∂P Iφ ∂c J = (cid:90) V N I (cid:18) ∂ ( g (cid:48) ( φ ) ξ + u ) ∂c (cid:19) d V , (49h) K IJφφ = ∂P Iφ ∂φ J = (cid:90) V G cr l ( B Iφ ) T B Jφ d V + (cid:90) V N I (cid:18) G cr l + g (cid:48)(cid:48) ( φ ) ξ + u (cid:19) N J d V . (49i)15or calculating ∂ ( g (cid:48) ( φ ) ξ + u ) /∂ u and ∂ ( g (cid:48) ( φ ) ξ + u ) /∂c in 49g and 49h, the readers are referred toAppendix B. The nodal contributions to the D IJ matrix result from D IJcc = ∂C Ic ∂c J = (cid:90) V N I N J d V , (50a) D IJφφ = ∂C Iφ ∂φ J = (cid:90) V χ N I N J d V . (50b)After forming element stiffness, damping and residuals, they are assembled into the global stiffnessmatrix, damping matrix and the right-hand side vector. The resulting set of nonlinear equation havebeen implemented in Matlab. The implementation of the isotropic model is presented in AppendixC.
Crack propagation is an irreversible process. However, the system of equations does not guaran-tee the irreversible evolution of non-conserved PF fracture field, where crack healing is possible. Toprevent crack from healing, irreversibility constraints are considered. Miehe et al. [29, 47] proposeda local history field of the maximum tension part of the elastic energy density function, which auto-matically deals with the irreversibility conditions. Nevertheless, this irreversibility condition cannotbe used in the present PFM since the boundedness cannot be guaranteed for non-quadratic energeticdegradation functions. Herein, a damage like formulation is utilized, where ∆ φ = −(cid:104)− ∆ φ (cid:105) ensuresthe irreversibility constraint.
4. Numerical example
To study the crack propagation under diffusion, the hybrid and also isotropic PF simulation of a SiNW electrode is explored. Two-dimensional circular plate is considered for modeling the cross sectionof NW electrode as shown in Fig. 1. The radius of the circle plate is equal to R = 60 nm and itsthickness is h = 3 µ m [54, 55]. Also, there is an initial crack with a length of a , which is located in thecenter of the specimen plate. The plane strain assumption is conducted and the edges of the models arefree from any displacement constraints. As shown a constant Li concentration, c max = 88 .
67 Kmol/m [8] is applied on boundaries of the electrode. Also, the initial value of concentration in the inner partof electrode is c = c ( t = 0) = 1 . [49]. As the time progresses, Li-ions start to diffuse fromouter boundaries into the inner region of the electrode particle.16 c max c R a x y Figure 1: Geometry of and loading conditions on the NW cross section plate.
The material properties of Si electrode [49, 56, 57, 58] are listed in Tab. 1. Unless otherwise stated,the fracture properties are set to be G cr = 7 N/m and l = 10 nm, and the length of pre-existingcrack is considered to be a = 60 nm. Table 1: Material properties of Si electrode [49, 56, 57, 58].
Parameter Symbol Value UnitYoung’s modulus of Si E Si
80 GPaPoisson’s ratio of Si ν Si E Li–Si
41 GPaPoisson’s ratio of Li–Si ν Li–Si . × − m /molMolecular mobility M
500 m /JsBoltzmann constant k B . × − J/KAbsolute temperature T N A . × χ . × − m /JsCritical energy release rate G cr . − . l −
30 nm
Elastic properties of Li–Si alloys strongly depend on the Li concentration [30, 59]. Herein, basedon a linear rule of mixtures, E = E Li–Si + (1 − c/c max )( E Si − E Li–Si ), the dependency of Young’s17odulus on concentration of lithated Si is considered. Also, with a similar relationship, ν = ν Li–Si +(1 − c/c max )( ν Si − ν Li–Si ), the dependency of Poisson’s ratio on Li concentration is considered.
Choosing a suitable mesh size is critical to obtain sound results at low cost and appropriate timespan. For meshing the model, 4-node isoparametric quadrilateral elements with the refined mesharound the vicinity of the crack are employed. Elaborate mesh convergence analyses are conductedto ensure the full convergence of nonlinear FE solution. Besides, time step convergence analysesshow that for ∆ t ≥ .
006 s the solution process shows divergence and for ∆ t ≤ .
003 s the resultsare virtually converged in time. The ∆ t = 0 . The evolutions of c and g ( φ ) σ p distributions in the electrode during the charging process for theisotropic and hybrid model are depicted in Fig. 2. It can be seen that at the beginning of process, c increases radially from the boundaries into the center of the plate and the interface gradually movesto the center of the NW. As a consequence, volumetric expansion of the electrode is observed whichinduces stress in the electrode. Despite the NW electrode without a pre-existing crack [60], here thevolumetric expansion is not isotropic and the electrode expands more in the direction perpendicularto the crack length. At early stages of lithiation, when higher concentrations are confined to the outerregion of the model, the corresponding volume expansion in this region is hindered by the interiordomain which causes compressive hydrostatic stress in the outer region and tensile hydrostatic stressin the inner region of the NW. It is also observed that hydrostatic stress concentrates at the vicinityof the initial crack tips, then it reduces by time which makes the crack more stable. The aboveoutcomes are in agreement with the previous studies in [11, 26, 10, 14, 15]. In addition, by the arrivalof concentration front to the vicinity of the crack tips, the concentration distribution in electrodebecomes nonuniform and a very high accumulation of c is observed at the tips. This is mainly becausethat Li-ions transfer from lower hydrostatic stress regions towards higher hydrostatic stress regionsbased on (12), which was reported in [7, 27, 13]. Contrary to the c distribution, the distributionof hydrostatic stress becomes uniformed by time evolution. For comprehensive descriptions on thedynamics of lithiation in Si NW one can refer to [60].18 sotropic Hybrid Isotropic Hybrid Figure 2: The evolutions of c and g ( φ ) σ p distributions for isotropic and hybrid models. Fig. 2 does not disclose a remarkable distinction between isotropic model and hybrid model forthe evolutions of c and g ( φ ) σ p . However, as it was shown, the outer regions of the electrode undergocompressive hydrostatic stress. Therefore, it is expected that when the initial crack grows and thecrack tips reach to the regions under compression, hybrid model shows less tendency to crack growththan the isotropic one. Fig. 3 shows the crack propagation during the charging process for isotropicand hybrid models. It is observed that at early stages of the process, both cracks start to propagatewith a similar trend. But as predicted, at the final stages, crack growth in the hybrid model isimpeded, unlike the isotropic model. Also, the grown crack in the isotropic model is evidently thicker,which indicates that the hybrid model has less crack growth tendency toward outer regions whichare under compression. As mentioned before, the presence of a crack at the center of the electrodeplate leads to anisotropic expansion. The length of the electrode in x direction at t = 6 s is 142 . . y direction is 161 . . = 0 s t = 2 s t = 3 s t = 5 s t = 6 s Figure 3: Crack propagation for isotropic and hybrid models.
Fig. 4 compares the crack length and crack thickness of NWs with different initial crack lengthsfor isotropic and hybrid models over the time. Before the start of crack propagation, the crack iscompacted due to anisotropic volumetric expansion which causes a small drop of the crack length.This phenomenon is more pronounced for NW electrodes with larger initial crack. In both models,cracks start to propagate with almost a constant slope until reaching certain values. Then the crackgrowth stops while crack length continues to increase slightly due to the volumetric expansion. Thecrack growth starts sooner and also ends sooner for the NWs with larger initial cracks. At the t = 6 s,the crack lengths are increased 62 . .
79% and 159 .
30% for isotropic models, and 57 . . .
38% for hybrid models, respectively for initial crack lengths of 30 nm, 60 nm and 90 nm.Obviously, the hybrid models show smaller final crack length since they prevent crack growth incompressive regions. The crack lengths for isotropic models at t = 6 s are 13 . .
64% and 4 . Time (s) C r a c k l e n g t h ( n m ) Isotropic, Initial length = 30 nmIsotropic, Initial length = 60 nmIsotropic, Initial length = 90 nmHybrid, Initial length = 30 nmHybrid, Initial length = 60 nmHybrid, Initial length = 90 nm
Time (s) -50510152025303540 C r a c k t h i c k n e ss ( n m ) Isotropic, Initial length = 30 nmIsotropic, Initial length = 60 nmIsotropic, Initial length = 90 nmHybrid, Initial length = 30 nmHybrid, Initial length = 60 nmHybrid, Initial length = 90 nm (a) (b)
Figure 4: (a) Crack length and (b) crack thickness versus time for isotropic and hybrid models.
The crack length during the lithiation process for two different values of G cr are plotted over thetime in Fig. 5 (a). Crack propagation in models with lesser critical energy release rate starts a littlesooner and ends later with significantly bigger crack lengths. Furthermore, the effect of G cr on thecrack length at t = 5 s for isotropic and hybrid models is shown in Fig. 5 (b). The increase of thecritical energy release rate reduces the crack length for both isotopic and hybrid models with a linearpattern, while, as expected, hybrid models show an average of 4 .
1% less crack length.Fig. 6 (a) illustrates the evolution of crack length versus time for different values of l . Similar tothe critical energy release rate, crack propagation for the models with lesser regularization constantstarts sooner and ends later with bigger crack lengths. Also, the increase of l reduces the cracklength for both isotopic and hybrid models with an exponential pattern, as shown in Fig. 6 (b). Theisotropic models show 3 . . .
19% and 6 .
53% larger crack length than the hybrid models for l = 5 nm, l = 10 nm, l = 20 nm and l = 30 nm, respectively.21 Time (s) C r a c k l e n g t h ( n m ) Isotropic, G cr =2.5 N/mIsotropic, G cr =14.9 N/mHybrid, G cr =2.5 N/mHybrid, G cr =14.9 N/m G c (N/m) C r a c k l e n g t h a t t = s ( n m ) Isotropic modelHybrid model (a) (b)
Figure 5: (a) Crack length versus time for different values of G cr and (b) Crack length at t = 5 s versus G c r . Time (s) C r a c k l e n g t h ( n m ) Isotropic, l =5 nmIsotropic, l =30 nmHybrid, l =5 nmIsotropic, l =30 nm l (nm) C r a c k l e n g t h a t t = s ( n m ) Isotropic modelHybrid model (a) (b)
Figure 6: (a) Crack length versus time for different values of l and (b) Crack length at t = 5 s versus l . . Conclusions In this contribution, a hybrid phase field model was proposed to investigate the fracture processof electrodes in Li-ion batteries. The framework links Li-ion diffusion, stress evolution and crackpropagation. The model takes account of the dependency of Elastic properties on the Li concentration.The variational approach along with the finite element method was conducted to obtain the discretizedgoverning equations of the system on space domain. The model was implemented in MATLAB, andall simulations were performed in the same software. Unlike the traditional isotopic models, theproposed hybrid formulation aims to prevent crack growth in the region under compression. Bothisotopic and hybrid models were applied to a silicon nanowire electrode particle, and the effects ofdifferent initial crack lengths and fracture properties were investigated. The numerical results shownthat the hybrid model shows less tendency to crack growth than the isotropic model. Moreover,the increase of the critical energy release rate and regularization constant reduces the crack lengthwith linear and exponential patterns, respectively. The proposed hybrid model has the potential toassist the understanding of the fracture behavior of Li-ion battery electrodes and provides insightfulguidelines for the design of next-generation electrodes in the future. It should be noted that the presentwork is only a preliminary study, and to demonstrate the capability and potential applicability of theproposed hybrid approach there are still different detailed aspects that must be explored. For example,electrodes with geometrical constraints which undergo huge compressive stresses [31, 61], delithiationprocess of an electrode, in which compressive stress is generated in the inner part of the electrode[35, 32, 34, 62], and electrode with multiple pre-existing cracks [7, 13, 35, 14, 25, 32]. Besides, differentelectrode particles with different shapes and sizes must be investigated.
Declarations of interest:
None.
References [1] K. A. Smith, Electrochemical control of lithium-ion batteries, IEEE Control systems magazine30 (2010) 18–25.[2] Y. Zhao, P. Stein, Y. Bai, M. Al-Siraj, Y. Yang, B.-X. Xu, A review on modeling of electro-chemo-mechanics in lithium-ion batteries, Journal of Power Sources 413 (2019) 259–283.[3] A. Mukhopadhyay, B. W. Sheldon, Deformation and stress in electrode materials for li-ionbatteries, Progress in Materials Science 63 (2014) 58–116.[4] S. Kumar, M. Nehra, D. Kedia, N. Dilbaghi, K. Tankeshwar, K.-H. Kim, Carbon nanotubes: Apotential material for energy conversion and storage, Progress in energy and combustion science64 (2018) 219–253. 235] M. Obrovac, L. Christensen, Structural changes in silicon anodes during lithium inser-tion/extraction, Electrochemical and solid-state letters 7 (2004) A93–A96.[6] M. Ebner, F. Marone, M. Stampanoni, V. Wood, Visualization and quantification of electro-chemical and mechanical degradation in li ion batteries, Science 342 (2013) 716–720.[7] P. Zuo, Y.-P. Zhao, A phase field model coupling lithium diffusion and stress evolution withcrack propagation and application in lithium ion batteries, Physical chemistry chemical physics17 (2015) 287–297.[8] I. Ryu, J. W. Choi, Y. Cui, W. D. Nix, Size-dependent fracture of si nanowire battery anodes,Journal of the mechanics and physics of solids 59 (2011) 1717–1730.[9] R. Grantab, V. B. Shenoy, Pressure-gradient dependent diffusion and crack propagation inlithiated silicon nanowires, Journal of the Electrochemical Society 159 (2012) A584–A591.[10] B. Chen, J. Zhou, R. Cai, Analytical model for crack propagation in spherical nano electrodesof lithium-ion batteries, Electrochimica Acta 210 (2016) 7–14.[11] X. Hu, Y. Zhao, R. Cai, J. Zhou, Surface effected fracture behavior of nano-spherical electrodesduring lithiation reaction, Materials Science and Engineering: A 707 (2017) 92–100.[12] H. Wang, E. Oterkus, S. Oterkus, Predicting fracture evolution during lithiation process usingperidynamics, Engineering Fracture Mechanics 192 (2018) 176–191.[13] H. Wang, E. Oterkus, S. Oterkus, Peridynamic modelling of fracture in marine lithium-ionbatteries, Ocean Engineering 151 (2018) 257–267.[14] H. Wang, E. Oterkus, S. Oterkus, Three-dimensional peridynamic model for predicting fractureevolution during the lithiation process, Energies 11 (2018) 1461.[15] Y. Gwak, Y. Jin, M. Cho, Cohesive zone model for crack propagation in crystalline siliconnanowires, Journal of mechanical science and technology 32 (2018) 3755–3763.[16] Y. Hu, X. Zhao, Z. Suo, Averting cracks caused by insertion reaction in lithium–ion batteries,Journal of Materials Research 25 (2010) 1007–1010.[17] I. Steinbach, Phase-field models in materials science, Modelling and simulation in materialsscience and engineering 17 (2009) 073001.[18] B. Bourdin, G. A. Francfort, J.-J. Marigo, Numerical experiments in revisited brittle fracture,Journal of the Mechanics and Physics of Solids 48 (2000) 797–826.2419] B. Bourdin, G. A. Francfort, J.-J. Marigo, The variational approach to fracture, Journal ofelasticity 91 (2008) 5–148.[20] I. Aranson, V. Kalatsky, V. Vinokur, Continuum field description of crack propagation, Physicalreview letters 85 (2000) 118.[21] A. Karma, D. A. Kessler, H. Levine, Phase-field model of mode iii dynamic fracture, PhysicalReview Letters 87 (2001) 045501.[22] V. Hakim, A. Karma, Laws of crack motion and phase-field models of fracture, Journal of theMechanics and Physics of Solids 57 (2009) 342–368.[23] R. Spatschek, E. Brener, A. Karma, Phase field modeling of crack propagation, PhilosophicalMagazine 91 (2011) 75–95.[24] X. Lu, C. Li, Y. Tie, Y. Hou, C. Zhang, Crack propagation simulation in brittle elastic materialsby a phase field method, Theoretical and Applied Mechanics Letters 9 (2019) 339–352.[25] C. Miehe, H. Dal, L.-M. Sch¨anzel, A. Raina, A phase-field model for chemo-mechanical inducedfracture in lithium-ion battery electrode particles, International journal for numerical methodsin engineering 106 (2016) 683–711.[26] P. Zuo, Y.-P. Zhao, Phase field modeling of lithium diffusion, finite deformation, stress evolutionand crack propagation in lithium ion battery, Extreme mechanics letters 9 (2016) 467–479.[27] P. Guan, L. Liu, Y. Gao, Phase-field modeling of solid electrolyte interphase (sei) cracking inlithium batteries, ECS Transactions 85 (2018) 1041–1051.[28] H. Amor, J.-J. Marigo, C. Maurini, Regularized formulation of the variational brittle fracturewith unilateral contact: Numerical experiments, Journal of the Mechanics and Physics of Solids57 (2009) 1209–1229.[29] C. Miehe, F. Welschinger, M. Hofacker, Thermodynamically consistent phase-field models offracture: Variational principles and multi-field fe implementations, International Journal forNumerical Methods in Engineering 83 (2010) 1273–1311.[30] Y. Zhao, B.-X. Xu, P. Stein, D. Gross, Phase-field study of electrochemical reactions at exte-rior and interior interfaces in li-ion battery electrode particles, Computer methods in appliedmechanics and engineering 312 (2016) 428–446.2531] X. Zhang, A. Krischok, C. Linder, A variational framework to model diffusion induced largeplastic deformation and phase field fracture during initial two-phase lithiation of silicon electrodes,Computer methods in applied mechanics and engineering 312 (2016) 51–77.[32] M. Klinsmann, D. Rosato, M. Kamlah, R. M. McMeeking, Modeling crack growth during liextraction in storage particles using a fracture phase field approach, Journal of the electrochemicalsociety 163 (2016) A102–A118.[33] M. Klinsmann, D. Rosato, M. Kamlah, R. M. McMeeking, Modeling crack growth during liinsertion in storage particles using a fracture phase field approach, Journal of the Mechanics andPhysics of Solids 92 (2016) 313–344.[34] M. Klinsmann, D. Rosato, M. Kamlah, R. M. McMeeking, Modeling crack growth during liextraction and insertion within the second half cycle, Journal of power sources 331 (2016) 32–42.[35] B.-X. Xu, Y. Zhao, P. Stein, Phase field modeling of electrochemically induced fracture in li-ionbattery with large deformation and phase segregation, GAMM-Mitteilungen 39 (2016) 92–109.[36] T. T. Nguyen, J. Bolivar, Y. Shi, J. R´ethor´e, A. King, M. Fregonese, J. Adrien, J.-Y. Buffiere,M.-C. Baietto, A phase field method for modeling anodic dissolution induced stress corrosioncrack propagation, Corrosion Science 132 (2018) 146–160.[37] M. Ambati, T. Gerasimov, L. De Lorenzis, A review on phase-field models of brittle fracture anda new fast hybrid formulation, Computational Mechanics 55 (2015) 383–405.[38] J.-Y. Wu, A unified phase-field theory for the mechanics of damage and quasi-brittle failure,Journal of the Mechanics and Physics of Solids 103 (2017) 72–99.[39] J.-Y. Wu, A geometrically regularized gradient-damage model with energetic equivalence, Com-puter Methods in Applied Mechanics and Engineering 328 (2018) 612–637.[40] J.-Y. Wu, V. P. Nguyen, A length scale insensitive phase-field damage model for brittle fracture,Journal of the Mechanics and Physics of Solids 119 (2018) 20–42.[41] C. Kuhn, A. Schl¨uter, R. M¨uller, On degradation functions in phase field fracture models,Computational Materials Science 108 (2015) 374–384.[42] C. Kuhn, R. M¨uller, A continuum phase field model for fracture, Engineering Fracture Mechanics77 (2010) 3625–3634.[43] S. Prussin, Generation and distribution of dislocations by solute diffusion, Journal of appliedphysics 32 (1961) 1876–1881. 2644] J. C.-M. Li, Physical chemistry of some microstructural phenomena, Metallurgical transactionsA 9 (1978) 1353–1380.[45] G. K. Singh, G. Ceder, M. Z. Bazant, Intercalation dynamics in rechargeable battery materials:General theory and phase-transformation waves in lifepo4, Electrochimica acta 53 (2008) 7599–7613.[46] F. Yang, Interaction between diffusion and chemical stresses, Materials science and engineering:A 409 (2005) 153–159.[47] C. Miehe, M. Hofacker, F. Welschinger, A phase field model for rate-independent crack prop-agation: Robust algorithmic implementation based on operator splits, Computer Methods inApplied Mechanics and Engineering 199 (2010) 2765–2778.[48] L.-Q. Chen, Phase-field models for microstructure evolution, Annual review of materials research32 (2002) 113–140.[49] Y. Xie, M. Qiu, X. Gao, D. Guan, C. Yuan, Phase field modeling of silicon nanowire basedlithium ion battery composite electrode, Electrochimica acta 186 (2015) 542–551.[50] J.-Y. Wu, Robust numerical implementation of non-standard phase-field damage models forfailure in solids, Computer Methods in Applied Mechanics and Engineering 340 (2018) 767–797.[51] N. M. Newmark, et al., A method of computation for structural dynamics, American Society ofCivil Engineers, 1959.[52] P. Wriggers, Nonlinear finite element methods, Springer Science & Business Media, 2008.[53] O. C. Zienkiewicz, R. L. Taylor, R. L. Taylor, R. Taylor, The finite element method: solidmechanics, volume 2, Butterworth-heinemann, 2000.[54] T. Toriyama, Y. Tanimoto, S. Sugiyama, Single crystal silicon nano-wire piezoresistors for me-chanical sensors, Journal of microelectromechanical systems 11 (2002) 605–611.[55] Y. Li, P. Gao, Q. Chen, J. Yang, J. Li, D. He, Nanostructured semiconductor solar absorberswith near 100% absorption and related light management picture, Journal of physics D: Appliedphysics 49 (2016) 215104.[56] M. Pharr, Z. Suo, J. J. Vlassak, Measurements of the fracture energy of lithiated silicon electrodesof li-ion batteries, Nano letters 13 (2013) 5570–5577.[57] L. A. Berla, S. W. Lee, Y. Cui, W. D. Nix, Mechanical behavior of electrochemically lithiatedsilicon, Journal of Power Sources 273 (2015) 41–51.2758] V. B. Shenoy, P. Johari, Y. Qi, Elastic softening of amorphous and crystalline li–si phaseswith increasing li concentration: a first-principles study, Journal of Power Sources 195 (2010)6825–6830.[59] J. Ratchford, B. Schuster, B. Crawford, C. Lundgren, J. Allen, J. Wolfenstine, Young’s modulusof polycrystalline li22si5, Journal of Power Sources 196 (2011) 7747–7749.[60] B. Eidel, M. Ahmadi, Insights into the interplay of diffusion and stress in the intercalationdynamics of a li-ion battery electrode by phase-field finite element modeling (2020).[61] S. W. Lee, H.-W. Lee, I. Ryu, W. D. Nix, H. Gao, Y. Cui, Kinetics and fracture resistance oflithiated silicon nanostructure pairs controlled by their mechanical interaction, Nature commu-nications 6 (2015) 7533.[62] X. H. Liu, S. Huang, S. T. Picraux, J. Li, T. Zhu, J. Y. Huang, Reversible nanopore formation inge nanowires during lithiation–delithiation cycling: An in situ transmission electron microscopystudy, Nano letters 11 (2011) 3991–3997. 28 ppendix A. The largest principal value of the effective stresses
The largest principal value of the effective stresses can be obtained by¯ σ = max (cid:8) T T σ , T T σ (cid:9) , (A.1)where T and T are the rotational transformation vectors which are defined as T = cos θ sin θ θ cos θ , T = sin θ cos θ − θ cos θ , (A.2)and the rotation angle θ is obtained as followstan(2 θ ) = 2 σ xy σ xx − σ yy . (A.3) Appendix B. Tension part of the elastic energy density function derivatives
Adopting the chain rule of differentiation, the tension part of the elastic energy density functionderivatives with respect to u and c respectively are derived as ∂ ( g (cid:48) ( φ ) ξ + u ) ∂ u = ∂ ( g (cid:48) ( φ ) ξ + u ) ∂ σ ∂ σ ∂ u = g (cid:48) ( φ ) ∂ ( ξ + u ) ∂ σ D B J u , (B.1) ∂ ( g (cid:48) ( φ ) ξ + u ) ∂c = ∂ ( g (cid:48) ( φ ) ξ + u ) ∂ σ ∂ σ ∂c = g (cid:48) ( φ ) ∂ ( ξ + u ) ∂ σ D N J . (B.2)From Appendix A, ξ + u = 12 E (cid:104) max (cid:8) T T σ , T T σ (cid:9) (cid:105) . (B.3)Then it obtains as ∂ ( ξ + u ) ∂ σ = 1 E (cid:104) max (cid:8) T T σ , T T σ (cid:9) (cid:105) (cid:104) max (cid:8) T T σ , T T σ (cid:9) σ (cid:105) . (B.4) Appendix C. Isotropic model implementation