Mechanical Stability of the Heterogenous Bilayer Solid Electrolyte Interphase in the Electrodes of Lithium–Ion Batteries

: Mechanical stability of the solid electrolyte interphase (SEI) is crucial to mitigate the capacity fade of lithium–ion batteries because the rupture of the SEI layer results in further consumption of lithium ions in newly generated SEI layers. The SEI is known as a heterogeneous bilayer and consists of an inner inorganic layer connecting the particle and an outer organic layer facing the electrolyte. The growth of the bilayer SEI over cycles alters the stress generation and failure possibility of both the organic and inorganic layers. To investigate the probability of mechanical failure of the bilayer SEI, we developed the electrochemical-mechanical coupled model with the core–double-shell particle/SEI layer model. The growth of the bilayer SEI is considered over cycles. Our results show that during charging, the stress of the particle changes from tensile to compressive as the thickness of bilayer SEI increases. On the other hand, in the SEI layers, large compressive radial and tensile tangential stress are generated. During discharging, the compressive radial stress of the bilayer SEI transforms into tensile radial stress. The tensile tangential and radial stresses are responsible for the fracture and debonding of the bilayer SEI, respectively. As the thickness ratio of the inorganic to organic layers increases, the fracture probability of the inorganic layer increases, while that of the organic layer decreases. However, the debonding probability of both layers is decreased. In addition, the SEI covering large particles is more vulnerable to fracture, while that covering small particles is more susceptible to debonding. Therefore, tailoring the thickness ratio of the inorganic to organic layers and particle size is important to reduce the fracture and debonding of the heterogeneous bilayer SEI.


Introduction
Lithium-ion batteries are considered a prominent energy storage device, ranging from consumer electronics to electric vehicles [1,2].The solid electrolyte interphase formation (SEI) inside the electrode is crucial to battery performance [3].The developed SEI layer leads to degradation and capacity loss by consuming cyclable lithium ions [4].Current research is focused on the development of a mechanically stable bilayer SEI, as during cycling, the rupture of the SEI layer additionally consumes the cyclable lithium content [5,6].
The decomposition of the electrolyte during initial charge-discharge cycles forms a passivating layer on the electrode surface termed the SEI layer [7].The SEI layer allows lithium ion transport and blocks electrons to prevent further electrolyte decomposition.The SEI layer grows over cycles causing capacity fade [8].Many studies on the bilayer structure of the SEI have been reported.For example, the bilayer SEI develops on the particle surface in the form of a shell, consisting of both inorganic and organic compounds [9,10].Smith et al. reported a heterogeneous bilayer SEI comprised of inorganic (Li 2 CO 3 ) and organic ((CH 2 OCO 2 Li) 2 ) species [11].Peled et al. concluded that the SEI comprises inner inorganic and outer organic layers [9,12].Lan et al. experimentally proved the double-layer hybrid SEI on a SnO 2 electrode [13].Ha et al. prepared a 3D-structured inorganic-organic hybrid bilayer SEI with enhanced mechanical stability [14].Lee et al. observed a bilayer SEI with a 3.7 nm inner inorganic layer and a 15.4 nm outer organic layer [15].Aspern et al. observed that after 1171 charge-discharge cycles, half of the developed SEI layer was inorganic, while the other half was the organic layer [16].Li et al. developed an in situ inner inorganic and outer organic bilayer SEI for zinc aqueous batteries [17].Zhao et al. synthesized an organic-inorganic SEI layer with sufficient mechanical strength [18].
In addition, SEI formation in solid-state batteries has been studied.Fitzhugh et al. developed an ab inito computational method to study the properties of the SEI layer formed at the interface of the solid-state electrolyte [19].Tu et al. investigated the mechanical fracture and debonding of the SEI layer formed on the surface of lithium metal inside solid-state batteries [20].In addition, many efforts have been made to enhance lithium conductivity and reduce the interfacial resistance of the solid-state electrolyte [21][22][23][24].Qin et al. added La 2 O 3 nanoparticles to a garnet-type solid electrolyte to enhance lithium conductivity [25].Using thermodynamic analysis, Qin et al. increased the lithium conductivity of a Ta-doped solid-state electrolyte [26].
The developed SEI layer varies in composition and physiochemical properties, which alters its mechanical stability [27][28][29].Since the continuous conversion of organic SEI to inorganic SEI occurs over several cycles, the actual mechanical properties of the developed SEI cannot be well predicted [30][31][32].In addition, different electrolyte additives also alter the actual thickness and porosity of the developed SEI [33,34].Zheng et al. reported that the SEI layer has an inhomogeneous multilayered structure with varied mechanical properties [35].Moeremans et al. studied the mechanical properties of the heterogeneous bilayer SEI in situ [36].They addressed that the mechanical properties of the developed layers can be controlled via the addition of additives or electrolyte formulations [36].Using the AFM topographical imaging technique, Zhang et al. confirmed that the SEI is inhomogeneous in morphology and mechanical properties [37].Shin et al. investigated that the stiffness of SEI layers widely varies between 0.2 and 80 GPa.Moreover, the inorganic LiF exhibited a peak value of 135.3 GPa [38].
Significant efforts are devoted to developing a mechanically stable bilayer SEI.One example is to develop an elastic SEI [7,39,40].The elastic polymeric films and inherently bonded lithium salts provided considerable mechanical strength [7].In addition, artificial single- [41][42][43] and double-layer SEIs are constructed to enhance mechanical stability, suppressing the fracture and debonding of the SEI layer [44,45].
Various attempts have been made to understand the mechanical failure of the heterogeneous bilayer SEI [46][47][48].For instance, Chen et al. studied the impact of SEI inhomogeneities on the fracture of the SEI layer in Si electrodes [49].Guo et al. examined the cracks in the outer SEI layer.The produced cracks stopped at the inorganic/organic interface [50].He et al. studied the stress in the heterogeneous SEI.They concluded that the peak tensile stress occurred at the active material/inorganic SEI layer interface.They further observed that the strength of SEI layers largely varies with the thickness of the inorganic layer, compared to the organic layer [51].Yuanpeng et al. studied the wrinkling and ratcheting of the SEI layer during lithiation with varying SEI thicknesses [52].
In most of the theoretical models, the bilayer SEI has been assumed to be a continuum with homogenized properties throughout the thickness.This idealized assumption cannot capture the interfacial debonding at the inorganic/organic interfaces [51].In addition, a constant thickness of the bilayer SEI has been assumed over cycles.The constant SEI thickness does not consider the effect of the increased mechanical constraint as the SEI layer grows.Furthermore, the impact of the inorganic/organic thickness ratio of the SEI layer on stress generation and electrochemical performance is rarely studied.In this article, we analyze the mechanical stability of a heterogeneous bilayer SEI, while considering the combined effect of a reduction in the state of charge (SOC) because of lithium consumption and an increase in the mechanical constraint as a result of the growing thickness of the bilayer SEI.The particle-SEI layer is modeled in the form of a core-double-shell, where the inner and outer shells represent the inorganic and organic SEI layers, respectively.A one-dimensional (1D) electrochemical model fully coupled with a two-dimensional (2D) core-double-shell model is developed to evaluate stress generation inside the particle while considering the growth of the bilayer SEI.The stresses are calculated for different inorganicto-organic-layer-thickness ratios.Furthermore, using the elastic strain energies, the fracture and debonding probability of both inorganic and organic SEI layers are analyzed.Finally, the impact of microstructural variation on the mechanical failure probability of the bilayer SEI is studied.

