Phase Field Modeling of Chemomechanical Fracture of Intercalation Electrodes: Role of Charging Rate and Dimensionality

TL;DR

This study uses phase-field modeling to analyze how charging rate and dimensionality influence crack formation in Li-ion battery cathodes, revealing different fracture behaviors in 2D and 3D geometries.

Contribution

It introduces a thermodynamically consistent phase-field approach to quantify chemomechanical fracture in battery particles, including a derived power-law for critical flaw size dependence on charging rate.

Findings

Critical flaw size depends on charging rate via a power-law.

2D and 3D crack propagation behaviors differ significantly.

3D cracks tend to bifurcate and limit inward penetration.

Abstract

We investigate the fracture of Li-ion battery cathodic particles using a thermodynamically consistent phase-field approach that can describe arbitrarily complex crack paths and captures the full coupling between Li-ion diffusion, stress, and fracture. Building on earlier studies that introduced the concept of electrochemical shock, we use this approach to quantify the relationships between stable or unstable crack propagation, flaw size, and C-rate for 2D disks and 3D spherical particles. We find that over an intermediate range of flaw sizes, the critical flaw size for the onset of crack propagation depends on charging rate as an approximate power-law that we derive analytically. This scaling law is quantified in 2D by exhaustive simulations and is also supported by 3D simulations. In addition, our results reveal a significant difference between 2D and 3D geometries. In 2D, cracks propagate deep inside the particle in a rectilinear manner while in 3D they propagate peripherally on the surface and bifurcate into daughter cracks, thereby limiting inward penetration and giving rise to complex crack geometries.

Tables

Table 1. Table 1: Material properties of \ce LiMn2O4 for numerical simulation of chemo-mechanical fracture  [ 17 ] .

Property	Symbol	Units	Value
Elastic Modulus	$E$	$N{m}^{-2}$	$2\times {10}^{11}$
Poisson Ratio	$\nu$	-	$0.3$
Fracture Toughness	$G_c$	$N{m}^{-1}$	$100$
Diffusivity	$D_0$	$m^2{s}^{-1}$	$2.2\times {10}^{-13}$
Maximum Concentration	$c_{\max }$	$\mathrm{mol}{m}^{-3}$	$2.37\times {10}^4$
Misfit Strain Constant	${ϵ}_0$	$m^3{\mathrm{mol}}^{-1}$	$1.09\times {10}^{-6}$
Density	$\rho$	$\mathrm{kg}{m}^{-3}$	$4.28\times {10}^3$
Temperature	$T$	$K$	$3\times {10}^2$
Dimensionless Expansion Coefficient	$\beta$	$m/m$	$0.025$
Griffith Length Scale	$l_G$	$\mu m$	$5\times {10}^{-4}$

Equations

$$\mathcal{F}(u,c,\Gamma )={\mathcal{F}}_{el}(u,c,\Gamma )+{\mathcal{F}}_c(c)+{\mathcal{F}}_{\Gamma }(\Gamma )$$

$${\mathcal{F}}_{el}(u,c,\Gamma )={\int }_{\Omega \setminus \Gamma }\mathcal{W}(u,c)dx$$

$${\mathcal{F}}_c(c)={\int }_{\Omega }{f}_c(c)dx$$

$$f_c=c_{\max }\mathcal{R}T\left[\frac{c}{c_{\max }}\ln \left(\frac{c}{c_{\max }}\right)+\left(1-\frac{c}{c_{\max }}\right)\ln \left(1-\frac{c}{c_{\max }}\right)\right]$$

$${\mathcal{F}}_{\Gamma }(\Gamma )=G_c{\mathcal{H}}^{n-1}(\Gamma )$$

$$\mathcal{F}(u,c,\varphi )={\mathcal{F}}_{el}(u,c,\varphi )+{\mathcal{F}}_c(c)+{\mathcal{F}}_{\varphi }(\varphi )$$

$${\mathcal{F}}_{el}(u,c,\varphi )={\int }_{\Omega }g(\varphi )\mathcal{W}(u,c)dx$$

$${\mathcal{F}}_{\varphi }(\varphi )=\frac{G_c}{4C_{\varphi }}{\int }_{\Omega }\left(\frac{w(\varphi )}{\xi }+\xi {\left|\nabla \varphi \right|}^2\right)dx$$

$$\frac{\partial c}{\partial t}=-\nabla \cdot J$$

$$J=-a(\varphi )M(c)\nabla \mu$$

$$\mu =\frac{\delta \mathcal{F}}{\delta c}$$

$$\frac{\delta \mathcal{F}}{\delta c}=\frac{d{f}_c}{dc}(c)-{ϵ}_{0}g(\varphi ){\sigma }_{kk}$$

$$\frac{\partial c}{\partial t}$$

$$\psi$$

$${J|}_{{\partial }_f\Omega }=\frac{i}{\mathcal{F}}\equiv \hat{J}$$

$$t_C=\frac{c_{\max }V}{\hat{J}A}$$

$${\mathcal{C}}_r=t_D{C}_r=\frac{t_D}{t_C}$$

$$\mathcal{J}={\mathcal{C}}_r\frac{\bar{V}}{\bar{A}}=\frac{t_D}{t_C}\frac{\bar{V}}{\bar{A}}=\frac{\hat{J}R}{c_{\max }{D}_0}$$

$${\int }_{\Omega }\left(\frac{c_t-c_{t-1}}{\delta t}\right)\tilde{c}dx+{\int }_{\Omega }a({\varphi }_{t-1}){\nabla }_{\Theta }c\cdot \nabla \tilde{c}dx$$

$$+{\int }_{\Omega }a({\varphi }_{t-1})M(c_t){\nabla }_{\Theta }\psi \cdot \nabla \tilde{c}dx+{\int }_{{\partial }_f\Omega }\mathcal{J}\tilde{c}ds=0$$

$${\int }_{\Omega }\left[{\psi }_t+\beta g({\varphi }_{t-1}){\sigma }_{kk}(u_{t-1},c_t)\right]\tilde{\psi }dx=0$$

$${\int }_{\Omega }g({\varphi }_t){\sigma }_{ij}(u,c)e_{ij}(\tilde{u})dx=0$$

$${\int }_{\Omega }\left[\frac{dg}{d\varphi }({\varphi }_t)\mathcal{W}(u,c)\right]\tilde{\varphi }dx$$

$$+\frac{G_c}{4C_{\varphi }}{\int }_{\Omega }\left[\frac{1}{\xi }\left(\frac{dw}{d\varphi }({\varphi }_t)\right)\tilde{\varphi }+2\xi \nabla \varphi \cdot \nabla \tilde{\varphi }\right]dx=0$$

$${\tau }_{\varphi }\frac{\partial \varphi }{\partial t}=\frac{1}{c_{\max }\mathcal{R}T}\frac{\delta \mathcal{F}}{\delta \varphi }$$

$$l_G=\frac{G_c}{\mathrm{E}}$$

$$a_{0,\min }\sim l_G(\beta {\mathcal{C}}_r{)}^{-2}.$$

$$l_c=\frac{l_G}{{\beta }^2}$$

