Fracture Toughness of Li x Si Alloys in Lithium Ion Battery

: Fundamental understanding of the fracture toughness of the Li x Si alloys is crucial for designing of Si based high-capacity and failure-resistant electrodes. In this study, molecular dynamics simulation informed continuum chemo-mechanical modelings with conservation integrals were conducted to derive fracture toughness of Li x Si alloys. Our modeling results show reasonable agreement with available experimental data, revealing that the fracture toughness of Li x Si alloys with low lithium concentration does not vary significantly with lithium concentration. In addition, we demonstrated that, if lithium redistribution caused by the stress gradient around crack tip needs to be considered, an appropriate chemo-mechanical path-independent J-integral should be used as the classic Rices J -integral is path-dependent. The obtained fracture toughness of the Li x Si alloys here


Introduction
2][3][4] However, the relative low capacity of currently commercialized LIBs has seriously hindered their extensive application in some emerging applications, such as electric vehicles and large energy grid. 5- 7 2][13][14][15][16][17][18] Therefore, a fundamental understanding of the degradation mechanisms in the Si anode during lithiation/delithiation cyclings, especially the damage initiation and propagation, is crucial for the rational design of next-generation failure-resistant Si based electrodes.
][16]19 For amorphous Si (a-Si), the volumetric change is ~280%.[20][21] However, the critical fracture size of a-Si particles can increase to ~870 nm.20 In order to investigate the fracture behavior of lithiated Si, Pharr et al. estimated the fracture energy of lithiated Si thin film by measuring the curvature changes of the film at different lithiation and delithiation stages, concluding that the fracture energy of lithiated Si is almost independent of Li concentration as the obtained fracture toughness is fracture with the increase of Li concentration in Si. 23 In consideration of the discrepancy between these two experimental studies, it is worthwhile to further investigate of the fracture toughness of LixSi as it ultimately dictates the cyclability of Si based electrodes.
0][31][32] Therefore, in this study, a MD simulation informed continuum chemomechanical model is presented for the calculation of fracture toughness of LixSi alloys.In particular, MD simulations are conducted for the material properties of LixSi alloys (i.e., Young's modulus, Poisson's ratio, and yield stress) and the critical fracture strain under which the crack in the pre-cracked specimen starts to propagate.Based on the obtained material properties and critical fracture strain, a continuum chemo-mechanical model is used to calculate the stress distribution in the same cracked specimen under the corresponding critical strain loading.With the stress distribution information, J-integral is directly calculated for the fracture toughness of LixSi alloys.

MD simulation
4] This 2NN-MEAM interatomic potential has been adequately demonstrated to provide the primary properties and lattice structures of both crystalline and amorphous LixSi alloy successfully. 33Based on this potential, MD simulations are conducted for the material properties of LixSi alloys and the critical fracture strain under which the crack in the pre-cracked specimen starts to propagate.
The atomic structures of different LixSi alloys are created by melting-quenching method. 35First, large enough pure c-Si block is generated by using LAMMPS MD simulator.Second, according to the ratio x between the number of Li and Si atoms of different LixSi alloys, Li atoms are randomly inserted into the Si block for the corresponding compositions.Then, with periodic boundary conditions (PBCs) prescribed in all three dimensions of the simulation cell and the time step set to be 1 fs, the temperature of the LixSi block is increased from 0 K to 2500 K during the first 100 ps, and maintained there for another 400 ps, and afterward decreased to 10 K with quenching rate ~5.0 K/ps during the following 500 ps.Finally, the obtained fully amorphized LixSi block is relaxed for another 500 ps using NPT ensemble.
To understand the mechanical behaviors of LixSi alloys, MD simulations at 10 K are conducted to obtain the uniaxial tensile stress-strain relationships.Columns (30 × 30 × 80 Å in x, y, and z directions) are cut from the fully relaxed LixSi blocks.Tensile simulations with NPT ensemble are carried out by applying a constant strain rate Based on the obtained stress-strain curves, the mechanical properties of LixSi alloys, such as Young's modulus and yield stress, are estimated.
For the study of the fracture behaviors of LixSi alloys, simulation cells with the size of 500 × 350 × 15 Å in x, y, and z directions are cut from the fully relaxed LixSi blocks.
PBCs are applied in both y and z directions.An edge crack with 0 a 90 = Å initial length and 4 Å width in x-y plane is created in the middle of the cell along x direction, as illustrated in Figure 1(a).After fully relaxed, the cell is stretched in y direction by changing the length of simulation box with a constant strain rate  of the simulation cell at the critical moment that the crack starts to propagate.