Electrochemical Model
Figure 1 shows the coupling between the 1D electrochemical model (Model 1) and the 2D core-double-shell model (Model 2).Model 1 consists of negative and positive electrodes, a separator, an electrolyte, and current collectors.The active particle and the bilayer SEI are modeled as a core-double-shell, where the active particle is modeled as a core, the inorganic SEI layer (SEI in ) as an inner shell, and the organic SEI layer (SEI or ) as an outer shell.The formation of the bilayer SEI is studied at the negative electrode in Model 1.Using direct projection coupling, the interfacial current density and the SEI current density, which are computed in Model 1, are applied to the particle/inorganic SEI interface (P/SEI in ) in Model 2. During lithiation/delithiation, the electrode region facing the separator experience large stresses.Therefore, the flux calculated from this region in Model 1 is input into Model 2.
constant thickness of the bilayer SEI has been assumed over cycles.The constant SEI thickness does not consider the effect of the increased mechanical constraint as the SEI layer grows.Furthermore, the impact of the inorganic/organic thickness ratio of the SEI layer on stress generation and electrochemical performance is rarely studied.In this article, we analyze the mechanical stability of a heterogeneous bilayer SEI, while considering the combined effect of a reduction in the state of charge (SOC) because of lithium consumption and an increase in the mechanical constraint as a result of the growing thickness of the bilayer SEI.The particle-SEI layer is modeled in the form of a core-double-shell, where the inner and outer shells represent the inorganic and organic SEI layers, respectively.A one-dimensional (1D) electrochemical model fully coupled with a two-dimensional (2D) core-double-shell model is developed to evaluate stress generation inside the particle while considering the growth of the bilayer SEI.The stresses are calculated for different inorganic-to-organic-layer-thickness ratios.Furthermore, using the elastic strain energies, the fracture and debonding probability of both inorganic and organic SEI layers are analyzed.Finally, the impact of microstructural variation on the mechanical failure probability of the bilayer SEI is studied.

Electrochemical Model
Figure 1 shows the coupling between the 1D electrochemical model (Model 1) and the 2D core-double-shell model (Model 2).Model 1 consists of negative and positive electrodes, a separator, an electrolyte, and current collectors.The active particle and the bilayer SEI are modeled as a core-double-shell, where the active particle is modeled as a core, the inorganic SEI layer (SEIin) as an inner shell, and the organic SEI layer (SEIor) as an outer shell.The formation of the bilayer SEI is studied at the negative electrode in Model 1.Using direct projection coupling, the interfacial current density and the SEI current density, which are computed in Model 1, are applied to the particle/inorganic SEI interface (P/SEIin) in Model 2. During lithiation/delithiation, the electrode region facing the separator experience large stresses.Therefore, the flux calculated from this region in Model 1 is input into Model 2. The SEI thickness computed in Model 1 is mapped to the thickness of SEI in and SEI or in Model 2. The thickness of the shells is increased every cycle to consider the growth of the bilayer SEI.Table 1 shows the cell-level parameters.