$$C_{r,\max }=\frac{1}{t_{c,\min }}=\frac{D_0}{R^2}{\left(\frac{A{l}_c}{a_{0,\min }}\right)}^{1/2}$$

$$\frac{{\mathcal{F}}_{\varphi }}{G_c{R}^2}\simeq \frac{{\mathcal{H}}^2(\Gamma )}{R^2}$$

$${\nabla }^4\mathcal{A}$$

$$J$$

$$\frac{\delta F}{\delta \bar{c}}$$

Full text

1 Introduction

With the demand for electric vehicles and hand-held electronics on the rise, research on rechargeable batteries and, specifically lithium-ion batteries, becomes increasingly important. The need to understand the failure mechanism of these batteries is essential for increasing their life span. Chemo-mechanical failure is one of the primary modes of degradation. The fracture of cathodic and anodic particles due to intercalation-induced stresses has been extensively studied experimentally [1, 2, 3]. The creation of new fracture surfaces impairs the performance of the batteries due to the loss of electrical contact [4, 5] and the creation of solid electrolyte interfaces (SEI) that promotes the irreversible loss of lithium (Li) ions [6, 7].

On the theoretical side, problems arising from the interplay of diffusion and mechanics have been long considered in the literature. Prussin [8] and Lawrence [9] were among the first to study the creation and motion of dislocations due to diffusion induced misfit strains. Subsequently, Liu et al. [10] studied the attraction of corrosive solutes to the crack tip. In the context of chemo-mechanical fracture, Huggins and Nix [11] studied the initiation of fracture due to the intercalation-driven misfit stresses using a simple 1D model of a thin film bilayer. In this bilayer geometry with a rigid substrate, a constant misfit strain caused by Li intercalation in the thin film is sufficient to create cracks by a mechanism similar to thermal expansion [12]. Therefore, using a classical Griffith criterion [13], Huggins and Nix were able to derive a critical film thickness for fracture. The extension to free-standing particles was subsequently considered in several studies [14, 15, 16, 17]. Unlike in a thin film constrained on a substrate, a uniform concentration does not create stresses in a free-standing particle. However, due to the finite time to diffusively homogenize the Li concentration inside a particle, a concentration gradient that produces stresses can nonetheless be created when the Li ion flux through the particle boundary, *i.e., *the charging rate (C-rate), is sufficiently large. By analogy with thermal shock, Woodford et al. [17, 18, 19] coined the term “electrochemical shock” to describe this mode of C-rate dependent fracture  [20]. A consistent model that considered the effect of intercalation induced stresses on diffusion was used by Christensen et al. [14, 15] to obtain a failure criterion solely based on the magnitude of the resulting stresses. Bhandakkar and Gao [16] investigated the initiation of a periodic array of equidistant cracks in a thin strip under an imposed constant galvanostatic flux. Using a cohesive zone model and neglecting the effect of fracture on the concentration field, they derived a scaling law relating the largest “safe” strip thickness, below which cracks do not initiate, and C-rate. Woodford et al. [17] investigated the propagation of fracture from an initial radial penny-shaped flaw on a spherical particle. Their calculation of the stress intensity factors at the crack tip was carried out in a simplified geometry, also neglecting the effect of fracture on the concentration field. Their results show that, at a given C-rate, both continuous and abrupt propagation modes are possible for different initial flaw sizes. They further show that for a given flux, there exists a largest safe particle size that does not fracture for any flaw size. Their numerical results indicate that this critical particle size scales as a power-law of charging rate (C-rate).

Even though those studies have yielded quantitative predictions of the dependence of safe particle size on C-rate [16, 17], they do not consider the full coupling between elasticity, fracture, and diffusion. In addition, penny-shaped cracks are assumed to remain coplanar as they penetrate a 3D spherical particle [17]. In practice, more complex non-coplanar crack patterns may develop that depend on C-rate. The goal of this article is to investigate the fracture of Li-ion battery cathodic particles using a thermodynamically consistent phase-field approach that captures the full coupling between elasticity, fracture, and diffusion, and that can describe arbitrarily complex crack paths. We exploit those advantages to quantify the relationship between crack propagation, flaw size, and C-rate, and to describe for the first time complex 3D crack patterns. Due to their variational formulation, phase-field models of fracture [21, 22, 23], offer a unique methodology to tackle the simulation of chemo-mechanical crack growth. These models have been validated by theoretical analyses [24] and comparisons of predicted and observed crack paths in non-trivial geometries [25]. They have been used to reproduce complex experimental observations in brittle fracture including thin-film fracture [26], thermal fracture [20], mixed mode fracture [27], dynamic fracture [27, 28, 29], fracture in colloidal systems [30], ductile fracture [31, 32, 33, 34], and fatigue fracture [35, 36, 37]. Given their potential over the past few years, researchers have extended the use of these models to chemo-mechanical fracture in battery particles [38, 39, 40, 41]. This approach has already been used [40, 41] to corroborate findings of Woodford et al. [17] such as the existence of unstable and stable crack propagation as a function of initial flaw size and to describe simple 3D crack patterns.

In this article, we extend and generalize the results of previous studies [16, 17, 40] to account for the effect of the crack length on the failure of 2D circular disks and 3D spherical particles. By performing an exhaustive series of 2D simulations for different flaw sizes and C-rates for a fixed particle radius, we identify three regimes of fracture propagation where (I) large flaws comparable to the particle size do not propagate due to insufficient driving stresses, (II) for intermediate flaw sizes, the critical flaw size scales as an approximate power-law function of C-rate with an exponent that we derive analytically, (III) for very small flaw sizes the C-rate required for propagation diverges resulting in a flux-independent minimum flaw size. Next, we obtain a scaling law relating the safe particle size, computed with a fixed flaw size to particle radius ratio, to the C-rate. Furthermore, we show that the topology of fracture changes profoundly in 3D compared to 2D. We find that, unlike in our 2D studies and previous 3D studies [17, 40], where coplanar cracks penetrate radially towards the center of the particle, 3D cracks remain mostly superficial and branch to tile complex crack patterns on the particle surface. Our results indicate that, despite this difference, the dependence of critical flaw size on C-rate follows a similar power-law scaling as in 2D.

This article is organized as follows. In Section 2, we outline a thermodynamically consistent formulation of chemo-mechanical concentration and stress evolution and fracture. In Section 3, we carry out a scaling analysis of the governing equations and define a subset of key dimensionless parameters. In Section 4.1, for a generic set of material parameters for \ceLiMn2O4, we present the results of an extensive set of numerical simulations in 2D circular disks and examine the propagation of radial flaws of different sizes. We investigate the influence of C-rate and initial flaw size on crack stability, generalizing the findings of [17, 40]. We extend our analysis to maximal C-rates in 4.2 and show that there exists a safe particle size regardless of the initial flaw size that can be predicted based on material properties including fracture energy, elastic modulus, and magnitude of misfit strain. We finally extend our analysis to 3D spherical particles with a single radial penny-shaped surface flaw in Section 4.3. Lastly, in Section 5, we summarize our main findings and point out possible future extensions.

2 Formulation

