Analytical computation of stress intensity factor for multi-physics problems

This work presents a methodology for the analytical calculation of the stress intensity factor when the stress distribution on the crack surfaces is non-homogeneous. At first, a polynomial function is used to express the non-homogenous stress distribution. Subsequently, the principle of superposition of effects is applied, and the stress intensity factor is computed by multiplying each polynomial term by its respective geometric factor. Finite element fracture model is used to compute the geometric factor of the single polynomial grade. To explain the method, a spherical body is considered, with central and superficial cracks. Each geometric factor depends on a normalized geometrical parameter (the ratio between the crack length and sphere radius). The proposed methodology is applied to determine the stress intensity factor in the case of a crack driving force caused by diffusive fields, such as the concentration gradient in particles of electrodes active material in lithium-ion batteries. The methodology allows to speed up the fracture computation, then it is used to give electrode design guidelines to limit the fracture likeliness and mechanical degradation in lithium-ion batteries, as well as it is the basis for the development of algorithms assessing the capacity loss and the remaining useful life of lithium-ion batteries in real-time.


Introduction
Fracture is one of the most common failure modes in several engineering applications.Indeed, cracks propagation and fatigue phenomena occur well below the yield strength of the material and may lead to the failure of the component.It is known that cracks initiate and propagate from geometrical discontinuities in structures under loading, such as micro-defects, edges, holes, and notches, acting as strain and stress raisers [1,2].However, the crack propagation size and the propagation rate depend on the geometry of the discontinuity and on the loading condition.
Different parameters can be used to assess the level of fracture, such as the energy release rate G, the J-integral, or the stress intensity factor (SIF) [3,4].Among these, the SIF is a wellknown parameter assessing the singularity of the stress field near the crack region, according to linear elastic fracture mechanics (LEFM) [5].Fracture criteria describing crack propagation and material failure are usually expressed as a function of SIF both for brittle [3] and fatigue fracture [6].For this reason, the accurate evaluation of SIF is of great importance.
Several analytical expressions for the SIF are available in the literature [7,8,9,10].However, finite element method (FEM) is often preferred, despite the ease of analytical computation, because analytical solutions are only applicable to a limited number of standard cases.In recent times, the fracture mechanics theory has been extended to fields beyond purely structural applications, such as the field of micro electro-mechanical systems [11,12,13] and lithium-ion batteries (LIBs) [14,15,16].
Nowadays, LIBs technology is experiencing rapid growth, driven by the decarbonization challenge and the spread of LIBs in new fields like automotive, and stationary energy storage [17,18,19,20,21,22].It is demonstrated that mechanical fracture is among the main causes of aging and performance decay in LIBs [23].As a result, fracture in LIBs is receiving significant attention from the scientific community, especially with the aim of developing LIBs with long life cycles and high performance.
The basic working principle of lithium-ion batteries is the movement of electrons in the external circuit (current flow), balanced by the movement of lithium ions between electrodes of opposite polarity and the insertion in their structure, driven by electrochemical reactions.The electrodes are composite materials, among other components, the most important material is the active material where electrochemical reactions take place.The active material is a particulate matter, so when lithium ions get into the electrode, they intercalate into these particles, occupying the interstices in the crystalline structure of the active material.Since lithium ions are much smaller than active material particles (≈ 10 μm in radius) they diffuse within the active material particles.Then, lithium concentration distribution within the particles is inhomogeneous, being higher in the particle core with respect to the surface during lithium extraction, and the opposite during insertion.The intercalation of lithium-ions in the crystalline structure of the electrodes causes the local deformation of the material, proportional to the concentration of lithium-ions.Then, strain mismatch arises in the particles due to the inhomogeneous swelling caused by the inhomogeneous distribution of lithium concentration.This differential strain results in the so-called diffusioninduced stress (DIS), leading to crack propagation in electrode microstructure ultimately [24,25,26,27,28,29,30].Finally, crack propagation causes the capacity fade and impedance rise in LIBs because it hinders the passage of lithium -ions due to the isolation of some portions of active material, and it triggers undesired side reactions on newly created crack surfaces (solid electrolyte interface growth), consuming lithium ions [15,16].
Fracture problem in LIBs is usually solved numerically with FEM in the literature because of two reasons: firstly, a coupling exists between the mechanical and electrochemical domains which makes the problem multi-physics, secondly, stress distribution on the crack surfaces is nonhomogeneous and no analytical solutions are available for SIF computation.However, numerical computation with FEM is often time-consuming, then it is unsuitable for integration into a damage model for real-time scenarios.
Two methodologies are reported to compute SIF in case of an arbitrary stress distribution over the crack surface, the weight factor method [3,15] and the Green's function method [31].The weight function method is particularly difficult to implement.Green's function method is easier to handle, but it still requires solving integrals.For this reason, a more easy and analytical way to compute SIF is pursued, on the basis of geometric factors.
This work provides a general analytical methodology for calculating the SIF resulting from any arbitrary stress profile on the surfaces of the crack.This methodology is effective for accurately and rapidly assessing the fracture level in LIBs electrodes.First, a polynomial function is employed to approximate the arbitrary stress profile.Then, the SIF is calculated by applying the principle of superposition of effects, involving a summation of products between each term of the polynomial stress function and their corresponding geometric factors.Geometric factors for spheres with cracks at both the center and the surface are obtained from a FEM model.Subsequently, the SIF in graphite particles of LIBs electrodes is computed using these geometric factors.The analytical results are validated against the results of the multi-physics FEM fracture 1306 (2024) 012009 IOP Publishing doi:10.1088/1757-899X/1306/1/0120093 model.Finally, the methodology is used to assess the SIF and the fracture likeliness as a function of different electrode design solutions, at different current rates.
The article is organized as follows.Section 2.1 provides the main Equations solving the mechanics of LIBs.The fundamentals of LEFM theory with special regard to multi-physics applications are provided in Section 2.2.Section 2.3 describes the analytical methodology for SIF computation due to arbitrary stress profile on the surfaces of the crack.The geometric factors for spheres with cracks located at the center and on the surface are provided in 3.1 and the results of the application of the analytical procedure to the case of graphite particles in LIBs are discussed in Section 3.2.Section 3.3 reports the influence of electrode design parameters on fracture in LIBs and Section 3.4 describes the correlation between mechanical and electrochemical (capacity loss) damage.