Cell-Level Model
The time-dependent lithium concentration inside the electrolyte is described by: where D l and ε l are the diffusivity and effective porosity, respectively, t + is the ion transport number, and j is the interfacial current density applied at the P/SEI in interface.The boundary conditions for the electrolyte current density are given in Equation ( 2) and illustrated in Figure 1: ( Inside the electrolyte, the ionic charge balance follows the governing equation [53]: with the following boundary conditions: The lithium diffusion in the particles is written by: where C(r,t) is the lithium concentration and D s represents the lithium ions diffusivity.The mass flux inside the particle is related to interfacial current density as [54]

Interfacial Kinetics
The lithium current density applied at the P/SEI in interface follows the Butler-Volmer equation: where k is the reaction rate constant.c surf and c l are the lithium concentration on the particle surface and electrolyte, respectively.The overpotential η on the P/SEI in interface is calculated as η = ϕ s − ϕ l − E eq − ∆ϕ SEI s , where ϕ s is the electrode potential and ϕ l represents the electrolyte potentials.E eq is the equilibrium potential and ∆ϕ SEI s represents the reduction in potential due to SEI layer resistance.The charge conservation inside the electrode obeys Ohm's law: The boundary conditions for charge balance are written as:

Effect of SEI Formation
During battery charging and discharging, the ethylene carbonate (EC) is reduced to lithium ethylene di-carbonate (CH 2 OCO 2 Li) 2 , consuming the cyclable lithium ions: [55] EC Li] 2 (10) The SEI formation follows the equation: [56] where j 0 SEI is the SEI exchange current density.η SEI = ϕ s − ϕ l − E SEI eq − ∆ϕ SEI s represents the overpotential of the SEI layer.In Model 1, the side reaction is coupled with the electrode reaction to simulate the formation of bilayer SEI and accompanied SOC variation.To solve this, the interfacial kinetics responsible for lithiation and SEI formation reaction are separately defined on the electrode/electrolyte interface [57].The volume fraction of lithium consumed in the bilayer SEI is given as: The temporal evolution of SEI layer volume fraction follows: where M SEI,In is the molecular mass of the SEI in , and ρ SEI,In is the SEI in density.As the volumetric fraction of the bilayer SEI increases, the electrode porosity is reduced.ε l | n is the electrode porosity at the nth cycle.In any cycle, the electrode porosity is calculated as: where ∆ε SEI | n depicts the change in the volume content of the bilayer SEI.The increase in the SEI thickness is: [58] ∂th SEI ∂t = ∂th SEI,In ∂t Over cycles, the thickness of the SEI layer increases as:

Reduction in State of Charge
The battery capacity is defined as: where ε s is the volumetric fraction of the active material inside the positive electrode and c m is the stoichiometric lithium concentration of the material.x max and x min represent the maximum and minimum SOC of the electrode.In this model, no degradation and side reactions are considered at the positive electrode.Part of the available lithium is used during SEI formation, which reduces the discharge capacity, ∆Q SEI .Q| n is the discharge capacity of the nth cycle.Over charge-discharge cycles, the net cell capacity can be obtained as: Substituting Equation ( 16) into Equation ( 17) and solving for x max : The decrease in x max of the positive electrode correspondingly reduces the x max of the anode.

Mechanical Model
During lithium intercalation deintercalation, the stress in Model 2 can be obtained as [59] where ε ij and σ ij represents the strain and stress tensors, respectively, Ω represents the change in volume of the electrode per one addition of lithium.In the current study, the radial σ r and tangential σ θ stress are studied inside Model 2. During charging and discharging, the stress-strain relationship is given as [60,61] du(r, t) In Equations ( 20) and ( 21), the last term on the right-hand side represents the diffusioninduced strain.The mechanical equilibrium is expressed: The governing equation for radial displacement can be obtained by plugging Equations ( 20) and ( 21) into Equation ( 22). where For the particle, the lower limit of integration is α = 0.For the inorganic and organic SEI layers, α = R p , and R p + th SEI,In , respectively, where R p is the particle radius.When there is no SEI layer, the boundary conditions become: [63] The boundary conditions during the SEI layer formation are: are the radial displacements and corresponding stresses of the particle and the bilayer SEI, respectively.The radial displacement u is plugged into Equations ( 20) and ( 21) to obtain the corresponding stresses inside the SEI shell and the core.

Stress inside the Bilayer SEI
Applying the appropriate boundary conditions shown in Equation ( 26), the σ r and σ θ in the SEI In (R p ≤ r ≤ R p + th SEI,In ) are calculated as: Similarly, the stresses inside the organic SEI layer (R p + th SEI,In ≤ r ≤ R t ), are: where In the presence of the SEI layer constraint, the σ r and σ θ inside the particle are calculated as: where m is: where ), and b = th SEI R t .When no SEI layer is formed, m =1 and Equations ( 31) and (32) transform into general analytical equations to calculate the σ r and σ θ [64].