We define the total free energy $\mathscr{F}(u,c,\Gamma)$ for a domain $\Omega\subset\mathbb{R}^{n}$, containing the crack set $\Gamma\subset\Omega$, for displacement $u$ and concentration $c$

[TABLE]

where $\mathscr{F}_{el}$ is the elastic energy, $\mathscr{F}_{\Gamma}(\Gamma)$ is the energetic cost of fracture and $\mathscr{F}_{c}$ is the free energy due to the intercalation of Li ions. We can write the elastic energy as

[TABLE]

where we define the elastic strain energy as $\mathcal{W}(u,c)=\sigma_{ij}\varepsilon_{ij}/2$ in which, the elastic strain is defined as $\varepsilon_{ij}(u,c)=e_{ij}(u)-\epsilon_{0}(c-c_{0})\delta_{ij}$ where $\epsilon_{0}$ is the volume expansion coefficient. Furthermore, we define the linear strain $e_{ij}(u)=(u_{i,j}+u_{j,i})/2$, and the Cauchy stress tensor $\sigma_{ij}(u,c)=\mathcal{C}_{ijkl}\varepsilon_{kl}$. Moreover, for Lame’s constants $\lambda,\mu$ the isotropic elasticity tensor is written as $\mathcal{C}_{ijkl}=\lambda\delta_{ij}\delta_{kl}+\mu(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})$.

We write the free energy of the Li ions intercalating in the host lattice as [42]

[TABLE]

where

[TABLE]

is the entropy of mixing for an ideal binary solution, $c_{\max}$ is the maximum concentration achievable when all accommodating sites are filled, $\mathcal{R}$ is the gas constant, and $T$ is the absolute temperature.

In the spirit of brittle fracture, we write the energetic cost of creating fracture surfaces as

[TABLE]

where $G_{c}$ is the energy required to create a unit area (unit length in 2D) of new cracks, $\mathcal{H}^{m}$ is the $m$–dimensional Hausdorff measure (*i.e., *$\mathcal{H}^{2}(\Gamma)$ is the aggregate area and $\mathcal{H}^{1}(\Gamma)$ is the aggregate length of cracks $\Gamma$ in three and two dimensions, respectively).

2.1 Phase-field model

We use the phase-field model to approximate the sharp interface free energy (1) by introducing a fracture phase field $\phi$ and an associated length-scale $\xi$. Roughly speaking, as $\xi\to 0$, the displacement field minimizing (6) converges to that of minimizing (1), the field $\phi$ converges to 1 almost everywhere and goes to zero near the cracks. In this article, we treat the length-scale $\xi$ as a regularization parameter to study the sharp-interface limit of the phase-field model that reduces to classical linear elastic fracture mechanics [24]. We write the approximate free energy replacing $\mathscr{F}_{el}(u,c,\Gamma)$ by $\mathscr{F}_{el}(u,c,\phi)$ and $\mathscr{F}_{\Gamma}(\Gamma)$ by $\mathscr{F}_{\phi}(\phi)$ as

[TABLE]

with the elastic energy

[TABLE]

and the energetic cost of fracture

[TABLE]

where $C_{\phi}=\int_{0}^{1}\sqrt{w(\phi)}\,d\phi$ is a scaling constant. In the past decade, there has been a growing trend in studying a broad class of rate independent gradient damage models in the form of (6) [43, 44]. In this article, we use the Karma-Kessler-Levine model (KKL) [22, 24] defined using $g(\phi)=4\phi^{3}-3\phi^{4}$, $w(\phi)=1-g(\phi)$. This model allows us to follow the propagation of a fracture from a single flaw by prohibiting the initiation of new cracks in undamaged material (*i.e., *$\phi=1$).

2.2 Diffusion equation

Following the classical argument of continuity (mass conservation), we write the diffusion equation for concentration as [42]

[TABLE]

where $J$ is the flux of Li ions. We define the flux as the product of the mobility and the gradient of chemical potential

[TABLE]

where the mobility of Li ions in the host lattice $M(c)=m_{0}\,(c/c_{\max})(1-c/c_{\max})$ first increases and then decreases as a function of relative concentration $c/c_{\max}$. We define the chemical potential as

[TABLE]

in which $\dfrac{\delta\mathscr{F}}{\delta c}$ is the Fréchet derivative of free energy $F$ with respect to the concentration $c$ and can be written as

[TABLE]

replacing from above in (10) we finally get:

[TABLE]

We note that (13) represents a Fickian diffusion with the second term coupling to the mechanical hydrostatic stresses.

As a first estimate, we make the crack surface completely permeable to Li ion diffusion similar to [40] by introducing $a(\phi)=1$. We motivate this choice by noticing that the electrolyte will leak inside the newly created fracture surfaces thus making them permeable to ion transfer. The precise choice of $a(\phi)$ should be determined by further study of specific material and the interaction of the electrolyte and fracture surfaces and is out of the scope of this article.

Furthermore, to model galvanostatic charging, we write the flow of Lithium ions from the boundary as a given imposed flux $\hat{J}$:

[TABLE]

where $i$ is the surface current density and $\mathcal{F}$ is Faraday’s constant.

The galvanostatic boundary condition of form (15) is a first order approximation of ion transfer on the particle boundary and. A more general study can be done by prescribing the value of the flux on the surface at a given voltage as using the Butler-Volmer equation to model reaction kinetics on the cathodic surface of the particle [45]. We should note that, as for any diffusion process of a bounded field (here $0\leq c\leq c_{\max}$) with a flux boundary, this boundary condition cannot be maintained for any arbitrary value of $\hat{J}$ for an infinite time. In particular, $\hat{J}\gg 1$ at $t\ll R^{2}/m_{0}\mathcal{R}T$ a depleted boundary layer is created where the concentration at the boundary reaches zero and the flux cannot be maintained any longer.

3 Numerical implementation

3.1 Dimensional analysis

For the flux boundary condition (15), it is pertinent to introduce the nominal charging time as the time required to fill the volume of the particle $V$ with a surface flux $\hat{J}$ acting on surface area $A$ i.e.,

[TABLE]

It is also customary in Li-ion literature to introduce the so-called charging rate $C_{r}=1/t_{C}$, which is usually measured in $\mathrm{hr}^{-1}$ units. To perform the numerical simulations, we adimensionalize the spatial dimensions by the particle radius $R$, the concentration by $c_{\max}$, the time by the diffusion time $t_{D}=R^{2}/D_{0}$ where $D_{0}=m_{0}\mathcal{R}T$ is the diffusion constant, and the stresses by energy per unit volume $c_{\max}\mathcal{R}T$. We write dimensionless charging rate $\mathcal{C}_{r}$ as

[TABLE]

and the dimensionless flux as

[TABLE]

where $\bar{V}=V/R^{n}$ is the dimensionless volume of the particle, and $\bar{A}=A/R^{n-1}$ is the dimensionless surface area of the flux boundary. The the dimensionless charging rate $\mathcal{C}_{r}$ can also be understood intuitively as a mechanical loading parameter noticing that the driving force for crack propagation is the gradient of the concentration field in (7) that is controlled by the flux $\hat{J}$. As a result, for low dimensionless charging rates $\mathcal{C}_{r}<1$ (where the nominal charging time is long compared to the diffusion time $t_{C}\gg t_{D}$) the concentration will homogenize and thus creates no misfit stresses.