Chemo-mechanical simulation
The previously developed finite strain chemo-mechanical model is adopted in this study. 18,32,36 I this model, the concurrent chemical reaction and diffusion processes are simulated in a unified manner by a non-linear diffusion model, despite the difference in interfacial reaction and bulk diffusion.Therefore, the classical diffusion equation is used to describe Li transport in the Si: where, c is the Li concentration, D is the diffusion constant, F J is the Li flux that related to the chemical potential of Li, µ , in Si defined as follows: 37 where, e d , p d , and c d are the elastic, plastic, and chemical stretch rates, respectively.
With the same-size edge-cracked LixSi samples as used in above MD simulation, the following equilibrium equation: is solved with corresponding boundary conditions.All the material properties obtained from the above MD simulations, such as Young's modulus, Poisson's ratio, and yield stress, are adopted as inputs for the modeling.Using the implicit coupled temperaturedisplacement procedure in ABAQUS/Standard with user subroutines, the corresponding Li concentration and stress-strain fields are updated incrementally.Meanwhile, contour integrations are conducted for the classic Rices J -integral: 38 ( ) where, as shown in Figure 2, i t is the traction vector on the contour ∂Ω that around the crack tip, i n is the unit outward normal on ∂Ω , i u is the displacement vector, dΓ is the length element along the contour, and e W is the elastic energy density.

Results and discussion
Upon successfully setting the simulation cell with corresponding boundary and loading conditions, uniaxial tensile simulations of LixSi samples, such as crystalline <111>-Si, a-Si, a-Li0.25Si,a-Li0.50Si,a-Li0.75Si,and a-Li1.00Si,are conducted at 10 K under the strain rate of 7 5 10 1 s × .For the amorphous LixSi samples, eight different configurations are simulated for each compositional case to average the obtained results.As shown in Figure 3, the obtained stress-strain curves of a-LixSi exhibit the typical characteristics of metallic materials. 35Material properties of pure Si and LixSi alloys, such as Young's modulus and yield stress, are estimated from the stress-strain curves.As listed in the Table 1, the Young's modulus and yield stress of LixSi alloys are highly dependent on the concentration of Li, while the Poisson's ratio shows much less sensitivity to the variation of Li concentration.Using the material properties obtained in the above MD simulations as inputs, continuum level chemo-mechanical calculations are conducted for the fracture toughness of pure Si and LixSi alloys.The same-size pre-cracked 2D plane strain samples are used, and the corresponding critical strain is applied as the boundary condition.As no Li redistribution is observed in the above MD simulations, Li distribution in the samples is assumed to be uniform throughout the chemo-mechanical simulation.Therefore, for simplicity of the model, no Li transport needs to be considered in the chemo-mechanical modeling.As a demonstration of the method presented in this study, we first calculated the fracture toughness of <111>-Si as , which is in a good agreement with the reported experimental measurements, namely =0.83~0.95MPa7] This calculation indicates that the approach presented here is appropriate for the derivation of fracture toughness.With confidence in the method, we next calculated the fracture toughness ( C J ) of a-Si and LixSi alloys, as listed in Table 2.The obtained C J for the low Li concentration cases investigated in this study matches the experimental results published by Pharr et al reasonably well. 22In addition, with the increment of Li concentration in the LixSi alloys, the C J decreases first to a minimum value around Li0.25Si and then increases.This trend also agrees with the previous experimental measurement conducted by Wang et al. 23 It should be noted that, as high strain rate is used in the above study, Li transport during the chemo-mechanical calculation of fracture toughness can be neglected.However, in reality, the loading time scale in quasi-static fracture toughness measurements is comparable to the Li diffusion time scale. 43Therefore, it is important to consider the Li redistribution caused by the stress gradient around the crack tip. 45ble 2.The critical strain and fracture toughness ( C J ) of pure Si and LixSi alloys obtained from MD simulation and chemo-mechanical modeling, respectively.In order to overcome the limitation of the classic Rices J -integral when there exists a Li concentration gradient, the newly developed chemo-mechanical J-integral needs to be used for the estimation of fracture toughness. 37,48 he chemo-mechanical J-integral is defined by: where, ψ is the grand potential of the system, T is the absolute temperature, and S is the entropy of mixing per volume of mixture giving by: max max max max max log( ) ( 1) log( 1) Using this chemo-mechanical J-integral, the crack driving forces for the case with Li transport are calculated, as shown in Figure 5.It should be noted that, as the contribution of Li is considered in the grand potential term in Eq. ( 6), the calculated chemo-mechanical J-integral is normalized by