Numerical Simulations
The models are numerically solved using COMSOL software.The Transport of diluted species (tds) module is adopted to solve the lithium concentration of the core, whereas the Solid mechanics module is used to compute the stresses in the core and bilayer SEI.After the mesh independence tests are conducted, the final mesh contained 257,830 triangular elements with 1,985,032 degrees of freedom.During the simulation, the charge-discharge process of the battery is controlled by varying the applied current.The battery is constantly charged at 1 C (cc_ch), until the maximum voltage of 4.25 V is achieved.The battery is further charged at the maximum voltage (cv_ch), while decreasing the current until it decreases to 0.05 C.Then, the battery is discharged at 1 C to 2.7 V (cc_dch).

Capacity Fading
Figure 2A depicts the discharge profiles for the first 2000 cycles.Over cycles, the capacity decreases to 80.2% of the initial capacity.This decrease in capacity is solely due to the formation of the bilayer SEI, which decreases the available cyclable lithium content.Figure 2B shows the reduction in SOC of the anode from the initial (33 to 27) % in 2000 cycles.The consumed ions are used in the formation of the bilayer SEI on the particle surface.Figure 2C shows the evolution of the heterogeneous bilayer SEI thickness.In 2000 cycles, the SEI in and SEI or increase to 0.308 and 0.462 µm, respectively.Figure 2D shows the thickness of SEI in and SEI or with different h SEI .The total thickness of the SEI layer is SEI total = SEI in + SEI or = h SEI SEI total + (1 − h SEI )SEI total .The h SEI is the ratio of SEI in to SEI total .

Stress Analysis
The developed bilayer SEI increases the mechanical constraint to the swelling and shrinkage of the particle during lithiation and delithiation.Over cycles, the bilayer SEI continuously grows on the particle surface, forming a core-double-shell structure.Figure 3a shows the radial stress contours inside Model 2 at 2000 cycles during lithiation.The zoomed view shows that the largest compressive radial stress occurs at the P/SEIin interface and decreases toward the center of the particle center and the surface of SEIor. Figure 3b depicts the tangential stress contour inside Model 2 at 2000 cycles during lithiation.The extended view shows the stress discontinuity at the interfaces.The maximum tangential stress occurred at the P/SEIin interface.Figure 3c shows the detailed contour maps of the tangential stress in the separate domain of the particle, inorganic, and organic SEI lay-

Stress Analysis
The developed bilayer SEI increases the mechanical constraint to the swelling and shrinkage of the particle during lithiation and delithiation.Over cycles, the bilayer SEI continuously grows on the particle surface, forming a core-double-shell structure.Figure 3a shows the radial stress contours inside Model 2 at 2000 cycles during lithiation.The zoomed view shows that the largest compressive radial stress occurs at the P/SEI in interface and decreases toward the center of the particle center and the surface of SEI or .Figure 3b depicts the tangential stress contour inside Model 2 at 2000 cycles during lithiation.The extended view shows the stress discontinuity at the interfaces.The maximum tangential stress occurred at the P/SEI in interface.Figure 3c shows the detailed contour maps of the tangential stress in the separate domain of the particle, inorganic, and organic SEI layers.The largest compressive stress inside the particle occurs at the surface as shown in Figure 3c(i).Figure 3c(ii) depicts that the highest tensile tangential stress in SEI in occurs at the P/SEI in interface and decreases along the thickness.Figure 3c(iii) illustrates that the tensile stress in SEI or is smaller than that in SEI in .Inside the organic layer, the largest stress happens at the SEI in /SEI or interface and decreases along the thickness.The compressive stress of the particle and bilayer SEI is caused by the simultaneous effect of diffusioninduced stress and mechanical confinement from the bilayer SEI.  Figure 4A shows the radial stress distribution inside Model 2 for hSEI = 0.4 at the end of lithiation.At the initial cycles, tensile radial stress occurred inside the particle.As the bilayer SEI becomes thicker during charge-discharge cycles, the corresponding mechanical constraint against the particle expansion increases, causing a reduction in the tensile behavior.After 800 cycles, the SEI layer confinement is sufficient to transform the tensile radial stress to compressive.The compressive stress of the particle increases outward from the particle center, reaches the maximum at the P/SEIin interface, and then decreases to zero along the SEI thickness.Figure 4A shows the radial stress distribution inside Model 2 for h SEI = 0.4 at the end of lithiation.At the initial cycles, tensile radial stress occurred inside the particle.As the bilayer SEI becomes thicker during charge-discharge cycles, the corresponding mechanical constraint against the particle expansion increases, causing a reduction in the tensile behavior.After 800 cycles, the SEI layer confinement is sufficient to transform the tensile radial stress to compressive.The compressive stress of the particle increases outward from the particle center, reaches the maximum at the P/SEI in interface, and then decreases to zero along the SEI thickness.Figure 4B depicts the magnified view of the compressive radial stress of the bilayer SEI, showing that the compressive stress is decreased in a linear manner.The gradient of the stress reduction is decreased at the interface of SEIin and SEIor.
Figure 4C illustrates the tangential stress inside Model 2. Compared to the particle, significantly large tensile tangential stress occurred in the bilayer SEI.The maximum tensile tangential stress occurs in SEIin and decreases along the SEI thickness.Inside the bilayer SEI, the tangential stress shows discontinuities twice at the P/SEIin and SEIin/SEIor interfaces.However, the tangential stress in the SEIor is smaller than that in SEIin.