Mechanics of lithium ion batteries
Active material particles of LIB electrodes undergo mechanical stress and strain due to lithium diffusion during battery operation.The problem is multi-physics because it involves computing the chemical concentration of lithium ions to get the stress.
The main equations of the electrochemical-mechanical model and their boundary conditions are reported in Table 1.The following assumptions are considered: the active material particle is spherical, the material is linear elastic, homogeneous, and isotropic [32].Then, equations are written as a function of the radial coordinate (r) thanks to the hypothesis of axisymmetry.
The chemical field is governed by the mass transfer equation (Equation 1), which describes the variation of lithium concentration c due to lithium diffusion, in analogy to the heat transfer equation.Stress-strain equations, i.e. constitutive (Equation 3), congruence (Equation 4), and equilibrium (Equation 5) equations, govern the mechanical field caused by the inhomogeneous lithium concentration, in analogy with the thermal-mechanical problem.
The constitutive equation shows that the total strain has two contributions, namely the elastic strain contribution (first term on the right-hand side of Equation 3) and the chemical strain contribution depending on the concentration c (second term on the right-hand side of Equation 3).It is pointed out that chemical strain ( Ω(c−c R ) 3 ), is equivalent to a thermal strain, where Ω 3 is replaced by the thermal expansion coefficient α and the concentration c is replaced by temperature T [26,28,33].
Once the diffusive problem (Equation 7) is solved and the lithium concentration c is got, the mechanical field is solved in turn.The equilibrium equation is written as a function of the displacement u using the constitutive (Equation 3) and congruence (Equation 4) equations.Then, the resulting equation is integrated twice and the displacement u (Equation 8) is got using the boundary conditions in Equation 6.The displacement u is replaced in the congruence equation (Equation 4) to compute radial and hoop strains, then, the constitutive equation (Equation 3) is used to obtain radial (Equation 9) and hoop stress (Equation 10).= 0 for t ≥ 0 ∂c ∂r r=R
In previous authors' work [26] it was shown that lithium transport is affected by DIS through the gradient of hydrostatic stress, which depends on the lithium concentration in turn.In this case, lithium transport is governed by Equation 11, which fully couples mechanical and chemical fields.
Where D is the diffusion coefficient, c R is the reference concentration, σ h = σr+2σc 3 is the hydrostatic stress, R g is the gas constant, and T is the temperature.
Equation 11 can be rewritten in the same form as the heat transfer equation and the solutions reported in Equation 7-10 can be employed.Then, an equivalent diffusion coefficient depending on lithium concentration substitutes the physical diffusion coefficient D in Equation 1, according to Equation 12 [26]. Where is the coefficient for the mechanical-chemical coupling, ν is the Poisson ratio, and E is the Young modulus.