Conclusions
In conclusion, we have presented a molecular dynamics (MD) simulation informed continuum chemo-mechanical model with conservation integrals for the investigation of fracture toughness of LixSi alloys.In this model, MD simulations are used to calculate the material properties and the critical fracture strain of LixSi alloys.Using the obtained MD simulation results as inputs, continuum level chemo-mechanical simulations are conducted for the stress field and thus the fracture toughness of LixSi samples with small scale yielding.
Using the presented model, we investigated the fracture toughness ( C J ) of LixSi alloys with low Li concentration.Our modeling results, which are reasonably consistent with the limited published results derived from different experimental measurements, revealed that the fracture toughness of LixSi alloys with low Li concentration does not vary significantly with Li concentration.In addition, we demonstrated that, when the loading time scale of fracture test is comparable to the Li diffusion time scale, Li redistribution caused by the stress gradient around crack tip needs to be considered, as stress gradient can drive Li to migrate from the far field to crack tip, leading to the aggregation of Li there, and this process is both thermodynamically and kinetically favorable.Under such circumstance, the classic Rices J -integral becomes path-dependent.Therefore, an appropriate chemo-mechanical path-independent J-integral needs to be used for the estimation of fracture toughness.The obtained fracture toughness of the LixSi alloys here provides guidance for the rational design of Si based electrodes, and the approach presented in this study also sheds light for the evaluation of the fracture toughness of other energy materials at different charging/discharging levels.
concentration (~Li2.8Si).22However, the nanoindentation based measurements conducted by Wang et al. revealed a highly Li concentration dependent fracture behavior of lithiated Si, namely, with the increment of Li concentration, the fracture toughness first decreases from 2 2.85 0.15 J m ± for unlithiated a-Si to a minimum value at Li0.31Si, and then increases substantially, 09Si, which indicates a brittle-to-ductile transition of axial direction (z direction) of the columns.In addition, PBCs are prescribed in all three directions of the simulation box and the pressure on the other two non-loading directions is relaxed to zero to ensure the uniaxial tension of the columns.

7 5
10 1 s × .All the MD simulations are conducted at 10 K in LAMMPS with 1 fs time step.The strain at the critical moment that the crack starts to propagate is recorded for the next step chemo-mechanical calculations.

Figure 1 .
Figure 1.A typical LixSi alloy cell for MD simulation with a pre-existed edge crack.(a)

V
ησ are the stress-independent and stress-dependent parts of the chemical potential, respectively.0 µ is a constant representing the chemical potential at a standard state, g R is the gas constant, max c is the concentration corresponding to the saturation state, m V is the molar volume, η is the coefficient of compositional expansion, and kk σ is the trace of stress tensor.For the mechanics part of the model, the finite-strain elasto-plastic framework is used, in which the total stretch rate tensor d is decomposed into three additive parts:

Figure 2 .
Figure 2. A crack tip contour used in the J-integral calculation.The planar crack is oriented

2 C 3 .
60 J m J = based on the contour integration of the classic Rices J -integral, leading to the critical stress intensity factor as initial Li concentration 0 c is loaded with a far field stress 0 σ which is much less than the yield stress Y σ along the vertical direction, the normalized effective stress eff Y / σ σ is shown in Figure 4(a).Without considering Li transport, the initial uniform Li distribution is unaltered in the sample, therefore the normalized changed amount of Li concentration due to the far field loading is 4(b)), and the normalized classic Rices J -integral (Figure 4(c)) obtained from three different integration contours (as shown in Fig. S1 in the Supporting Information) are overlapping with each other, indicating that the J-integral is path-independent.However, if stress gradient driving Li transportation effect is considered and the Li diffusion constant is set to be the same far field stress loading condition, the crack tip stress gradient can drive Li to aggregate around the crack tip (Figure 4(e)), causing more severe stress concentration at the crack tip (Figure 4(d)).In addition, Li concentration gradient around the crack tip can lead to the path-dependent behavior of the classic Rices J -integral, as shown in Figure 4(f).

Figure 4 .
Figure 4. Compared with the normalized effective stress (a), normalized changed amount the classic Rices J - integral is normalized by.With the Li flow caused by crack tip stress gradient, the chemomechanical J-integrals obtained from different integration contours differ slightly from each other.Once the Li redistribution reaches a steady state, the chemo-mechanical Jintegrals will finally converge to a single value, indicating the path-independent property of the chemo-mechanical J-integral.

Figure 5 .
Figure 5. Normalized chemo-mechanical J-integrals obtained from different integration

Table 2 .
It should be noted that, as all the MD simulations are for Li migration.Therefore, no Li redistribution around crack tip is observed, even though previous studies have already demonstrated that Li atoms are highly mobile and can easily