Effect of hSEI
In this section, we discuss how the variation in hSEI affects the stress distribution in Model 2. Figure 5A depicts the radial stress distribution inside Model 2 with different hSEI.The results are plotted for 2000 cycles.As hSEI increases, the tensile radial stress tends to become compressive stress because thick SEIin provides a larger mechanical constraint to the particle expansion than thin SEIin.Figure 5B shows the magnified view of the stress inside the bilayer SEI. Figure 5C illustrates the tangential stress inside Model 2 with different hSEI.Similar to the radial stress, an increase in hSEI transforms the tangential stress in the core from tensile to compressive.However, the tensile tangential stress of the SEIin  Figure 4C illustrates the tangential stress inside Model 2. Compared to the particle, significantly large tensile tangential stress occurred in the bilayer SEI.The maximum tensile tangential stress occurs in SEI in and decreases along the SEI thickness.Inside the bilayer SEI, the tangential stress shows discontinuities twice at the P/SEI in and SEI in /SEI or interfaces.However, the tangential stress in the SEI or is smaller than that in SEI in .

Effect of h SEI
In this section, we discuss how the variation in h SEI affects the stress distribution in Model 2. Figure 5A depicts the radial stress distribution inside Model 2 with different h SEI .The results are plotted for 2000 cycles.As h SEI increases, the tensile radial stress tends to become compressive stress because thick SEI in provides a larger mechanical constraint to the particle expansion than thin SEI in .Figure 5B shows the magnified view of the stress inside the bilayer SEI. Figure 5C illustrates the tangential stress inside Model 2 with different h SEI .Similar to the radial stress, an increase in h SEI transforms the tangential stress in the core from tensile to compressive.However, the tensile tangential stress of the SEI in is reduced, as shown in Figure 5D.These stresses cause the fracture and debonding of the bilayer SEI.

Fracture and Debonding Analysis
The tensile tangential stress that occurred in the bilayer SEI may cause a fracture of the bilayer SEI.In the core-double-shell structure, the elastic strain energy responsible for the cracking of the bilayer SEI can be calculated as [65]   where x = 'in' or 'or' represents the inorganic or organic portion of the bilayer SEI.The maximum tensile tangential stresses at the P/SEIin and SEIin/SEIor interfaces are chosen for the calculation of elastic strain energy.Figure 6A,B show the evolution of tensile tangential stress inside SEIin and SEIor, respectively, over 2000 cycles.In the initial cycles, the stresses are higher and decrease as the cycling continues.In addition, as hSEI decreases, the tangential stress increases.Figure 6C,D show the corresponding elastic strain energies responsible for SEI rupture.For small hSEI, the fracture possibility of the inorganic layer is