Fracture mechanics theory for multi-physics problems
Pre-existent defects in electrode microstructure mainly propagate due to mode-I caused by tensile hoop stress, despite the numerous fracture mechanisms reported in the literature according to the specific active material [15,32,33].
According to LEFM theory, the stress field near the crack tip is infinite and can be expressed as a function of the SIF K, according to Equation 13 [5].
Where σ ij is the stress tensor, r c and θ are the polar coordinates with the origin at the crack tip, and f ij (θ) is a dimensionless shape function.
SIF can be computed as a function of a dimensionless geometric factor Y (Equation 14), which considers both the geometry of the specimen and of the crack (size, location, and shape).
Where a is the crack size and σ is the constant far-field stress, which is the stress unaffected by the crack presence.Numerous works in the literature provide various geometric factors for different cracks and body geometries [7,8,34].However, only solutions for simple cases have been developed, such as those involving constant or linear stress on the crack surfaces (caused by tension, bending, or their combination).
Tabulated geometric factors cannot be used when the stress distribution on crack surfaces is non-homogeneous, such as for stress resulting from diffusive fields (thermal or chemical).Then, SIF is usually computed using FEM from the J-integral developed by Rice [35] according to Equation 15, when the hypotheses of LEFM hold.
plane strain (15) The J-integral for a 2D fracture problem is the contour integral reported in Equation 16, considering an arbitrary counterclockwise path Γ around the crack tip (Figure 1).This contour integral is path-independent, meaning that its value is independent of the particular choice for the path Γ.
Where W = ij 0 σ ij d ij (Einstein notation) is the strain energy density, σ ij and ε ij are the components of the stress and strain tensors, u is the displacement vector, t is the traction vector acting on the path Γ and its components are t i = σ ij n j , n j is the versor normal to Γ, x 1 and x 2 are the coordinate directions, and ds is the element length along the path Γ.
The standard form J-integral does not keep the path-independence property in case of multiphysics fracture problems, such as thermal-mechanical [36,37] or mechanical-diffusive fracture problems [38,39,40].The modified version of path-independent J-integral for these multi-physics problems is reported in Equation 17.
Where Λ is the area enclosed within the path Γ, σ P = [σ 11 , σ 22 , σ 33 ] is the principal stresses vector, ε MF is the strain vector caused by the multi-physics fields, i.e. the chemical ( Ω(c−c R ) 3 ) or thermal field (α(T − T R )).The additional area integral (the second term on the right-hand side of Equation 17) is added to the standard Rice's expression of J-integral (the first term on the right-hand side of Equation 17) to ensure that the energy balance is satisfied and the J-integral is path-independent.