We should also highlight the important dimensionless numbers that uniquely define the simulations performed, namely the relative strength of the elastic energy compared to the chemical energy $\mathrm{E}/c_{\max}\mathcal{R}T$, Poisson’s ratio $\nu$, maximum misfit strain $\beta=c_{\max}\epsilon_{0}$, and the relative domain geometry *i.e., *radius and initial flaw size, compared to the Griffith length scale $R/(G_{c}/\mathrm{E})=R/l_{G}$ and $a_{0}/l_{G}$.

3.2 Governing equations

To implement our numerical simulations using the Galerkin finite element method we introduce the weak forms of the governing equations. The governing equations for the concentration diffusion is derived from its flow rule (13)-(14) by multiplying both sides with test functions and integrating by parts. We also incorporate an implicit time integration scheme to ensure the accuracy and stability of the integration.

[TABLE]

where $n$ is the surface normal to $\partial_{f}\Omega$, subscripts denote time steps with $\delta t$ as the time step size, and we define $\nabla_{\Theta}\{\circ\}=(1-\Theta)\nabla\,\{\circ\}_{t}+\Theta\nabla\,\{\circ\}_{t-1}$ as the implicit gradient operator associated with time-fraction $\Theta$. In all calculations in this paper we used $\Theta=0.5$ which corresponds to the midpoint method and results in a second order accurate and unconditionally stable time integration for concentration field $c$.

Moreover, since in practical systems the time-scale of elasticity and fracture propagation are orders of magnitude smaller than that of diffusion, we assume that they are instantaneous. In this setting, we seek the minimizers for the displacement field $u$ and the fracture phase-field $\phi$ for each time step $t_{i}$. Hence, the governing equations for displacement (*i.e., *elasticity) and fracture phase-field, are written as Euler-Lagrange equations of the total energy (1) for displacement field $u$ and phase field $\phi$

[TABLE]

where $c_{t}$ is the concentration given by the solution of (19)-(20).

3.3 The solution algorithm

The phase-field fracture method requires that the spatial resolution of discretization to resolve the characteristic approximation length $\xi$. The resulting problems are often very large and necessitate the use of a parallel programming paradigm and the complex numerical tools therein. Our implementation relies on the distributed data structures provided by libMesh [46] and for linear algebra on PETSc [47, 48]. On the other hand, we assume that elasticity and fracture are instantaneous and write their governing equations as the weak forms of Euler-Lagrange equation for minimizers of (6) with respect to displacement field $u$ and phase field $\phi$ respectively (see 3.2 for details). This is roughly equivalent to the limit of vanishing relaxation time $\tau_{\phi}\rightarrow 0$ assuming that the phase field $\phi$ follows Ginsburg-Landau gradient dynamics:

[TABLE]

We use a classical alternate minimization algorithm 1 [23] since the governing equations for elasticity and phase-field are only convex in either $u$ or $\phi$ when the other is kept constant [23]. It is also worth mentioning that to enforce irreversibility of fracture and ensure boundedness of phase field $0\leq\phi\leq 1$ and the relative concentration $0\leq c\leq 1$, we use a bounded reduced space Newton minimization scheme for the discrete energy provided in PETSc [47, 48].

4 Numerical results

In the following section, we focus on the numerical simulation of a cathodic particle at the time of charging. We first present our two-dimensional results for circular particles with a preexisting radial flaw on its surface under galvanostatic and potentiostatic boundary conditions. Subsequently, we analyze the fracture of spherical particles with penny-shaped cracks in three dimensions.

4.1 Chemo-mechanical fracture of circular particles: (I) galvanostatic (flux) boundary condition

The misfit strains generated during charging and discharging processes can lead to the creation and propagation of cracks in Li-ion battery particles. For a preexisting flaw on the surface of a cathodic particle, the removal of Li-ions during the charging process causes the outer layer of the particle to contract faster than its inner core. Therefore, for fast enough charging rates, the region of tensile stresses created in the outer periphery can activate surface defects creating cracks that will then propagate through the particle. In this section, we present the results of numerical simulation for the fracture of circular particles induced by the removal of Lithium ions during charging. Our goal is two fold: (i) to understand the activation and propagation of a preexisting flaw in a circular cathodic particle, (ii) 2D simulations also enable us to combine the results of many such simulations to give insight into critical parameters for the design of these particles.

Fig. 1 shows a schematics for this problem. We assume a preexisting radial crack $\Gamma_{0}$ of length $a_{0}$ and impose a dimensionless galvanostatic flux $\mathcal{J}$ (corresponding to the dimensionless charge-rate $\mathcal{C}_{r}=t_{D}/t_{C}$ according to (18)). As stated before, we treat phase-field length scale $\xi$ as a regularization of Griffith brittle fracture. Therefore for the Griffith length scale defined as

[TABLE]

we use $\xi>l_{G}$. For these numerical simulations we use a constant relative process zone size $\xi/R=1.25\times 10^{-2}$ for relative flaw size $a_{0}/R>0.1$ and use $\xi/a_{0}=5$ for $a_{0}/R<0.1$ for optimal use of computational resources. Table 1 summarizes the material properties corresponding to \ceLiMn2O4 used in our simulations.

To highlight the mechanism and modes of radial crack penetration in these particles, we first study three sample cases in Figs. 2–4. These sample results correspond to fracture of a $R/l_{G}=4.2\times 10^{4}$ and $a_{0}/l_{G}=2\times 10^{3},10^{4}$ initial radial flaws (corresponding to a $R=21\,\mu\mathrm{m}$ particle with $a_{0}=1,5\,\mu\mathrm{m}$ for material properties in Table 1) and using $\mathcal{C}_{r}=5.57,8.35$ dimensionless charging rates (circled gray in Fig. 7). In these simulations, we first compare two cases with different initial flaw lengths at the same dimensionless charging rate $\mathcal{C}_{r}$ and then study two cases where we keep $a_{0}$ constant and change the $\mathcal{C}_{r}$. As we stated before, our simulation results show that tensile hoop stresses are created on the periphery of these particles that can then drive the surface flaws to penetrate radially inside the particle (top row in Figs. 2–4). Figs. 2 shows that the crack propagation for the larger initial flaw $a_{0}/l_{G}=10^{4}$ under lower dimensionless charging rate $\mathcal{C}_{r}=5.57$ is continuous. However, the smaller initial flaw under the same charging rate $\mathcal{C}_{r}=5.57$ propagates abruptly jumping many process zone sizes (see third columns in Fig. 3). The abrupt propagation occurs in the context of Griffith fracture where the crack releases more energy as it propagates (*i.e., *for the energy release rate defined as $G=-\partial\mathscr{F}_{el}/\partial a$ at a frozen concentration (constant load), the crack is unstable if ${dG}/{da}>0$). Analogous results on abrupt propagation due to misfit strains are also reported in the context of thermally driven cracks. For example, Bahr et al. [49] explicitly calculate the energy release rate as a function of crack length for thermal-quenching-induced cracks.