Fracture and Debonding Analysis
The tensile tangential stress that occurred in the bilayer SEI may cause a fracture of the bilayer SEI.In the core-double-shell structure, the elastic strain energy responsible for the cracking of the bilayer SEI can be calculated as [65] where x = 'in' or 'or' represents the inorganic or organic portion of the bilayer SEI.The maximum tensile tangential stresses at the P/SEI in and SEI in /SEI or interfaces are chosen for the calculation of elastic strain energy.Figure 6A,B show the evolution of tensile tangential stress inside SEI in and SEI or , respectively, over 2000 cycles.In the initial cycles, the stresses are higher and decrease as the cycling continues.In addition, as h SEI decreases, the tangential stress increases.Figure 6C,D show the corresponding elastic strain energies responsible for SEI rupture.For small h SEI , the fracture possibility of the inorganic layer is low in all cycles, although the stress is large, as shown in Figure 6A.This is because the increase in thickness has a stronger effect than the stress in Equation (34).On the other hand, the cracking possibility of the organic layer is high for small h SEI after 800 cycles, as shown in Figure 6D.As h SEI increases, the fracture probability of the SEI in increases, while that of the SEI or decreases.Therefore, for small h SEI , a fracture in SEI is more likely to occur at the organic layer, while for large h SEI , it occurs at the inorganic layer.The simulation results suggest that having a similar thickness of the inorganic and organic layers is preferred to avoid a fracture of the SEI layers.
Mathematics 2023, 11, x FOR PEER REVIEW 13 of 20 low in all cycles, although the stress is large, as shown in Figure 6A.This is because the increase in thickness has a stronger effect than the stress in Equation (34).On the other hand, the cracking possibility of the organic layer is high for small hSEI after 800 cycles, as shown in Figure 6D.As hSEI increases, the fracture probability of the SEIin increases, while that of the SEIor decreases.Therefore, for small hSEI, a fracture in SEI is more likely to occur at the organic layer, while for large hSEI, it occurs at the inorganic layer.The simulation results suggest that having a similar thickness of the inorganic and organic layers is preferred to avoid a fracture of the SEI layers.In addition to the fracture, debonding of the bilayer SEI can occur when tensile radial stress develops at the interface.The tensile radial stress develops during delithiation as the shrinkage of the particle is prevented by the surrounding bilayer SEI.The debonding energy release rate for the inorganic SEI layer is given as: [61] Similarly, for the organic layer, Equation ( 35) becomes: In addition to the fracture, debonding of the bilayer SEI can occur when tensile radial stress develops at the interface.The tensile radial stress develops during delithiation as the shrinkage of the particle is prevented by the surrounding bilayer SEI.The debonding energy release rate for the inorganic SEI layer is given as: [61] Similarly, for the organic layer, Equation ( 35) becomes: Figure 7A shows the radial stress distribution along the particle and bilayer SEI at the 500th cycle with h SEI = 0.4.Figure 7B shows the zoomed view near the bilayer SEI.During the charging process (cc_ch), the radial stresses at the P/SEI in and SEI in /SEI or interfaces are compressive, as the swelling of the particle is confined by the bilayer SEI.After the constant current charging (cc_ch) is changed to the constant voltage charging (cv_ch), the compressive radial stress at the interfaces further increases.However, in subsequent discharging (cc_dch), the radial stress at the interfaces transforms from a compressive to a tensile state, because the particle shrinks upon delithiation.This tensile radial stress inside the bilayer SEI increases over cycles as shown in Figure 7C. Figure 7D depicts the debonding energy release rate for whole cycles.Over cycles, the debonding probability of both SEI in and SEI or increases.Figure 7A shows the radial stress distribution along the particle and bilayer SEI at the 500th cycle with hSEI = 0.4.Figure 7B shows the zoomed view near the bilayer SEI.During the charging process (cc_ch), the radial stresses at the P/SEIin and SEIin/SEIor interfaces are compressive, as the swelling of the particle is confined by the bilayer SEI.After the constant current charging (cc_ch) is changed to the constant voltage charging (cv_ch), the compressive radial stress at the interfaces further increases.However, in subsequent discharging (cc_dch), the radial stress at the interfaces transforms from a compressive to a tensile state, because the particle shrinks upon delithiation.This tensile radial stress inside the bilayer SEI increases over cycles as shown in Figure 7C. Figure 7D depicts the debonding energy release rate for whole cycles.Over cycles, the debonding probability of both SEIin and SEIor increases.Figure 8A,B show the contour plot of radial stress for the whole range of cycles and h SEI at the P/SEI in and SEI in /SEI or interfaces, respectively.The maximum stress occurred at later cycles and small h SEI .As the h SEI increases, the stress decreases.Moreover, the radial stress at the P/SEI in interface is compressive in low cycles and at large h SEI .As the cycling proceeds, the radial stress transforms from the compressive to the tensile state.The stress of the SEI in /SEI or interface is greater than that of the P/SEI in interface.
Mathematics 2023, 11, x FOR PEER REVIEW 15 of 20 Figure 8A,B show the contour plot of radial stress for the whole range of cycles and hSEI at the P/SEIin and SEIin/SEIor interfaces, respectively.The maximum stress occurred at later cycles and small hSEI.As the hSEI increases, the stress decreases.Moreover, the radial stress at the P/SEIin interface is compressive in low cycles and at large hSEI.As the cycling proceeds, the radial stress transforms from the compressive to the tensile state.The stress of the SEIin/SEIor interface is greater than that of the P/SEIin interface.Figure 8C,D show the debonding strain energy of the P/SEIin and SEIin/SEIor interface, respectively.The contour map indicates that debonding is more likely to occur as the cycle increases and as hSEI decreases.The white area in Figure 8C depicts the no-debonding region, where the stress is compressive.In our simulation, debonding of the P/SEIin interface does not occur up to 800 cycles for hSEI = 0.8.

Effect of Particle Size on Fracture and Debonding
A series of simulations are carried out to study the effect of particle size on the fracture and debonding of the bilayer SEI.The particle size is varied from 5 to 25 um. Figure 9 shows the fracture and debonding energy release rate of the inorganic and organic SEI layers.As the particle size increases, the fracture is more dominant, while debonding is more likely to occur with the decrease in particle size.Since the fracture energy release Figure 8C,D show the debonding strain energy of the P/SEI in and SEI in /SEI or interface, respectively.The contour map indicates that debonding is more likely to occur as the cycle increases and as h SEI decreases.The white area in Figure 8C depicts the no-debonding region, where the stress is compressive.In our simulation, debonding of the P/SEI in interface does not occur up to 800 cycles for h SEI = 0.8.