Analytical computation of SIF
The principle of superposition of effects can be exploited to analytically calculate the SIF as a function of geometric factors even when the stress distribution on crack surfaces is arbitrary.First, a polynomial of grade n is used to approximate the arbitrary nominal stress distribution (the stress unaffected by the crack presence) at the location of the crack, as reported in Equation 18.
Where σ i are the polynomial coefficients, and x is the coordinate along the crack surface as reported in Figure 2).
Then, Equation 18 is combined with Equation 13, and the principle of superposition of effects is applied, resulting in the generic expression of SIF as reported in Equation 19.
x a (x) = 0 +...+ n x n Where Y i is the geometric factor corresponding to the i-th term of the polynomial function approximating the stress for the specific crack configuration.
Referring to Figure 3, the procedure used in this work for the analytical computation of SIF in active material particles with central and superficial cracks, and for the validation of the results, is summarized below.• Geometric factor computation.The computation of each geometric factor Y i for each ith grade of the polynomial function approximating the stress profile is performed considering different ratios a/R, where R is the radius of the sphere.First, a FEM model of a sphere with cracks located at the center and the surface is built in Ansys Mechanical APDL.Then, a pressure load equal to the single grande of the polynomial stress function (σ i x i ) is applied on the nodes of the crack surface.Then, the SIF (K F EM i ) is computed from the J-integral (Equation 15) using the CINT command in Ansys.Finally, the geometric factor Y i is computed according to Equation 20.
This procedure is carried out just one time, once geometric factors are computed, they can be applied to different stress distributions and crack/body sizes.• Analytical SIF computation in case of DIS.The nominal hoop stress distribution over the crack region is computed according to the equations of the electrochemical-mechanical model reported in Table 1.Then, the hoop stress is fitted with a polynomial, and SIF is computed using Equation 19, using the geometric factors computed at the beginning.• Numerical SIF computation in case of DIS.SIF is computed using a FEM model built in Ansys Mechanical APDL.Lithium concentration in the active material particle is computed according to the equations of the electrochemical-mechanical model reported in Table 1.Then, the thermal analogy [32,33] is used to convert the concentration into an equivalent temperature.Finally, the equivalent temperature is mapped on structural nodes of the FEM model, and SIF is obtained from the J-integral through the CINT command (Equation 15).• Validation.The results of the analytical computation of SIF are compared with the ones resulting from numerical computation with FEM.
A detailed discussion on the strategy used to build the FEM model and mesh is not reported here but a more comprehensive explanation is present in previous authors' works [33,32].

Results and Discussion
This section provides the results of the computation of the geometric factors for a sphere with a crack located at the center and surface.Subsequently, the analytical computation of SIF due to the hoop stress caused by lithium diffusion in active material particles is performed.Finally, a comparison between the analytical and FEM model results is performed.

Geometric factor results
In this work, a polynomial stress distribution of grade 6 is considered, leading to seven geometric factors.The geometric factors computed are general and applicable to any polynomial stress distribution.In different scenarios, when the stress distribution can adequately be represented by a lower-grade polynomial, higher-order geometric factors can be neglected accordingly.
Figures 4 and 5 report the geometric factors for sphere with central and superficial cracks, respectively.The geometric factors exhibit quadratic dependency of the normalized crack length a/R.The analytical expressions for these geometric factors are provided in Tables 2 and 3.