A similar contrast between continuous and abrupt propagation is also predicted by Woodford et al. [17] where they explicitly calculate the mode-I stress intensity factor $K_{I}$ (and, by symmetry since $K_{II}\equiv 0$, the energy release rate $G$) for a radial penny-shaped crack in a spherical particle. Their calculations show that for some choices of particle size and initial flaw size ${dG}/{da}>0$; therefore, for such flaw sizes, propagation is unstable and abrupt. Concurring with our simulations, their calculations (figure 5. in [17] where their results correspond directly to Figs. 2–3) show that the flaws smaller than $a_{0}/l_{G}<4\times 10^{3}$ will propagate abruptly given our choice of parameters and vice versa.

The transition from abrupt to continuous propagation can also be understood in analogy with mechanical loading in standard fracture mechanics configurations. At lower charging rates where the concentration field has to penetrate on the scale of the particle size before the energy release rate reaches the fracture energy, the activation of the initial notch is analogous to a crack in half plane under constant far field opening stress that results in an unstable propagation. On the other hand, at high charging rates where the concentration field penetrated on the scale of the initial notch only, the problem resembles a compact specimen with the crack opening from the back that results in stable propagation. Our hypothesis is further verified, comparing Fig. 3 to Fig. 4. We can see that at the higher $\mathcal{C}_{r}$ (higher flux), the initial abrupt crack propagation is shorter. While surprising at first glance, we can explain the longer abrupt fracture propagation at lower $\mathcal{C}_{r}$, noticing that the flaw is activated earlier for the higher $\mathcal{C}_{r}$. Thus, at the time of the initial jump, the hoop stresses penetrate farther inside the particle for the lower flux providing more elastic energy to be converted to new fracture surfaces.

Different modes of fracture propagation are further demonstrated and quantified in Figs. 5–6 where we show the evolution of relative crack lengths $a/R$ and the dimensionless hoop stress in front of the crack at $r=R$ versus the charging time fraction ${t}/{t_{C}}$ for the simulations of a $R/l_{G}=4.2\times 10^{4}$ particle with $a_{0}/l_{G}=2\times 10^{3},~{}10^{4}$ preexisting radial flaws respectively and for different dimensionless charging rates $\mathcal{C}_{r}$ (including cases highlighted in Figs. 3–4). Similar to Fig. 2, the simulations for the larger $a_{0}/l_{G}=10^{4}$, depicted in Fig. 5, show a clear trend whereby the crack is activated $t/t_{C}\simeq 0.2$ and propagates smoothly. Unlike results presented in Figs. 5, the initial crack propagation in Fig. 6 is abrupt and decreases with higher $\mathcal{C}_{r}$. We should also note, in Figs. 5–6, that while the abrupt initial propagation is bigger for the smaller flux, the cracks extend farther into the particle for higher fluxes in line with our intuitive understanding that higher fluxes provide more energy. Similar observations were also made in [40] (see figure 15 in the reference) where for a $R/l_{G}=4.65\times 10^{10}$ particle containing $a_{0}/l_{G}=9.3\times 10^{8}$ initial crack they observed a larger initial abrupt propagation for lower fluxes and vice versa.

Our self-consistent simulations also allow us to observe the interaction of the stresses with the concentration field. We note that the tensile crack-tip stresses attract ions from its vicinity and results in crack tip enrichment. Many semi-analytic simulations, currently available in the literature, are based on the radial approximation of the concentration field [16, 17] (*i.e., *$c=c(r)$) which is calculated for an un-cracked particle. We study the crack-tip enrichment further in the A where in a simple setting we identify an enrichment length scale where its ratio to Griffith length scale, scales as the ratio of maximum misfit stresses to the chemical energy squared (see (37)).

With the propagation mechanism elucidated, we now turn our attention to obtaining design parameters for these particles. Fig. 7 shows a combined diagram for results of our simulations for $R/l_{G}=4.2\times 10^{4}$ particles containing initial flaws of different sizes. In this phase diagram, the circles mark activated cracks and cross marks depict those not activated at a given charging rate $\mathcal{C}_{r}$. As expected, longer radial cracks need lower charging rate to activate, but very long initial flaws do not propagate since the hoop stresses around the crack tip never grow large enough. We duplicated the simulations for a larger $R/l_{G}=2\times 10^{5}$ particle in Fig. 8.

We can make three important observations from Figs. 7–8; First the simulations show that there exists a safe charging rate $\mathcal{C}_{r,\min}\simeq 2,1$ for $r/l_{G}=4.2\times 10^{4},2\times 10^{5}$ particles respectively where flaws regardless of their size do not propagate. Secondly, for moderate dimensionless charging rates $\mathcal{C}_{r}=O(1)$ ($t_{D}\sim t_{C}$), the minimum flaw size activated $a_{0,\min}$ decreases as the inverse of the dimensionless charging rate squared *i.e., *$a_{0,\min}\sim\mathcal{C}_{r}^{-2}$. Thirdly, in Figs. 7–8, there exists a minimum flaw size $a_{0,\min}^{\infty}$ that flaws smaller, regardless of the dimensionless charging rate do not propagate.

The inverse square law can be understood by noting that the at moderate fluxes the concentration needs to penetrate at the particle length scale before there is enough energy for the flaw to activate (see Figs. 2–4). It is worth noting that the propagation at these critical fluxes is always abrupt only for smaller initial notches as detailed previously in this section (see Fig. 6). Therefore the time to activate an initial flaw is similar to the diffusion time of ions in particle size $t_{0}\sim t_{D}=R^{2}/D$. Using the mass conservation, we can write that the mass accumulated in the particle is equal to the mass of the ions inserted through its surface *i.e., *$R^{2}\Delta c\sim R\hat{J}\,t_{0}\sim R^{3}\hat{J}/D$. Rearranging the terms, we find the characteristic variation of Li ion’s concentration across the particle scales as $\Delta c\sim\hat{J}R/D$. This variation generates maximum hoop stress $\sigma_{\theta\theta}\sim\mathrm{E}\beta\Delta c/c_{\max}$ at the particle surface. According to the standard Griffith criterion, this stress can activate a flaw of size $a_{0,\min}\sim G_{c}\mathrm{E}/\sigma_{\theta\theta}^{2}$. Combing the above expressions for $\sigma_{\theta\theta}$, $\Delta c$, $t_{D}$, $t_{C}$, and using (24) we obtain the prediction

[TABLE]

In other words, the ratio of flaw size to the Griffith length-scale $a_{0,\min}/l_{G}$, scales as the inverse square of dimensionless charging rate times the maximum misfit strain *i.e., *$a_{0,\min}/l_{G}=A(\beta\mathcal{C}_{r})^{-2}$ where $A$ is a scaling constant. We can now identify the “misfit length scale”

[TABLE]

which takes into account that the magnitude of maximum misfit stresses generated $\mathrm{E}\beta$ is the appropriate measure of stresses in diffusion driven fracture. We should highlight that similar length scale was also used in [20] for the study of thermally driven cracks.