Effect of Particle Size on Fracture and Debonding
A series of simulations are carried out to study the effect of particle size on the fracture and debonding of the bilayer SEI.The particle size is varied from 5 to 25 um. Figure 9 shows the fracture and debonding energy release rate of the inorganic and organic SEI layers.As the particle size increases, the fracture is more dominant, while debonding is more likely to occur with the decrease in particle size.Since the fracture energy release rate of the SEI in is higher than that of the SEI or , fracture highly occurs at the inner SEI layer (inorganic layer) of large particles.On the other hand, the debonding energy release rate of the SEI or is greater than that of the SEI in , and debonding highly occurs at the SEI in /SEI or interface of small particles.Our simulation results suggest that tuning the particle size to approximately 13 µm is a better choice to minimize both SEI fracture and debonding.
Mathematics 2023, 11, x FOR PEER REVIEW 16 of 20 rate of the SEIin is higher than that of the SEIor, fracture highly occurs at the inner SEI layer (inorganic layer) of large particles.On the other hand, the debonding energy release rate of the SEIor is greater than that of the SEIin, and debonding highly occurs at the SEIin/SEIor interface of small particles.Our simulation results suggest that tuning the particle size to approximately μm is a better choice to minimize both SEI fracture and debonding.

Conclusions
In this paper, we developed a 1D electrochemical model fully coupled with a coredouble-shell particle/SEI model to investigate the mechanical stability of the heterogeneous bilayer SEI over multiple cycles.The SEI layer is considered a double-layer shell consisting of the inorganic layer as the inner shell and the organic layer as the outer shell.Our simulation results show that the increase in the mechanical constraint due to the growth of the bilayer SEI transforms the tensile stress inside the particle into compressive stress.Tensile tangential stress occurs at the particle/SEI interfaces, which leads to the initiation of a fracture inside the inorganic SEI layer, and the total fracture of both inorganic and organic SEI layers.As the thickness ratio of the inorganic layer increases, the compressive radial stress at the interface increases, while tangential stress decreases.At the end of discharge, the compressive radial stress at the interface converts to tensile stress, which leads to debonding of the interfaces.
As the thickness ratio of the inorganic layer increases, the fracture probability of the inorganic layer increases, while that of the organic layer decreases.On the other hand, the debonding probability of both inorganic and organic layers increases as the thickness ratio of the inorganic layer decreased.In addition, for the effect of particle size, the simulation results suggest that, in a multiparticle electrode, fracture of the SEI layers is more likely to occur for large particles, while the debonding of SEI layers is more likely to occur for small particles.Therefore, tailoring the thickness ratio of the inorganic layer and particle size is important to reduce the fracture and debonding of the heterogeneous bilayer SEI.

Conclusions
In this paper, we developed a 1D electrochemical model fully coupled with a coredouble-shell particle/SEI model to investigate the mechanical stability of the heterogeneous bilayer SEI over multiple cycles.The SEI layer is considered a double-layer shell consisting of the inorganic layer as the inner shell and the organic layer as the outer shell.Our simulation results show that the increase in the mechanical constraint due to the growth of the bilayer SEI transforms the tensile stress inside the particle into compressive stress.Tensile tangential stress occurs at the particle/SEI interfaces, which leads to the initiation of a fracture inside the inorganic SEI layer, and the total fracture of both inorganic and organic SEI layers.As the thickness ratio of the inorganic layer increases, the compressive radial stress at the interface increases, while tangential stress decreases.At the end of discharge, the compressive radial stress at the interface converts to tensile stress, which leads to debonding of the interfaces.
As the thickness ratio of the inorganic layer increases, the fracture probability of the inorganic layer increases, while that of the organic layer decreases.On the other hand, the debonding probability of both inorganic and organic layers increases as the thickness ratio of the inorganic layer decreased.In addition, for the effect of particle size, the simulation results suggest that, in a multiparticle electrode, fracture of the SEI layers is more likely to occur for large particles, while the debonding of SEI layers is more likely to occur for small particles.Therefore, tailoring the thickness ratio of the inorganic layer and particle size is important to reduce the fracture and debonding of the heterogeneous bilayer SEI.

Figure 1 .
Figure 1.Simulation model showing the coupling between the 1D electrochemical model and the 2D core-double-shell model.

r
u SEI,In r r=R p +th SEI,In = u SEI,Or r r=R p +th SEI,In σ p r r=R p = σ SEI,In r r=R p , σ SEI,In r r=R p +th SEI,In = σ SEI,Or r r=R p +th SEI,In σ , u SEI r , and σ p r , σ SEI r

Figure 2 .
Figure 2. (A) Discharge curves over cycles.(B) Decrease in SOC of the negative electrode due to SEI formation.(C) Evolution of SEI layer thickness.(D) Thickness of the inorganic and organic SEI layer according to h SEI .

Figure 2 .
Figure 2. (A) Discharge curves over cycles.(B) Decrease in SOC of the negative electrode due to SEI formation.(C) Evolution of SEI layer thickness.(D) Thickness of the inorganic and organic SEI layer according to h SEI .

Figure 3 .
Figure 3. Stress contours in the particle and SEI layer at the end of 2000 cycles.(A) Radial stress with magnified view in the inorganic and organic SEI layers.(B) Tangential stress with magnified view in the inorganic and organic SEI layers.(C) Detailed tangential stress distribution (i) inside the particle, (ii) inorganic SEI layer, and (iii) organic SEI layer.The results are plotted for hSEI = 0.4.