Comparison between analytical and numerical SIF results in case of DIS
Analytical and numerical SIF results caused by lithium diffusion in active material particles of LIB electrodes are compared in this section.20, and the red dots correspond to the fitting using quadratic functions (p( a R ) 2 + q( a R ) + r).The equations in Section 2.1 are used to analytically solve the electrochemical-mechanical problem, including the coupling mechanical-chemical coupling.The galvanostatic operation is considered, corresponding to a constant flux on the particle surface (Equation 2), which depends on the current I delivered by the battery.The case study is a graphite particle, considering the following physical parameters: Young modulus E = 15 GPa, Poisson ratio ν = 0.3, partial molar volume Ω = 4.2×10 −6 m 3 /mol, diffusion coefficient D = 2×10 −14 m 2 /s, maximum concentration c max = 2.9155 × 10 4 mol/m 3 , temperature T = 298 K. Tensile hoop stress occurs at the particle centre during lithiation, vice-versa it occurs on the surface during delithiation, then the particle with central crack is simulated during lithiation, on the other hand, particle with superficial crack is simulated during delithiation.Normalized current rate is considered, then the C-rate equal to 1 corresponds to the current enables complete filling (or emptying) of the active material particle within one hour.In addition, SOC is the percentage fraction of lithium concentration with respect to the maximum concentration (SOC = c cmax ).Figures 6b,d compare SIF values computed analytically and numerically.The results obtained with the two methodologies are in good agreement, demonstrating the validity of the proposed analytical procedure for SIF computation due to a generic stress distribution over the crack.
SIF depends on hoop stress and crack size, according to Equation 19.The distribution of hoop stress within the particle at 50% SOC and different C-rates is shown in Figure 6a during lithium insertion and in Figure 6c during lithium extraction.As also demonstrated in previous authors' works [26,28,27], higher C-rates cause higher concentration gradients within active material particles, giving rise to higher stress.Similarly, SIF values increase with the increase of C-rate, as shown in Figure 6b,d.Then, higher C-rate is more detrimental because it increases fracture and accelerates LIB degradation, as confirmed by previous authors' work [32] and experimental measurements [41,42].On the other hand, as the crack length increases, the SIF increases accordingly, for the crack length below a/R = 0.5 for the particle with central crack and a/R = 0.15Then, the SIF decreases as the crack length increases.This occurs because the hoop stress distribution varies along the radial coordinate.Indeed, the positive (tension) hoop stress decreases and becomes negative (compressive) when moving from the particle core to the surface during lithium insertion (Figure 6a), and from the surface to the core during extraction (Figure 6c).Then, the crack tip of longer cracks undergoes lower hoop stress because it is closer to the less tension/compressive region.Then, referring to Equation 19, the SIF increases with the increase in crack length as long as the increase in crack length balances the decrease in the hoop stress.
Several works in the literature compute the SIF analytically in order to avoid multi-physics FEM simulations with high computational cost, especially when the aim is to integrate the fracture results into a degradation model for the estimation of capacity fade during battery operation [43,44,45,46,47].In these works, SIF is directly computed as K = Y σ √ a, making two approximations: (a) constant far-field stress is used, neglecting that the stress distribution on crack surfaces is not uniform; (b) the geometric factor of the plate undergoing uniaxial stress is used (Y 0 = 1.12 √ π [3]) instead of using the one relative to the spherical case.Figure 7a provides an insight into the error regarding the first approximation.The SIF for spherical particle with central crack computed with the analytical procedure proposed in this work is compared with the SIF computed as K = Y 0 σ √ a, where σ is equal to the stress in the uncracked particle at the crack tip location, and Y 0 is the geometric factor for constant stress distribution on the crack surfaces (Y 0 in Figure 4).The comparison demonstrates that the non-homogeneous stress distribution less affects the SIF for shorter On the other hand, SIF obtained considering uniform stress distribution decrease faster as the crack length increases, eventually reaching zero when the stress at the crack tip becomes compressive.Figure 7b compares SIF values obtained with the analytical computation methodology proposed in this work and the one computed with the approximation of the plate's shape factor and constant superficial stress, which is common in the literature.A particle with a superficial crack is considered, and the far-field stress σ is the hoop stress on the surface of the particle without the crack.The results show a great difference between SIF values, which still increases as a/R increases.Indeed, the decrease in SIF as the crack length increases is not observed when only focusing on the superficial hoop stress, thereby neglecting the complete stress distribution along the crack surfaces.Additionally, the geometric factor of the plate significantly differs from that of the sphere.
The approximations usually made in the literature affect the accuracy of the SIF computation,