Using our phase-field simulations presented in Figs. 7–8 we can identify the scaling constants for the two particle sizes as $A=35.93,26.62$ for $R/l_{G}=4.2\times 10^{4},2\times 10^{5}$ particles respectively. Perhaps not surprisingly, since equation (25) has many simplifications and does not encode all particle size dependencies. Most notably it ignores the effect of the relative initial flaw size compared to the particle radius $a_{0}/R$, where changing the relative size of the initial flaw will change the evolution of the concentration field for the moderate charging rates considered. Furthermore, (25) also ignores the effects of the crack tip enrichment. As alluded to before, tensile stresses at the tip attract concentration, this introduces another length scale (see (37)) into the system that, in principle, can introduce dependency on the particle radius. Therefore, the scaling constant in the case of two different particle sizes are different. With the scaling constants extracted from the phase-field simulations we can carry the analysis further and obtain the maximum safe charging rate for a given particle size. For these practical charging rates, we can rewrite the maximum safe $C_{r,\max}$ below in which no flaws can be activated in a particle of radius $R$ as

[TABLE]

This scaling law predicts the most conservative charging rate (minimum charging time $t_{C}$) in terms of basic material properties.

Furthermore, to demonstrate the particle size dependency, it is easy to rearrange the power-law in equation (25) for a given dimensionless flaw size $\bar{a}_{0}=a_{0}/R$ as $R\sim\hat{J}^{-2/3}$. Fig. 9 depicts such a power-law emerging from the combined results of a series of simulations for particles of different size with a $a_{0}/R=0.1$ initial radial flaw on their surface. A similar power-law was independently derived by Bhandakkar and Gao [16] for initiation of a periodic array of cracks in a thin film using a cohesive zone model. There, authors study the initiation of an array of equidistant cracks such that the maximum stress in the system is equal to the cohesive strength of the material under study. They then, given the fracture energy of the material, investigate whether the displacement opening for the potentially initiated cracks will exceed the critical displacement required to maintain them. They find that regardless of the cohesive strength of the material, there exists a critical film thickness $H_{c}\sim\hat{J}^{-2/3}$ below which no fracture is initiated in the thin film.

Following our third observation, we see that the minimum flaw size activated $a_{0,\min}$ for small initial flaws deviates from the scaling law (25) and approaches a constant value $a_{0,\min}^{\infty}\simeq 10^{3}\,l_{G}$ for both particle sizes. At very large $\mathcal{C}_{r}$ required to activate these minimal flaws, the concentration reaches its minimum physically allowed value $c=0$ at the particle surface in a time $t_{0}\ll t_{D}$. In the next Section 4.2, we study the limit $\mathcal{C}_{r}\rightarrow\infty$ where $t_{0}/t_{D}\rightarrow 0$ by imposing the potentiostatic (Dirichlet) boundary condition $c=0$ at the particle surface.

4.2 Chemo-mechanical fracture of circular particles: (II) potentiostatic (Dirichlet) boundary condition

As discussed before for large fluxes $\mathcal{C}_{r}\rightarrow\infty$, the concentration field reaches its minimum $c=0$ at time $t_{0}\ll t_{C}$ as a result of which a depleted boundary layer is created on the surface of the particle. Therefore, it is more convenient to study this limit using Dirichlet boundary conditions. In this section, we present the results of numerical simulation for fracture of circular particles using $c(R)=0,~{}t\in[0,t_{\max}]$ Dirichlet boundary condition that is analogous to the maximum flux attainable for this system. Using this boundary condition, we can find the minimum flaw size activated for different particle radii. Similar to the previous section, we chose phase field length scale $\xi$ such that the initial flaw is well resolved (*i.e., *$a_{0}/\xi\geq 4$).

Fig. 10 shows the activation of a $a_{0}/l_{G}=400$ radial flaw in a $R/l_{G}=10^{5}$ particle under Dirichlet boundary conditions. Unlike the simulations analyzed in the previous section, the initial flaw is activated at $t_{a}\ll t_{D}$ in these simulations. We observe that the fracture propagation stems from the creation of an ion-depleted boundary layer of thickness $h\sim\sqrt{D_{0}t_{a}}$ with a size comparable to the minimum flaw size $a_{0}$ but much smaller than the particle radius *i.e., *$h\ll R$.

Fig. 11 shows the combined results of these numerical simulations for six particle sizes with different initial flaw sizes where we can make two observations. Firstly, our numerical results show that for our choice of parameters, there exists a maximum safe particle size $R_{\max}\simeq 10^{4}l_{G}$ that no flaw would propagate in it. Simply put, since the minimum activated flaw size decreases with the particle size, it becomes comparable to the particle size for small particles which cannot produce high enough deriving forces to propagate them. This size is analogous to critical thickness derived in [11] for a simple 1D bilayer. Secondly, we notice that as a function of growing particle radius $R\rightarrow\infty$, the smallest flaw activated asymptotically approaches a constant value $a_{0,\min}^{\infty}/l_{G}\rightarrow 5\times 10^{4}$. Thus the smallest flaw activated for a large particle becomes independent of its radius $R$. We can elucidate this observation noticing that in the absence of cracks, the maximum hoop stress generated under Dirichlet boundary conditions $\sigma\sim\mathrm{E}\beta$ is independent of the particle radius. Therefore, analogous to a flaw on the boundary of a half-space, the flaw size scales $a_{0,\min}^{\infty}\sim G_{c}\mathrm{E}/(\sigma)^{2}=l_{c}$ where the dimensionless prefactor is a function of the particle geometry. Also, we can relate our observations of minimum flaw size $a_{0,\min}$ at a given radius to results presented in the previous section 4.1 for large fluxes $\mathcal{C}_{r}\rightarrow\infty$. For example, for a $R/l_{G}=4.2\times 10^{4}$ particle presented in Fig. 7 the minimum flaw activated is predicted to be $a_{0}/l_{G}\simeq 4\times 10^{2}$ consistent with the results in Fig. 11.

4.3 Fracture of spherical cathodic particles with penny-shaped radial flaws

Although insight gained from the two-dimensional numerical simulations in the previous section is invaluable, only true 3D calculations can hope to capture all essential aspects of chemo-mechanical fracture in these particles. In this section, we demonstrate the similarities and differences between the simplified 2D and more realistic 3D simulations. Following the previous section, we model a spherical $R/l_{G}=2\times 10^{4}$ particle with a penny-shaped flaw on its surface (see 14). Unlike many similar calculations in 3D (*e.g., * [41, 40]), in this article, we simulate the complete sphere without explicit use of any symmetries. To this end, the elastic null-space (*i.e., *translation and rotational modes) was calculated and removed prior to the elastic sub-iteration in the alternate minimization algorithm (*i.e., *solving (21)). Since the computational cost of a uniform fine mesh was prohibitive, we chose a static adaptive meshing scheme, where for $r/l_{G}>1.6\times 10^{4}$ a fine mesh with an average edge length of $200l_{G}$ was generated and gradually coarsened to a coarse mesh with an average edge length of $500l_{G}$ for $r/l_{G}<1.2\times 10^{4}$. This meshing scheme, for different initial flaw sizes, resulted in the computational domain discretized into $\sim 13\textup{--}15\,\mathrm{M}$ tetrahedral elements (resulting overall in roughly the same number of degrees of freedom). Following the 2D simulations, we set the phase-field length scale to $\xi/R=3\times 10^{-2}$ and use the material properties as presented in table 1. Moreover, due to the prohibitive cost of the 3D simulations, we limit our investigation to three flaw sizes $a_{0}/l_{G}=10^{3},1.5\times 10^{3},2\times 10^{3}$ and charging rates $15\leq\mathcal{C}_{r}\leq 45$.