Figure 3 .
Figure 3. Stress contours in the particle and SEI layer at the end of 2000 cycles.(A) Radial stress with magnified view in the inorganic and organic SEI layers.(B) Tangential stress with magnified view in the inorganic and organic SEI layers.(C) Detailed tangential stress distribution (i) inside the particle, (ii) inorganic SEI layer, and (iii) organic SEI layer.The results are plotted for h SEI = 0.4.

Mathematics 2023 , 20 Figure 4 .
Figure 4. Evolution of radial and tangential stress inside the particle and SEI layer during lithiation.(A) Radial stress inside the particle and SEI layer, with (B) magnified view inside the inorganic and organic SEI layer.(C) Tangential stress inside the particle and SEI layer, with (D) magnified view inside the inorganic and organic SEI layer.

Figure 4 .
Figure 4. Evolution of radial and tangential stress inside the particle and SEI layer during lithiation.(A) Radial stress inside the particle and SEI layer, with (B) magnified view inside the inorganic and SEI layer.(C) Tangential stress inside the particle and SEI layer, with (D) magnified view inside the inorganic and organic SEI layer.

Figure
Figure 4B depicts the magnified view of the compressive radial stress of the bilayer SEI, showing that the compressive stress is decreased in a linear manner.The gradient of the stress reduction is decreased at the interface of SEI in and SEI or .Figure4Cillustrates the tangential stress inside Model 2. Compared to the particle, significantly large tensile tangential stress occurred in the bilayer SEI.The maximum tensile tangential stress occurs in SEI in and decreases along the SEI thickness.Inside the bilayer SEI, the tangential stress shows discontinuities twice at the P/SEI in and SEI in /SEI or interfaces.However, the tangential stress in the SEI or is smaller than that in SEI in .

Mathematics 2023 ,
11, x FOR PEER REVIEW 12 of 20 is reduced, as shown in Figure 5D.These stresses the fracture and debonding of the bilayer SEI.

Figure 5 .
Figure 5.Effect of hSEI on radial and tangential stress inside the particle and SEI layer during lithiation.The results are plotted for 2000 cycles.(A) Radial stress inside the particle and SEI layer, with (B) magnified view inside the inorganic and organic SEI layer.(C) Tangential stress inside the particle and SEI layer, with (D) magnified view of inside the inorganic and organic layers.

Figure 5 .
Figure 5.Effect of h SEI on radial and tangential stress inside the particle and SEI layer during lithiation.The results are plotted for 2000 cycles.(A) Radial stress inside the particle and SEI layer, with (B) magnified view inside the inorganic and organic SEI layer.(C) Tangential stress inside the particle and SEI layer, with (D) magnified view of inside the inorganic and organic layers.

Figure 6 .
Figure 6.Evolution of tangential stress at the (A) P/SEIin interface and (B) SEIin/SEIor interface, with different hSEI.Contour plot of fracture energy release rate Gf as a function of hSEI and cycle numbers for the (C) inorganic and (D) organic SEI layer.

Figure 6 .
Figure 6.Evolution of tangential stress at the (A) P/SEI in interface and (B) SEI in /SEI or interface, with different h SEI .Contour plot of fracture energy release rate G f as a function of h SEI and cycle numbers for the (C) inorganic and (D) organic SEI layer.

Figure 7 .
Figure 7. (A) Radial stress profile along the radius of particle at the end of constant current charge (cc_ch), constant voltage charge (cv_ch), and discharge (cc_dch).(B) Magnified view of the particle/SEI interface and SEI layers.The magnified area is shown in (A).(C) Evolution of radial stress in the inorganic and organic SEI layer interface (hSEI = 0.4).(D) Evolution of debonding energy release rate in the SEI layers (hSEI = 0.4).

Figure 7 .
Figure 7. (A) Radial stress profile along the radius of particle at the end of constant current charge (cc_ch), constant voltage charge (cv_ch), and discharge (cc_dch).(B) Magnified view of the particle/SEI interface and SEI layers.The magnified area is shown in (A).(C) Evolution of radial stress in the inorganic and organic SEI layer interface (h SEI = 0.4).(D) Evolution of debonding energy release rate in the SEI layers (h SEI = 0.4).

Figure 8 .
Figure 8. Contour plots of radial stress and debonding energy release rate Gd as a function of cycle number and hSEI.Radial stress inside the (A) inorganic and (B) organic SEI layers.Debonding energy release rate Gd inside the (C) inorganic and (D) organic SEI layers.In (C), the white area indicates no debonding onset, due to the compressive stresses.

Figure 8 .
Figure 8. Contour plots of radial stress and debonding energy release rate G d as a function of cycle number and h SEI .Radial stress inside the (A) inorganic and (B) organic SEI layers.Debonding energy release rate G d inside the (C) inorganic and (D) organic SEI layers.In (C), the white area indicates no debonding onset, due to the compressive stresses.

Figure 9 .
Figure 9. Fracture and debonding energy release rate of the inorganic and organic SEI layers as a function of particle size.

Figure 9 .
Figure 9. Fracture and debonding energy release rate of the inorganic and organic SEI layers as a function of particle size.

Table 1 .
Parameters for the Cell level.

Table 1 .
Parameters for the Cell level.