Central crack
Super cial crack and then, the results of the degradation model assessing the decay of LIB performance during its usage.On the contrary, the analytical procedure developed in this work is as accurate as performing multi-physics FEM fracture model simulation.The advantage is that analytical computation is significantly faster and requires less computational resources as compared to the FEM model (fractions of a second for the analytical computation compared to minutes/hours for the FEM model) because it just a sum of products to compute the SIF.Then, the proposed analytical methodology for SIF computation is the optimal compromise to assess the impact of fracture on battery damage, balancing accuracy with computational cost.

Impact of fracture on the electrode design
The effect of the electrode design parameters on fracture is quantified in this section, such as the size of the active material particles and the thickness of the electrode layers.In this way, it is possible to give electrode design guidelines to limit fracture and mechanical degradation in LIBs.
Electrode design parameters affect the electrochemistry of the LIB, namely the current flowing into each elementary cell (sequence of current collectors, cathode, separator, anode), then the lithium flux on the surface of active material particles, as well as the resulting lithium concentration and stress distributions.First, lithium flux on particle surfaces (J flux ) is computed by solving the electrochemical model of the entire LIB as described in previous authors' works [32], varying the electrode design parameters.Then, equations in Section 2.1 are solved using the computed lithium flux in the boundary condition reported in Equation 2.
The effect of the active material particle size on the fracture behavior is reported in Figure 8.The concentration gradient within smaller particles is lower than bigger particles, then the former are less affected by fracture because of the lower stress.However, electrodes with smaller particles result in LIB with lower tap density and energy density.In addition, the manufacturing cost for electrodes with smaller particles is higher.Then, a trade-off is necessary between lifespan (which increases as active material particle size is reduced), electrochemical performance, and cost.Figure 9 reports the influence of electrode thickness on the SIF.Thicker electrodes allow to store greater amount of energy per unit of volume (or mass) because the number of elementary cells connected in parallel is lower, and then the amount of inactive components, i.e. separator and current collectors, is reduced.In this condition, considering the battery capacity and the current delivered by the battery as constant, a lower number of elementary cells results in a higher current flowing into each cell.This higher current results in higher lithium flux (J flux ) on active material particles surface, then steeper concentration gradient, higher stress, and SIF.On the other hand, thinner electrodes are generally used for high-power applications thanks to the lower internal resistance and greater current capabilities.Furthermore, thinner electrodes are preferred for higher C-rates than thicker electrodes because they are less affected by fracture, as shown in Figure 9.

Electrochemical-mechanical damage model
Fracture of electrode microstructure is among the main causes of LIB degradation [15].As summarized in the scheme reported in Figure 10, fatigue phenomena and cracks growth occurring during battery operation lead to capacity fade because of two reasons: (a) isolation of active material, which becomes inactive and no longer able to host lithium ions; (b) creation of new surfaces which are exposed to the electrolyte, giving rise to side reactions such as the formation of solid electrolyte interphase (SEI) layer.These reactions consume lithium ions, which cannot be cycled further.

Active material isolation and decrease of electronic conductivity
Hindering of the passage of lithium ions

Capacity fade
New surfaces exposed to the electrolyte

Repeated charge and discharge cycles
Active material Primary SEI New SEI The rate of capacity fade caused by the formation of the SEI layer on the cracks surfaces created during repeated charge/discharge cycles is computed according to Equation 21 [15].
Where Q is the capacity, N the number of cycles, A c is the crack surface and X is a constant taking into account the amount of lithium ions consumed during the SEI formation reaction and the rate of the reaction.The increase of the crack surface with the number of cycles ( dAc dN ) is computed according to Equations 22a-b [15].
Where a is the crack length, C and m are material coefficients and Equation 22b is the Paris' law, correlating the increase of crack length ( da dN ) to the SIF range (ΔK).Equations 21-22 show that the degradation in LIB is a multi-physics phenomenon because the mechanical damage (the increase in crack length in Equation 22) is correlated with the chemical damage (SEI formation and loss of lithium inventory in Equation 21).
Since the capacity loss depends on the SIF, the analytical computation of SIF developed in this work will significantly speed up the capacity loss computation.Then, this work is the basis for the development of the algorithms assessing the capacity loss and the remaining useful life of LIB during usage, which will be the focus of the future research of the authors.