Fig. 12 shows the complex fracture topology that results from the activation of the penny-shaped flaw in 3D. The 3D fracture pattern highlights the role of dimensionality and follows from the fact that hoop stresses ($\sigma_{\theta\theta}$, $\sigma_{\phi\phi}$) reach their maximum values on the surface. Consequently, the tessellation of the particle surface by the crack releases the stresses and inhibits the inward propagation of the crack. These peripheral cracks only alleviate these stresses perpendicular to the crack surface, thereby causing new cracks to be initiated with different orientations than the plane of the initial penny-shaped crack. Therefore, we can hypothesize that for smaller charging rates where the opening stresses (Li ions) need to penetrate farther inside, the radial propagation is augmented compared to higher charging rates where the stresses generated are more superficial. This hypothesis is confirmed by the results of 3D simulations presented in Fig. 12. In that figure for all different initial flaw sizes simulated in 3D, crack propagation is abrupt and the added freedom for the cracks to release the stresses by tessellating the surface results in two dominant crack topologies: (I) cracks that propagate coplanar to the initial flaw under higher $\mathcal{C}_{r}$ (b,d in Figs. 12–17), and (II) cracks with initial coplanar propagation that tip split and result in a more complex topology under lower $\mathcal{C}_{r}$ (a,c,e,f in Figs. 12–17).

As explained before, we account for this transition using an argument similar to the one presented in Section 4.1 for abrupt versus continuous propagation in 2D. Unlike 2D radial cracks that can only release energy by penetrating towards the center of the particle, 3D penny-shaped cracks can both propagate radially and peripherally. The radial fracture in 3D simulations is akin to the radial propagation in 2D, *i.e., *the bulk elastic energy is released due to the crack opening up in the back. On the other hand, the peripheral propagation is analogous to the creation of (mod) cracks in a biaxially stretched thin-films [50, 12] or formation of imperfect polygonal patterns due to thermal quenching [20]. Therefore, since the highest opening stresses are always created on the surface of the particle, the initial propagation is always unstable in the peripheral direction. To clarify these two different fracture modes, we analyze two 3D topologies designated as cases a and b in Figs. 12–17. Due to the complex fracture topology in 3D, we use the dimensionless surface energy $\mathscr{F}_{\phi}/(G_{c}\,R^{2})$ as a measure of the surface area of the cracks, noticing that following equations (5) and (8):

[TABLE]

Fig. 13 depicts the dimensionless surface energy for $a_{0}/l_{G}=2\times 10^{3}$ at two different dimensionless charging rates: $\mathcal{C}_{r}=21$ (blue line in Fig. 13) and $\mathcal{C}_{r}=30$ (red line in Fig. 13). Fig. 14 shows the initial penny-shaped crack of size $a_{0}/l_{G}=2\times 10^{3}$ for cases a-b. Similar to the arguments presented for the 2D simulations at lower $\mathcal{C}_{r}$ (cases a,c,e,f in Figs. 12–17 and those depicted using red circles in Fig. 17), the flaw is only activated when the concentration has penetrated on the scale of the particle size. As seen, for example, in a-1 in Figs. 13 and 15, the propagation of the initial flaw is first planar which then tip splits due to high biaxial stresses (due to the symmetry of the problem far from the initial flaw $\sigma_{\theta\theta}=\sigma_{\phi\phi}$). At higher $\mathcal{C}_{r}$ (cases b,d in Figs. 12–17 and those shown using orange diamonds in Fig. 17) the initiation is faster and creates a coplanar crack with the initial flaw as seen, for example, in b-1 in Figs. 13 and 16. Upon further Li ion depletion, a secondary pair of cracks are initiated perpendicular to the initial circumferential crack as depicted in b-2–3 in Figs. 13 and 16. We should highlight that the radial penetration of the 3D penny-shaped crack is similar to radial propagation in 2D simulations. As a result, for case a at the lower $\mathcal{C}_{r}=21$ after the initial activation the crack abruptly penetrates radial distance of $\simeq 12\times 10^{3}\,l_{G}$ compared to $\simeq 6\times 10^{3}\,l_{G}$ for the case b at the higher $\mathcal{C}_{r}=30$.

Fig. 17 depicts the aggregate results of a series of 3D numerical simulations for a $R/l_{G}=2\times 10^{4}$ particle for $a_{0}/l_{G}=10^{3},1.5\times 10^{3},2\times 10^{3}$ initial penny-shaped radial flaws. Despite the major difference in crack propagation path (*i.e., *penetrating cracks in two-dimensional circular particles vs. surface cracks in the three-dimensional sphere), our 3D results suggest that critical flux to activate a surface flaw follows the inverse square power-law $a_{0,\min}\beta^{2}/l_{G}\sim\mathcal{C}_{r}^{-2}$ for moderate charging rates as in the 2D simulations. This is not surprising since the same arguments presented in Section 4.1 to justify the power-law still applies for spherical particles exposed to moderate fluxes. Furthermore, the results in Fig. 17 also suggest that, like 2D simulations (see Figs. 3–4), the transition of the inverse square power-law occurs at $a_{0}/l_{G}\simeq 10^{3}$.

We also should note that the tiling of the sphere surface is of particular theoretical interest. The polygonal tiling and its number of defects is prescribed by Euler’s celebrated theorem [51]. In contrast, in many physical systems, the number of defects on the curved surface goes beyond the minimum number necessary and is assigned by the local energetic minima. Over the past decade, significant progress has been made in closely connected areas of crystal formation on spherical surfaces [52, 53, 54] and pattern formation as the result of buckling [55]. Although the mechanism of surface tilings generated in this section is an attractive subject for further research, in this article, we limit ourselves to the general topology of the cracks generated.

5 Conclusions

In this article, we developed a thermodynamically consistent framework by combining the phase-field fracture method and diffusion to model chemo-mechanical fracture. We presented our formulation in Section 2 and detailed our implementation of it in Section 3.

As our first case study, we investigated in Section 4.1 the fracture of 2D circular disks. Using different initial flaw sizes, we showed how the steep gradient created as a result of the charging rate could cause a surface flaw to propagate and fracture the particle. Our numerical results show that for a given particle size, there exists a maximum flaw independent charging rate that can be used as a conservative limit in practice. Furthermore, motivated by our simulation results, we showed how the activation of the surface flaws follows an inverse square law $a_{0,\min}\sim\mathcal{C}_{r}^{-2}$ over intermediate dimensionless charging rates $\mathcal{C}_{r}=O(1)$ ($t_{D}\sim t_{C}$). We justified this power-law behavior based on a Griffith type analysis of the stresses generated far from the crack-tip and showed how it could be used to calculate a maximum safe charging rate $C_{r,\max}$ given the elastic and fracture properties as well as an estimate of the flaw sizes in the system. We should note that although the activation of surface flaws follows this simple power-law expression, the precise flux to activate a flaw is dictated by a non-trivial concentration profile around the crack-tip. Since the scaling law analysis ignores the ratio of flaw size to the radius of the particle, as well as the crack-tip enrichment, the scaling constant can only be derived from the numerical simulations, especially for small particles where the initial flaw size plays a more significant role. Furthermore, we showed that depending on the particle and flaw size, the initial propagation could be abrupt or continuous for low and high fluxes, respectively. While puzzling at first, we described how the abrupt fracture propagation length decreases for increasing fluxes for moderate initial flaw sizes due to smaller bulk energy available at the fracture onset.

We then extended our study to high fluxes that are necessary for the activation of very small flaws ($a_{0}<10^{3}\,l_{G}$ for our choice of parameters). Our results show that for these small flaws, the safe charging rate deviates from the previously obtained scaling law. We explained our observation, noting that the high charging rate creates a depleted layer on the periphery of the particle and thus loses its effectiveness in creating a steep enough gradient to activate these flaws. As a result, we found out that there exists a minimum safe flaw size $a_{0,\min}$ for a given particle size that does not propagate under any charging rate. To effectively address these maximal charging rates, in Section 4.2, we examined the activation of surface flaws using potentiostatic (Dirichlet) boundary conditions. Our simulation results show that there exists a C-rate independent, safe particle size that no flaw of any size will propagate in it. In addition, they show that in large particles the minimum activated flaw size approaches a constant value (*e.g., *$a_{0,\min}^{\infty}\simeq 200l_{G}$ for our choice of parameters). In other words, our numerical simulations suggest that particles (no matter how large) containing flaws smaller than $a_{0,\min}^{\infty}$ do not crack due to diffusion-driven misfit stresses.

Finally, in Section 4.3, to investigate the role of dimensionality, we performed a series of 3D simulations on spherical particles with penny-shaped flaws. Using our numerical observations, we showed that, unlike in 2D and assumptions [17] and results [40] of previous studies, the crack topology changes from a coplanar penetrating mode to a surface tiling mode. These full (*i.e., *without any symmetries assumed) 3D calculations show that the initial mechanical mode of failure in three-dimensional particles during charging is due to the fracture on their surface. While all the propagations from the initial penny-shaped crack in 3D were abrupt, we showed how the change of the fracture topology could be explained using arguments akin to those used to justify the length of abrupt propagation in 2D. Furthermore, our admittingly limited 3D results suggest that $a_{0,\min}\sim\mathcal{C}_{r}^{-2}$ scaling law is still valid in 3D for flaws larger than $a_{0}>10^{3}\,l_{G}$.

Lastly, it is crucial to highlight that, in this article, we only model chemo-mechanical fracture due to Lithium diffusion with no phase change or discontinuity in expansion. As highlighted, for example, in [18, 19], coherency stresses generated at the phase and grain boundaries can result in charging rate independent fracture in Li-storage materials.

6 Acknowledgments

Acknowledgments: A.M. and A.K. acknowledge the support of Grant No. DE-FG02-07ER46400 from the U.S. Department of Energy, Office of Basic Energy Sciences. The majority of the numerical simulations were performed using resources of the Extreme Science and Engineering Discovery Environment (XSEDE) under the resource allocation TG-MSS160013. Additional numerical simulations were also performed on the Northeastern University Discovery cluster at the Massachusetts Green High Performance Computing Center (MGHPCC).

Appendix A Concentration enrichment around crack-tip due to mechanical loads

Examining equations (13)-(14), it is easy to notice that the ions flow toward the regions with higher hydrostatic pressures; therefore, it is not surprising that in the presence of a crack, a higher concentration will accumulate at the crack-tip (see Fig. 18 for example). More specifically, for small eigen-strains (*i.e., *$\beta\ll 1$) the stresses become independent of the concentration field (*i.e., *the diffusion equation would be driven by the magnitude of hydrostatic stress).

To derive the enrichment at the crack-tip, we can rewrite the coupled equations of elasticity and concentration in terms of Airy stress function $\mathcal{A}$ in 2-D:

[TABLE]

where $\mathrm{E}^{*}=\dfrac{\mathrm{E}}{1-\nu^{2}}$ for plane-stress. Therefore, for small eigen-strains *i.e., *$\beta\ll 1$ the stresses become independent of the concentration field (*i.e., *the diffusion equation would be driven by the magnitude of hydrostatic stress). Thus, for steady-state conditions and in absence of a surface flux, one can write

[TABLE]

A similar solution to (31) can be obtained for the dilute approximation where $f(c)=c\ln(c)-c$ as:

[TABLE]

The above expression was also derived by the direct solution of the diffusion equation in [10].

To find the concentration profile around the crack-tip, we can replace the expression of $\mathrm{Tr}(\sigma)$ from the asymptotic plane-stress solution of mode-I fracture:

[TABLE]

[TABLE]

where $K_{I}$ is the stress intensity factor. After some algebra $\mathrm{Tr}(\sigma)$ can be written as:

[TABLE]

where $K_{I}$ is the mode-I stress intensity factor. Using (35) we can write the concentration around the crack-tip at the time of fracture as:

[TABLE]

where

[TABLE]

can be identified as the intrinsic length scale for the concentration of ions around the crack-tip. Equation (37) shows that the ratio of the enrichment length scale to Griffith length scale scales as the square ratio of maximum misfit stresses to chemical energy. In (36) one can find the steady-state chemical potential $\mu_{0}$, from far field concentration as

[TABLE]

As we showed in the Section 4.1, crack-tip enrichment is a common occurrence in diffusion-driven fracture of Lithium-ion battery particles. We can easily calculate the length scale $r_{c}/l_{G}\simeq 7591.75$ for \ceLiMn2O4 at room temperature where the crack-tip concentration is captured approximately by (36). Fig. 19 shows a comparison between the results of the numerical simulation for a $R/l_{G}=4.2\times 10^{4}$ particle (Fig. 2) and (36). The simulation is performed in a circular geometry of radius $R$ containing a sharp $1^{\circ}$ notch from $r=-150\,\xi$ to $r=0$ at $\theta=\pi$. To simulate near tip stress fields, the displacement fields associated with (33)–(34) were imposed on the boundary of the domain. The concentration is initially uniform $c/c_{\max}=0.5$ everywhere and the value of $\mu_{0}$ was calculated based on the resulting concentration at $t/t_{D}=1$ and $r=R$.

We should note that the radial crack, driven by charging the cathodic particle, can stop propagating in the middle of the particle. In this situation the enrichment carried by the crack-tip can be shielded from the depleting flux by chemo-mechanical force exerted at the crack-tip. The remaining concentration then can change the dynamics of the charging process. Furthermore, while the main focus of this article is on the diffusion of Li-ions in battery particles, crack-tip enrichment can play an important role in other systems where diffusion and fracture happen concurrently such as corrosive cracks, crack-tip embrittlement, and fracture in poroelastic media [56, 57].