Conclusions
A generic methodology for the analytical computation of SIF in case of arbitrary stress profile on crack surfaces is presented in this work.The methodology can be effectively applied to assess fracture resulting from multi-physics phenomena, such as fracture arising in active material of LIBs electrodes due to lithium diffusion.In this case, an analytical solution for SIF is not available in the literature because the stress distribution on the crack surface is nonhomogeneous.Firstly, the arbitrary stress profile is expressed using a polynomial function, then the principle of superposition of effects is applied, and SIF is computed as a summation of products between each term of the polynomial stress function and their corresponding geometric factors.A FEM model is built to compute the geometric factors for a sphere with a crack at the center and surface.The expression of the geometric factors is generic and can be used for any kind of stress profile and crack/radius geometrical sizes.The developed methodology is applied for the analytical computation of SIF caused by diffusion induced stress arising in graphite particles during lithium insertion and extraction.The analytical results are compared to the numerical results, showing good agreement.SIF due to DIS is usually computed in the literature using multi-physics FEM model, which is time-consuming (the simulation time can range between minutes to some hours), or analytically.In the latter case, the computation is inaccurate because the non-homogeneous stress distribution on crack surfaces is neglected, and the geometric factor of plates is often used, which is far from the case of spherical geometry.On the contrary, the proposed analytical computation of SIF is as accurate as performing FEM fracture model simulations, but it is less time-consuming because the SIF computation consists of just a sum of products, running in fractions of seconds.The methodology is effectively applied to assess the influence of electrodes parameters on fracture, in order to give electrode design guidelines limiting fracture and then LIB degradation.It is observed that thinner electrodes with smaller active material particles are preferred from a fracture mechanical point of view because they result in lower SIF values.Finally, the developed analytical methodology can be effectively used to describe fracture in active material particles with LIB cycling, being eventually integrated into a degradation model for the estimation of LIB performance decay in real-time.

Figure 1 .
Figure1.Scheme of the crack area, where Γ is an arbitrary counterclockwise path surrounding the crack tip, Λ is the area enclosed by the path Γ, x 1 , x 2 are the coordinate directions and n is the versor normal to Γ.

Figure 2 .
Figure 2. Crack surface with polynomial stress distribution, where x is the coordinate along the crack surface and a is the crack length.

Figure 3 .
Figure 3. Methodology for the analytical computation of SIF, where the validation stage is highlighted in green.

Figure 6 .
Figure 6.Results of DIS and fracture model.The graphite particle is simulated at different C-rates and results correspond to 50% SOC.Distribution of hoop stress (a) during lithiation and (c) during delithiation.SIF values for different a/R for particle (b) with a crack at the center during lithiation and (d) a crack at the surface during delithiation.

Figure 7 .
Figure 7.Comparison between SIF values resulting from literature and obtained with the procedure presented in this work.(a) Graphite particle with central crack during lithiation at 1C and 50% SOC, where σ c,tip is the hoop stress in the uncracked particle at the crack tip location.(b) Graphite particle with superficial crack during delithiation at 1C and 50% SOC, where σ c,sup is the hoop stress on the particle surface.

Figure 8 .
Figure 8. Influence of active material particle size on SIF during (a) charge (lithium insertion in graphite particle) and (b) discharge (lithium extraction from graphite particle).

Figure 9 .
Figure 9. Influence of the electrode thickness on SIF during (a) charge (lithium insertion in graphite particle) and (b) discharge (lithium extraction from graphite particle).

Figure 10 .
Figure 10.Sketch of the damaging mechanisms triggered by fracture during LIB operation.

Table 2 .
Coefficients of the quadratic function p( a R )2+ q( a R ) + r for geometric factors of sphere with central crack.

Table 3 .
Coefficients of the quadratic function p( a R ) 2 + q( a R ) + r for geometric factors of sphere with superficial crack.