# Capturing the stress evolution in electrode materials that undergo phase transformations during electrochemical cycling

Bo Wang, Julien Réthoré, Katerina E Aifantis

► To cite this version:

Bo Wang, Julien Réthoré, Katerina E Aifantis. Capturing the stress evolution in electrode materials that undergo phase transformations during electrochemical cycling. International Journal of Solids and Structures, 2021, 10.1016/j.ijsolstr.2021.03.019. hal-03183921

HAL Id: hal-03183921

https://hal.science/hal-03183921

Submitted on 29 Mar 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Capturing the stress evolution in electrode materials
that undergo phase transformations during
electrochemical cycling

Bo Wang<sup>a</sup>, Julien Réthoré<sup>b</sup>, Katerina E. Aifantis<sup>a,*</sup>

<sup>a</sup>Department of Mechanical and Aerospace Engineering, University of Florida, Gainesville, Florida, 32608, USA

<sup>b</sup>Research Institute in Civil and Mechanical Engineering (Gem), CNRS UMR 6183 CNRS / Ecole Centrale de Nantes / Université de Nantes, F-44 321, Nantes, France.

**Abstract**

The present work sheds light on the stresses generated in a spherical particle subjected to phase transformations during ion-insertion. In order to account for the physical process that occur during electrochemical cycling, the models used are those of small deformation and account for the effects of phase trans- formation, chemo-mechanical coupling and concentration-dependent material properties. The two-phase lithiation is modeled by the Cahn-Hilliard equation. It is found that the DISs arise from the inhomogeneous volume expansions re- sulting from Li concentration gradients and the hydrostatic stress facilitates the diffusion of Li-ions under elastic deformation while it hinders diffusion in the plastic case. When the elastic modulus is reduced the magnitude of the diffusion-induced stress decreases but the strain increases under elastic defor- mation whereas the opposite occurs for the plastic case. Furthermore, if the electrode is assumed to undergo strain softening during plastic deformation, smaller stresses and higher plastic strains are predicted than when strain hard- ening is assumed. The models highlight the importance of chemo-mechanical coupling effects, concentration-dependent material properties and plastic defor- mation on diffusion-induced stresses. These findings render prospective insights

Email address: kaifantis@ufl.edu (Katerina E. Aifantis)

<footnote>*Corresponding author</footnote>
for designing next-generation mechanically stable phase transforming electrode materials.

*Keywords:* Electrodes, Phase transformation, Chemo-mechanical coupling, Concentration-dependent elastic modulus

# 1. Introduction

Lithium ion batteries (LIBs) are currently widely used in portable electronics and will continue to play a significant role in next generation electric vehicles, utility grids and electric/hybrid-electric airplanes. To fulfill the practical requirements, LIBs should have an optimal electrochemical and mechanical performance, i.e, high energy density, high power density, and long cycle life. The operation of LIBs involves the continuous diffusion of Li-ions between the anode and cathode and therefore the performance of LIBs depends largely on the electrode active materials (materials that react with Li⁺).

The cathode is the source of lithium ions (Li⁺), and the anode should have the ability to host the Li⁺ during charge and vice versa for the discharge process. Lithium insertion and de-insertion are accompanied by changes in lattice spacing. As a result, the lattice structure can be altered, and phase transformations can occur, resulting in volume changes. Consequently, diffusion-induced stresses (DISs) are generated. The energy density, power density and cycle life of LIBs are closely related to DISs. Excessive stresses result in pulverization of the electrode materials and their detachment from the current collector. As a result, fracture of the electrode and capacity decay occurs.

For instance, commercially used cathode materials (LiMn₂O₄ and LiCoO₂) for LIBs have achieved only half of the theoretical capacity, due to the mechanical stresses that arise from phase transformations and changes in the crystal structures/lattice parameters during the charging/discharging processLiu et al. (2004); Venkatraman et al. (2000); Amatucci et al. (1996); Aifantis & Hackney (2003). Similarly, high-capacity anode materials (Si and Sn) cannot retain their initial capacity (900-4000 mAh/g) due to the fracture that results from
the lithiation-induced volume expansions (300%-400%) He et al. (2012); Wang et al. (2013). Therefore, understanding the stress evolution process is crucial for revealing the degradation mechanisms of LIBs and thus, improving their performance. Initial models considered the stress evolution Aifantis & Hackney (2003) and fracture Aifantis et al. (2005, 2006) in active electrode particles by using a purely mechanical approach. Since then, multiple models have been developed to determine the DIS in material systems. The phenomenon of DIS was first studied by employing an analogy between thermal stress and DIS to analyze the transverse stresses introduced by solute lattice contraction of boron and phosphorus in a thin silicon plate during mass transport Prussin (1961). A number of analytical solutions to DIS problems in elastic media of simple spherical and cylindrical geometries were provided in Li (1978), while the DISs arising from inhomogeneous concentration in materials was in investigated in Larcht'e & Cahn (1982); Larche & Cahn (1984). The effect of DISs on diffusion in hollow cylinders Lee et al. (2000); Wang et al. (2002) and square sandwich composites Ko et al. (2005) has also been studied. Later on, numerous researchers followed this analogy approach to investigate the DISs in nanowire electrodes Deshpande et al. (2010a), spherical electrode particles Verbrugge & Cheng (2009); Cheng & Verbrugge (2010, 2008), cylindrical electrodes Deshpande et al. (2010b); Song et al. (2012) and layered electrode plates Zhang et al. (2012).

For intercalation-induced stresses, one-dimensional (1D) models have been developed by to investigate the stress generation due to Li⁺ insertion into a spherical carbon anode Christensen & Newman (2006b), but instead of thermal analogy, the governing equations were derived based on lattice deformation analysis. Given that analytical solutions to diffusion equations under specific conditions are available Crank (1979); Jaeger & Carslaw (1959), analytical expressions have been obtained to study the stress evolution in spherical electrode particles Cheng & Verbrugge (2009, 2010).

Most of these studies are centered on single-phase reactions, assuming that the electrode particle phase does not change during Li⁺ insertion/de-insertion. However, it's well known that lithium-ion insertion often leads to the forma-
tion of different phases. Phase transformations have been experimentally ob-
served during Li⁺ insertion/de-insertion and could play an important role in the
stress evolution during the discharging/charging process Shin & Pyun (1999a,b);
Zhang & White (2007). Multiple models have been constructed to capture the
effect of phase transformations on DISs. The “shrinking core” model Srinivasan
& Newman (2004) considers that the core and shell would consist by lithium-rich
and lithium-poor phases, respectively, during Li⁺ insertion and vice-versa for
de-insertion. This model leads to an abrupt change (jump) in the concentration
at the phase boundary which also results in a stress discontinuity at the inter-
face between the two phases. By using this model and extending the theoretical
framework for a single-phase reaction Christensen & Newman (2006b), the stress
generation in a spherical LiMn₂O₄ particle was estimated Christensen & New-
man (2006a). This “shrinking core” model has also been applied to study other
spherical electrode particles and other configurations Deshpande et al. (2011);
Liu et al. (2014); Esmizadeh et al. (2019). Due to the sharp interface, concen-
tration and stress discontinuities appear at the interface. It is noted that the
moving phase boundary needs to be tracked and specific boundary conditions
need to be applied during insertion and de-insertion. Therefore, this sharp in-
terface model poses numerical difficulties for practical applications. A possible
alternative approach would be phase field modeling of the diffusion. By em-
ploying a diffusive interface model, the ion distributions can be modeled with a
continuous concentration field. As a result, the reaction front does not need to
be explicitly tracked. Such phase field models have been successfully applied to
study the stress generation of LiFePO₄ Zhang & Kamlah (2018) and LiMn₂O₄
Huttin & Kamlah (2012) particles.

Phase field models, aside from the work in Hu et al. (2019a) focus on the
DISs of cathode materials. Due to the generally small volume expansions that
cathodes experience during lithiation, infinitesimal/small strain continuum me-
chanics formulations can be used. For example, LiFePO₄ shrinks approximately
7% after full delithiation Yamada et al. (2005); Wang et al. (2005), and therefore
cathodes are assumed to undergo linear elastic deformation during the interca-
lation/deintercalation process Woodford et al. (2010); Hao & Fang (2013a,b). However, all high capacity anodes, such as Si, are well known for large volume expansions during lithiation, which can reach 420% upon maximum lithiation (Si₄.₄Li). The deformation is so large that plastic deformation can take place. Although linear elastic studies of the DISs Jagannathan & Chandran (2014); Tsagrakis & Aifantis (2018); Golmon et al. (2010) in these anodes can provide some insights on the stress evolution, plastic models Hu et al. (2010); Zhao et al. (2011b) would be more reasonable since the volume change is large enough that yielding can occur. Both in-situ measurements Sethuraman et al. (2010) and first-principle studies Zhao et al. (2011a) of the stress evolution in thin-film amorphous Si anodes during lithiation/delithiation have revealed that yielding and plastic flow occur. Elastoplastic and viscoplastic models have been adopted to study the DISs in Si and Sn anodes during lithiation Chen et al. (2017); Huang et al. (2013); Liu et al. (2012); Chen et al. (2014). An important finding in these aforementioned studies is that the tangential stress transits from compressive to tensile during lithiation when plasticity is taken into account, which is believed to be the cause of the experimentally observed fracture of Sn nanoparticles Aifantis et al. (2012), Si spherical particles Liu et al. (2012) and Si nanowires Ryu et al. (2011) during lithiation. Elastic models predict a compressive hoop stress on the outer surface at all times while plastic models can capture this transition (compression in the core to tension in the shell) since the hoop tension is the crack driving force of the observed surface cracking Huang et al. (2013); Liu et al. (2012); Chen et al. (2014); Chang et al. (2018).

Despite the fact that there have been a lot of models on DISs of electrode materials, a comprehensive understanding has yet to be achieved. One aspect is the chemo-mechanical (CM) coupling effect. Initial works, particularly those using the “shrinking core” model, did not consider CM coupling. Subsequent models considered the coupling effect between chemical diffusion and mechanical deformation Christensen & Newman (2006b); Liu et al. (2014); Zhang & Kamlah (2018); Hao & Fang (2013a,b); Tsagrakis & Aifantis (2018); Zhao et al. (2011a); Chen et al. (2014). However, few works have taken into considera-
tion the change that occurs in the elastic modulus of electrode materials during Li⁺ insertion. A decrease or increase in the elastic modulus during Li-insertion (lithiation), as a function of Li-concentration, have been observed experimentally Sethuraman et al. (2010); Hertzberg et al. (2011) and calculated through first principle calculations Maxisch & Ceder (2006); Qi et al. (2014, 2010). Table 1 shows the Young’s modulus of the lithiated and delithiated phases of some electrode materials. The effect of a concentration dependent elastic modulus on DIS was systematically examined using linear and gradient elasticity Natarajan & Aifantis (2020), showing that there was a negligible difference in DISs when the variation of the elastic modulus due to Li insertion was only 10% and a single-phase reaction during Li-insertion occurred. In such cases, constant moduli may be used for modeling DISs in a single-phase battery system. However, very few studies examined the effects of the concentration-dependent material properties on DISs in phase-transformation electrode materials, and there isn’t much work comparing plastic and elastic models. Herein we present a new formulation that can account for chemomechanical coupling, concentration dependent elastic moduli for two phase transitions and large volume expansions, for both electrodes that undergo either pure elasticity or also plasticity during lithiation. By performing such a comprehensive study of the DIS, more accurate predictions of the stress fields in the electrodes can be provided, allowing for a better selection of the materials and microstructure that will result in improved electrochemical performance. It should be noted that this work can be applied to other electrochemical systems that operate under the diffusion of ions, such as Na-ion batteries.

## 2. Phase field modeling of DISs

As mentioned earlier, a core-shell model is often used to model the diffusion process of electrode materials that exhibit phase segregation into Li-rich and Li-poor phases, respectively, during Li-ion insertion and de-insertion. It’s often assumed that during insertion the phase boundary moves in the radial direction
<table><caption>Table 1: Young's modulus of some active electrode materials in lithiated and delithiated phases</caption><thead><tr><th>Delithiated phase</th><th>Method</th><th>Young's modulus(GPa)</th><th>Lithiated phase</th><th>Young's modulus(GPa)</th></tr></thead><tbody><tr><td>CoO<sub>2</sub></td><td>Computation<br>Qi et al.<br>(2014)</td><td>59.8</td><td>LiCoO<sub>2</sub></td><td>264</td></tr><tr><td>C(graphite)</td><td>Computation<br>Qi et al.<br>(2010)</td><td>32</td><td>LiC<sub>6</sub></td><td>109</td></tr><tr><td>Si(amorphous)</td><td>Computation<br>Shenoy et al.<br>(2010)</td><td>96</td><td>Li<sub>3.75</sub>Si</td><td>41</td></tr><tr><td>Si(crystaline)</td><td>Experiment<br>Hertzberg<br>et al. (2011)</td><td>90</td><td>Li<sub>3.75</sub>Si</td><td>12</td></tr><tr><td>Sn(crystaline)</td><td>Computation<br>Stournara<br>et al. (2012)</td><td>51</td><td>Li<sub>3.5</sub>Sn</td><td>24.7</td></tr></tbody></table>

and the expansion is isotropic. To model the sharp phase boundary, the Cahn-Hilliard type diffusion equation can be used. The effect of chemo-mechanical coupling on DISs can be incorporated into the model via a variational formulation approach. As for the electric field, the conductivities of the electrode materials are high enough in most cases so that the electric potential inside the electrode particles are regarded as uniform, therefore, the influence of electric field on DIS can be neglected Zuev & Tsvetkov (2017).

Both, the concentration and temperature fields are governed by the diffusion equation. Therefore, it's naturally to deduce that a change in concentration will cause DIS in a similar manner as a change in temperature would lead to
thermal stresses, as has been done in most of the aforementioned articles, and the constitutive relations from thermo-elasticity can be adopted without too much adjustment. For plastic models, either viscoplastic or elastoplastic models can be used. Viscoplastic models for thermal-induced stress have been implemented into commercial finite element software such as ABAQUS and COMSOL. When using such software, the DIS due to (de-)insertion can be modeled easily just by replacing the temperature with the concentration and the thermal expansion coefficient with the chemical expansion coefficient. For implementation in open-source finite element codes such as deal.II and FEniCS, elastoplastic models are much more convenient and are therefore implemented in FEniCS in this work.

## 2.1. DISs for elastic solids

### 2.1.1. Phase field model for elastic solids

Based on a variational formulation, a theoretical framework for DIS can be established. The basic idea is energy minimization of the total potential energy. For the cathode with small volume variations, the deformation is assumed to remain in the elastic stage. The system (total) free energy of the domain $\Omega$ under consideration is

$$ \Psi = \int_{\Omega} \psi d\Omega \qquad (1) $$

where $\Psi$ is the free energy density, i.e., free energy per unit reference volume and is comprised of two parts:

$$ \psi = \psi(\varepsilon, c, \nabla c) = \psi^e + \psi^c \qquad (2) $$

where $\psi^e$ is the bulk elastic energy density solution, $c$ is the actual species concentration, given in terms of moles per unit reference volume (the molar concentration per unit volume in the reference configuration), and $\psi^c$ is the energy contribution due to changes in the concentration solution. The strain and displacement are related as:

$$ \varepsilon = \frac{1}{2} [\nabla u + (\nabla u)^T] \qquad (3) $$
where **u** is the displacement field. The elastic strain energy density is

$$ \psi^e (\varepsilon^e) = \frac{1}{2} \varepsilon^e : \mathcal{C} : \varepsilon^e \qquad (4) $$

where the elastic strain tensor is $\varepsilon^e = \varepsilon - \varepsilon^c$. $\varepsilon^c$ is the chemically induced strain tensor and is given by:

$$ \varepsilon^c = \alpha (c^* - c_0^*) \mathbf{I} \qquad (5) $$

where **I** is a 2<sup>nd</sup>-order identity tensor, $\mathbb{C}$ is a 4<sup>th</sup>-order elasticity tensor, given by $\mathbb{C} = \lambda\delta_{ij}\delta_{kl} + \Lambda(\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk})$, $\delta_{ij}$ describes the Kronecker delta, $\lambda$ and $\Lambda$ are the Lamé constants, and $c^*$ is the relative concentration. The normalized concentration $c^*$ is used to characterize the relative saturation level of Li in an electrode. It's normalized by the maximum concentration $c_{max}$ in the Li/electrode system, ranging from 0 (no ions) to 1 (full ion insertion), given by: $c^* = c/c_{max}$. $c_0^*$ is a reference concentration, usually set as 0. $\alpha$ is a dilatation coefficient (analogous to the thermal expansion coefficient), given by:

$$ \alpha = \sqrt[3]{1 + \Xi c_{\max}} - 1 \doteq \frac{1}{3} \Xi c_{\max} \qquad (6) $$

where $\Xi$ is the partial molar volume of Li-ions. Due to the presence of a sharp interface, the energy contribution due to a concentration gradient on diffuse interfaces cannot be neglected. It should be noted that for phase transformation materials, there is a sharp interface, however, phase field modeling uses diffusive interfaces with a finite thickness (but very thin) to approximate it, resulting in a continuous concentration. In traditional sharp interface models, the thickness is assumed to be zero, and as a result, the concentration is discontinuous across the interface. The chemical energy density can be split into two parts:

$$ \psi^c = \psi^c(c, \nabla c) = \psi^{ch} + \psi^{cg} \qquad (7) $$

where $\psi^{ch}$ is the homogeneous chemical free energy density (free energy density of a homogeneous system of uniform concentration) and $\psi^{cg}$ is the chemical gradient energy density due to the concentration gradient $\nabla c$. The chemical stored free energy (chemical contribution to the stored energy) is assumed to
have the form:

$$ \psi^{ch} = \psi^{ch}(c) = RTc_{\max} [c^* \ln c^* + (1-c^*) \ln(1-c^*) + \chi c^* (1-c^*)] + \mu_0 c_{\max} c^* \quad (8) $$

Here R is the gas constant, T is the absolute temperature, and RT denotes the specific molar energy. $\mu_0$ is a reference value of the chemical potential of the diffusing species, set as 0 since it will not affect the diffusion behavior. $\chi$ is the constant partial molar volume (the volume of a mole of the species, parameter $\chi$ indicates the convexity of the energy) which is dimensionless and controls the shape of the double-well energy function characterizing the Li-rich and Li-poor phases. This form of the homogenous energy density favors a separation of the system towards $c^* = 0$ and $c^* = 1$, i.e., minimization of the Gibbs free energy. To ensure phase separation, $\chi > 2$. The chemical gradient energy density is given by:

$$ \psi^{cg} = \frac{\kappa}{2} |\nabla c|^2 = \frac{\bar{\kappa}}{2} |\nabla c^*|^2 \quad (9) $$

depending on the concentration gradient and representing an interfacial energy. $\bar{\kappa} = \kappa c_{\max}^2$ and $\kappa$ is the gradient energy coefficient, with units of energy per unit volume times a length squared. The parameter $\kappa$ modulates the contribution of the large concentration gradient at the two-phase interface to the free energy. When the elastic property of the material is strongly influenced by concentration, based on rule of mixtures from composites, it can be approximated by a simple linear relation:

$$ E(c^*) = E_0 + kc^* \quad (10) $$

where E is the concentration-dependent Young's modulus, and $E_0$ is the Young's modulus of the delithiated phase (i.e. initial elastic modulus), and k is a positive or negative parameter allowing for a reduction or increase in the modulus during lithiation. The chemical potential $\mu$ is defined as the variational derivative of
the system free energy:

$$
\begin{align}
\mu &= \frac{\delta\psi}{\delta c} = \frac{\partial\psi}{\partial c} - \nabla \cdot \left(\frac{\partial\psi}{\partial\nabla c}\right) \nonumber \\
&= \mu_0 + RT \ln \left(\frac{c^*}{1-c^*}\right) + RT\chi(1-2c^*) - \Xi\sigma_h \tag{11} \\
&\quad - c_{\max}\nabla \cdot (\kappa\nabla c^*) - \frac{1}{2c_{\max}} \frac{dS_{ijkl}}{dc^*} \sigma_{ij}\sigma_{kl} \nonumber
\end{align}
$$

where $\sigma_h = \frac{1}{3}\text{tr}(\boldsymbol{\sigma})$ is the hydrostatic stress and $S_{ijkl}$ is the concentration-dependent compliance tensor. Therefore, the diffusional potential is related to the Li concentration and the local stress state, as well as the concentration gradient and the changing material properties due to intercalation. For linear elastic and isotropic materials, the compliance tensor is given by:

$$S_{ijkl} = \frac{1-2v}{3E} \delta_{ij} \delta_{kl} + \frac{1+v}{2E} \left( \delta_{ik} \delta_{jl} + \delta_{il} \delta_{jk} - \frac{2}{3} \delta_{ij} \delta_{kl} \right) \quad (12)$$

where $v$ is Poisson's ratio. As a result,

$$\frac{dS_{ijkl}}{dc^*} \sigma_{ij} \sigma_{kl} = \frac{1}{E^2} \frac{dE}{dc^*} (v\sigma_{ii}\sigma_{kk} - (1+v)\sigma_{jl}\sigma_{jl}) \quad (13)$$

The Li transport in the host materials can be described by the diffusion equation, as

$$\frac{\partial c^*}{\partial t} + \nabla \cdot \mathbf{j} = 0 \quad \text{in } \Omega \quad (14)$$

where $\mathbf{j}$ is the flux, and $\Omega$ is the bulk of the electrode. The gradient of the chemical potential is the chemical driving force for diffusion, i.e. the gradient of the diffusional potential of Li $\mu$ serves as the driving force for Li diffusion:

$$\mathbf{j} = -M\nabla\mu \quad (15)$$

where $M$ is the mobility tensor, which is generally a function of Li concentration $c$. For isotropic materials, it's assumed to have the form:

$$M = c^* (1 - c^*) M I \quad (16)$$

where $M$ is the solute mobility and is independent on time. The diffusion coefficient $D$ is given by:

$$D = MRT \quad (17)$$
To solve for the diffusion, suppose the electrode is under galvanostatic operation and the particle is charged galvanostatically, with a uniform lithium flux through its free surface according to

$$Mc^*(1-c^*) \frac{\partial \mu}{\partial r} \Big|_{r=\rho} = i \qquad (18)$$

Here $i$ is the applied current density (per unit area), and $\rho$ is the radius of the electrode particle. The variational boundary condition is given by:

$$\kappa \nabla c \cdot \mathbf{n}_c = 0 \quad \text{on} \quad \partial \Omega^c \qquad (19)$$

with $\mathbf{n}_c$ being the outward unit normal to the concentration boundary $\partial \Omega^c$. The weak form of the diffusion equation can be given by:

$$\int_{\Omega} \frac{\partial c^*}{\partial t} v dV = \int_{\partial \Omega} iv dS - \int_{\Omega} Mc^*(1-c^*) \nabla \mu \nabla v dV \qquad (20)$$

$$\int_{\Omega} \mu w dV = \int_{\Omega} R T \left[ \ln \left( \frac{c^*}{1-c^*} \right) + 2\chi (1-c^*) \right] w dV - \int_{\Omega} \kappa c_{\max} \nabla c^* \nabla w dV \\ - \int_{\Omega} \Xi \sigma_h w dV - \int_{\Omega} \frac{1}{2c_{\max}} \frac{dS_{ijkl}}{dc^*} \sigma_{ij} \sigma_{kl} w dV \qquad (21)$$

where $v$ and $w$ are arbitrary test functions. Considering a spherical particle during the lithiation process, the mechanical equilibrium equation under spherically symmetrical deformation is given by:

$$\frac{d\sigma_r}{dr} + \frac{2(\sigma_r - \sigma_\theta)}{r} = 0 \qquad (22)$$

where $\sigma_r$ is the radial stress, $\sigma_\theta$ is the hoop stress and $r$ is the radial distance. The constitutive law (stress-strain relation) for an isotropic spherical electrode during lithiation is given by:

$$\left. \begin{aligned} \sigma_r &= \frac{E}{(1+\nu)(1-2\nu)} \left[ (1-\nu) \left( \frac{du}{dr} - \alpha c^* \right) + 2\nu \left( \frac{u}{r} - \alpha c^* \right) \right] \\ \sigma_\theta &= \frac{E}{(1+\nu)(1-2\nu)} \left[ \left( \frac{u}{r} - \alpha c^* \right) + \nu \left( \frac{du}{dr} - \alpha c^* \right) \right] \end{aligned} \right\} \qquad (23)$$
Here $\epsilon_r = \frac{du}{dr}$ is the radial strain and $\epsilon_\theta = \frac{u}{r}$ is the tangential strain and $u$ is the radial displacement. The boundary conditions for this problem are: zero radial displacement at the center and zero radial stress on the outer surface, i.e., $u|_{r=0} = 0$ and $\sigma_r|_{r=\rho} = 0$. As a result, when $E$ is not dependent upon $c$, then the stress field can be given analytically by Hu et al. (2019b):

$$ \left. \begin{aligned} \sigma_r &= \frac{2\alpha E}{1-\nu} \left( \frac{1}{\rho^3} \int_0^\rho c^* r^2 dr - \frac{1}{r^3} \int_0^r c^* r^2 dr \right) \\ \sigma_\theta &= \frac{2\alpha E}{1-\nu} \left( \frac{1}{\rho^3} \int_0^\rho c^* r^2 dr + \frac{1}{2r^3} \int_0^r c^* r^2 dr - \frac{c^*}{2} \right) \end{aligned} \right\} \quad (24) $$

For a general case, the Young's modulus may not be a constant. In such cases, the analytical solution cannot be obtained, and finite element or other numerical methods need to be employed. The weak form of the mechanical equilibrium equations can be given by:

$$ \int_0^\rho q \left( \frac{d\sigma_r}{dr} + \frac{2(\sigma_r - \sigma_\theta)}{r} \right) r^2 dr = 0 \quad (25) $$

where $q$ is a test function. Integrating by parts and applying boundary conditions, one can finally obtain:

$$ \int_0^\rho \left( \sigma_r \frac{dq}{dr} + 2\sigma_\theta \frac{q}{r} \right) r^2 dr = 0 \quad (26) $$

By substituting Equation 23 into Equation 26, we have the weak form in terms of the displacement as follows:

$$ \int_0^\rho q \left[ (\lambda + 2\Lambda) r^2 \frac{d^2 u}{dr^2} + 2(\lambda + \Lambda) \left( r \frac{du}{dr} - u \right) - (3\lambda + 2\Lambda) \alpha r^2 \frac{dc^*}{dr} \right] dr = 0 \quad (27) $$

### 2.1.2. Numerical implementation of the phase field model

In order to seek the solutions, the finite element method is used. Due to spherical symmetry, the problem can be reduced to a one-dimensional problem. The weak form cannot be solved directly. Instead, we have to discretize in both the temporal and spacial space. Therefore, Equation 20 and Equation 21 can be rewritten as:

$$ \int_{\Omega} \frac{c_{n+1}^* - c_n^*}{\Delta t_n} v dr = - \int_{\partial \Omega} iv dS - \int_{\Omega} M c_n^* (1 - c_n^*) \nabla \mu_n \nabla v dr \quad (28) $$
$$ \int_{\Omega} \mu_{n+1} w dV = \int_{\Omega} R T \left[ \ln \left( \frac{c_n^*}{1 - c_n^*} \right) + 2\chi (1 - c_n^*) \right] w dV - \int_{\Omega} \kappa c_{\max} \nabla c_n^* \nabla w dV \\ - \int_{\Omega} \Xi (\sigma_h)_n w dV - \int_{\Omega} \frac{1}{2c_{\max}} \frac{dS_{ijkl}}{dc^*} (\sigma_{ij})_n (\sigma_{kl})_n w dV \qquad (29) $$

where the subscripts n + 1 and n are the corresponding values at the time step t$_{n+1}$ and t$_{n}$ respectively, with Δt$_{n}$ = t$_{n+1}$ − t$_{n}$ being the n$^{th}$ time step. The hydrostatic stress at t$_{n}$ is given by:

$$ (\sigma_h)_n = \frac{1}{3} [(\sigma_r)_n + 2(\sigma_\theta)_n] \qquad (30) $$

And Equation 13 at t$_{n}$ can be given by:

$$ \frac{dS_{ijkl}}{dc^*} (\sigma_{ij})_n (\sigma_{kl})_n = \frac{1}{E^2} \frac{dE}{dc^*} \left[ [(\sigma_r)_n + 2(\sigma_\theta)_n]^2 - [(\sigma_r)_n^2 + 2(\sigma_\theta)_n^2] \right] \quad (31) $$

By using this kind of mixed formulation, the problem becomes to find c$_{n+1}$* and μ$_{n+1}$ given that c$_{n}$* and μ$_{n}$ are known. The built-in DOLFIN Newton solver in FEniCS Dol is used to solve the nonlinear equations. As for the mechanical equilibrium equation, the built-in NonlinearVariationalSolver in FEniCS is employed. The basic idea is to calculate the Gateaux derivative of the nonlinear form F(u; q) = 0 (corresponding to Equation 26). Thus, the Jacobian and the trial function will be automatically computed and solved. By solving the diffusion equation (Equation 28 and Equation 29), the concentration and the chemical potential can be obtained. Then the radial displacement can be obtained by solving the weak form of the mechanical equilibrium equation (Equation 27). The unknown field variables (c*, μ, u) can be solved by using this staggered scheme at every time step.

### 2.1.3. Numerical results

Some typical phase-transformation cathode materials are LiFePO$_{4}$, LiCoO$_{2}$ and LiMn$_{2}$O$_{4}$. Although all these materials undergo phase transformations, it should be noted that LiFePO$_{4}$ has a clear concentration jump during lithiation between the different phases and parameter χ in Equation 8 must be greater
than 2 to ensure phase separation. However, LiCoO₂ and LiMn₂O₄ undergo a structural transformation and there is no clear concentration jump. In such cases, χ must be less than 2 and in particular, it can be set as 0.

<table><caption>Table 2: Material properties of LiFePO₄</caption><thead><tr><th>Name</th><th>Symbol and unit</th><th>Value</th></tr></thead><tbody><tr><td>Diffusion coefficient</td><td>D[m²/s]</td><td>10<sup>-14</sup></td></tr><tr><td>Young's modulus</td><td>E[GPa]</td><td>124.5</td></tr><tr><td>Poisson's ratio</td><td>ν</td><td>0.25</td></tr><tr><td>Partial molar volume</td><td>Ξ [m³/mol]</td><td>2.9 × 10<sup>-6</sup></td></tr><tr><td>Gradient coefficient</td><td>κ̄ [J/m]</td><td>5.02 × 10<sup>-10</sup></td></tr><tr><td>Expansion coefficient</td><td>α</td><td>2.21 × 10<sup>-2</sup></td></tr><tr><td>Maximum concentration</td><td>c<sub>max</sub>[mol/m³]</td><td>2.29 × 10<sup>4</sup></td></tr></tbody></table>

In the present simulations, the following dimensionless/normalized variables are used as shown in Table 3, where $t_0 = \frac{\rho^2}{D}$ is the characteristic time, F is Faraday's constant, and z is the valence number of a Li-ion ($z=1$). As an example, the stress evolution of a spherical LiFePO₄ electrode is studied. The material properties of LiFePO₄ electrode are given in Table 2 Zhang & Kamlah (2018). The initial radius of the electrode is $\rho = 100$ nm. Normalized radius is 1 and total number of the mesh elements of the unit interval is 100. (It follows that the model results can be applied to any electrode material since the variables are normalized. The general trend is of interest here and not the specific values.) The parameters are as follows: $E_0=124.5$ GPa, $\nu=0.24$, $M^* = 1$, $i^* = 1$, $\kappa^* = 0.0004$, $k=0$, $\Xi=3.49 \times 10^{-6}$ m³/mol, $\chi = 2.6$, $c_{max}=2.29 \times 10^4$ mol/m³, $T=300$ K, $\Delta t=0.0005$, $\alpha=0.0221$. These parameters are used to study the evolution of the concentration, radial stress, hoop stress, hydrostatic stress, von Mises stress, as well as the evolution of the hoop stress on the outer surface, as shown in Figure 1. Only the DISs during the lithium insertion process are considered and a constant Young's modulus is assumed for the two phases during lithiation.
<table><caption>Table 3: Normalized/dimensionless variables</caption><thead><tr><th>Name</th><th>Symbol</th><th>Normalized variable</th></tr></thead><tbody><tr><td>Normalized gradient co-efficient</td><td>κ*</td><td>κc<sub>max</sub> / (RTρ<sup>2</sup>)</td></tr><tr><td>Normalized mobility</td><td>M*</td><td>MRTt<sub>0</sub>/ρ<sup>2</sup></td></tr><tr><td>Normalized time</td><td>t*</td><td>t/t<sub>0</sub></td></tr><tr><td>Normalized radial distance</td><td>r*</td><td>r/ρ</td></tr><tr><td>Normalized stress components</td><td>σ<sub>ij</sub>*</td><td>σ<sub>ij</sub>/E<sub>0</sub></td></tr><tr><td>Normalized chemical potential</td><td>μ*</td><td>μ/(RT)</td></tr><tr><td>Normalized ion flux</td><td>i*</td><td>iρ/(zFDc<sub>max</sub>)</td></tr><tr><td>Normalized partial molecular volume</td><td>Ξ*</td><td>ΞE/(RT)</td></tr></tbody></table>

It takes approximately 6000 $\Delta t$ to complete the lithium insertion process, i.e. reach 100% lithiation of the active particle. Figure 1(a) shows the radial distributions of the normalized Li concentration at lithiation times 500$\Delta t$, 2000$\Delta t$ and 4000$\Delta t$, respectively. A sharp interface is present between the Li-rich and Li-poor phases. At the initial stage (500$\Delta t$), the differences in the concentration profiles are not distinct, and the stresses and strains (Figure 1(b)-(h)) have similar radial distributions. As lithiation continues 2000$\Delta t$ and 4000$\Delta t$, the Li-poor phase tends to have a higher concentration when compared with the case without CM coupling, and the phase interface tends to move faster in the coupling case.

The radial stress (Figure 1(b)), is almost constant in the core (Li-poor) region and gradually decreases to zero on the outer surface. The hoop stress (Figure 1(c)) is also constant in the core-region and transits from tensile in the core to compressive in the shell. The hydrostatic stress (Figure 1(d)) shows a similar trend. The Li-poor region is under constant tension while the Li-rich shell is in compression. From Equation (21), it can be seen that the compressive
<figure><img src="image_2.png" /><figcaption>Figure 1: Radial distribution of (a) Concentration profiles, (b) Normalized radial stress, (c) Normalized hoop stress, (d) Normalized hydrostatic stress, (e) Normalized von Mises stress at time steps 500Δt, 2000Δt, and 4000Δt with/without chemo-mechanical coupling effects. The evolution of (f) Normalized hoop stress on the outer surface. All the stress components are normalized by the Young's modulus of the delithiated electrode material. The unit for the time step is Δt. The solid and dashed lines are for with and without chemo-mechanical coupling effect respectively.</figcaption></figure>

stress leads to a higher chemical potential and the tensile stress leads to a lower chemical potential. Since the chemical potential is the driving force for diffusion, this hydrostatic stress state facilitates the diffusion process. It's reflected on the concentration profiles, showing that the phase-boundary interface moves faster in the CM coupling case. For the von Mises stress (Figure 1(e)), it's
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

Figure 2: Radial distribution of (a) Concentration profiles, (b) Normalized radial stress, (c) Normalized hoop stress, (d) Normalized hydrostatic stress, (e) Normalized von Mises stress, and (f) Radial strain, (g) Hoop strain and (h) Equivalent strain at time step 2000 $\Delta t$ with chemo-mechanical coupling effects under elastic deformation. Constant modulus, reduced modulus and increased modulus during ion insertion are considered. All the stress components are normalized by Young's modulus of the delithiated electrode material. The dotted and dashed lines are for a 10% reduction or increase in the elastic modulus during ion-insertion respectively. As for the solid lines, a constant elastic modulus is assumed.
almost zero in the core region, and there's an abrupt change across the interface. The maximum von Mises stress occurs close to the interface due to the volume mismatch between the two phases. Away from the interface towards the outer surface, the von Mises stress gradually decreases. It should be noted that higher radial, hoop, hydrostatic, and von Mises stresses are predicted for the coupling cases when compared with the ones without CM coupling at the same time step. This is due to the fact that equilibrium concentrations are almost the same for the Li-rich phase while a higher equilibrium concentration is predicted for the coupling case. The coupling model predicts a smaller miscibility gap (concentration differences between the Li-rich phase and Li-poor phase), thus, less volume change and DISs are induced.

In order to study the influence of material property changes due to lithiation, the parameter *k* is set to be -0.1*E*₀ and +0.1*E*₀ (corresponds to a 10% reduction or increase in modulus with increasing ion concentration respectively). All the other parameters are the same and CM coupling is considered. Since the stress and strain profiles have the similar trends, only the profiles at at time step 2000 Δ*t* are plotted and analyzed as shown in Figure 2. It should be noted that when a concentration dependent modulus is considered, the concentration profile (Figure 2(a)) is almost the same. Similarly, the distributions of the stresses (Figure 2(b)-(e)) are the same but the magnitudes are different. This is due to the fact that a higher stress is induced for materials with an increasing modulus as a function of Li-concentration when compared to the cases where the modulus is either not dependent on concentration or it decreases as a function of it. Also, the stresses predicted from the coupling case are smaller than those from the non-coupling case when the same material parameters are used. This is because the coupling model predicts a smaller miscibility gap (concentration differences between the Li-rich phase and Li-poor phase).

To gain a better understanding of the physics behind the stress evolution, let's consider the lithiation process. During Li- insertion, the outer shell contains more Li-ions than the inner core and, as a result, the core expands to a lesser degree. The expansion is constrained by the exterior region, leading to
<figure><img src="image_11.png" /><figcaption>Figure 3: Radial distribution of (a) Normalized hoop stress in the particle center, and (c) Normalized hoop stress at the outer surface at time step 2000 $\Delta t$ with chemo-mechanical coupling effects under elastic deformation. Both of the modulus increasing or decreasing during ion insertion are considered. All the stress components are normalized by the Young's modulus of the delithiated electrode material. The dotted and dashed lines are for a 10% reduction or increase in the elastic modulus during ion-insertion respectively. As for the solid lines, a constant elastic modulus is assumed.</figcaption></figure>

a compressive hoop stress in this region. The traction-free boundary condition requires that the radial stress vanishes at the outer surface at all times. The volume mismatch between the core and shell leads to tensile stresses in both the radial and hoop directions. Due to a constant Li concentration and spherical symmetry, the radial stress equals the hoop stress in the core regions. As lithiation proceeds, the radial and hoop stresses in the core continue to grow in the elastic case due to the growing Li concentration in Li-poor regions. However, the hoop stress near the outer shell decreases as insertion progresses, due to the movement of the phase interface towards the center. Consequently, the Li-poor phase restricts the swelling of the Li-rich phase near the exterior to a lesser degree.

The strain distributions are also investigated. The radial strain (Figure 2(f)) is almost constant in the core region and gradually increases towards the outer surface, while the hoop strain (Figure 2(g)) and equivalent strain (Figure 2(h)) show similar trends as those of the von Mises stress distribution. An interesting finding is that the magnitudes of the strains are higher when the modulus de-
creases during ion insertion, as opposed to when it increases with concentration.

Temporal distributions of the stresses in the center (Figure 3(a)) and on the outer surface (Figure 1(f) and Figure 3(b)) are recorded. Radial stress in the particle center is equal to hoop stress in the particle center due to spherical symmetry. It can be seen that the radial/hoop stress in the center increase faster in the early lithiation stages than in the later stages. While for the hoop stress on the outer surface, the magnitude increases quickly during early lithiation stages and then gradually decreases. The radial and hoop stresses in the center are tensile while the hoop stress on the outer surface is always compressive. Due to the traction-free boundary conditions on the outer surface, the radial stress on the outer surface vanishes (i.e., is always zero). As a result, cracks are not likely to occur and propagate during lithiation.

## 2.2. DISs for elasto-plastic materials

### 2.2.1. Phase field model for elasto-plastic materials

In alloy-based electrode systems, such as Li/Si, a large volumetric change may be induced (over 400% for Li/Si), which can give rise to plastic deformation. If infinitesimal changes occur in the stress and chemical composition, a unique infinitesimal change of strain is produced, and the change is independent of the speed applied. In this case, we can say that the chemo-plasticity model is rate-independent. The model can be formulated from an energy perspective. Generally speaking, the formulation can be either in terms of stress space or strain space. The constitutive relation (flow rule) for the plastic strain rate is usually given in terms of the stress space. As a result, the stress-based formulation is employed here.

The Gibbs free energy $\psi = \psi(\sigma, c^*, \xi)$ is a function of the stress tensor ($\sigma$), chemical composition ($c^*$) and internal (intrinsic) variable ($\xi$). In this notation, the chemoelastic part of the strain tensor is given by:

$$ \epsilon^{ec} = \frac{\partial \psi}{\partial \sigma} \qquad (32) $$
The yield condition in stress space can be expressed as:

$$f(\sigma, c^*, \xi) = 0 \qquad (33)$$

And the chemoplastic strain rate is given by:

$$\dot{\epsilon}^p = \dot{\gamma} \frac{\partial g}{\partial \sigma} \qquad (34)$$

where $g = g(\sigma, c^*, \xi)$ is the plastic flow potential and $\dot{\gamma}$ is the rate of the plastic multiplier. Since plastic deformation is an irreversible process, the multiplier must be non-negative. For the associative flow rule, $g = f$.

Kuhn-Tucker conditions must be satisfied during the deformation process, which are stated as:

$$f \le 0, \dot{\gamma} \ge 0, \dot{\gamma}f = 0 \qquad (35)$$

This consistency condition is also equivalent to:

$$\dot{f} \le 0, \dot{\gamma} \ge 0, \dot{\gamma} \dot{f} = 0 \qquad (36)$$

The hardening behavior of the material is determined by $h$ and it's a function of $c^*$, and the history variable $(\vartheta)$. The hardening rule can be expressed as:

$$\dot{\xi} = \dot{\gamma} h(c^*, \vartheta) \qquad (37)$$

and the hardening parameter $H_o$ is defined as:

$$H_o = -\frac{\partial f}{\partial \xi} h(c^*, \vartheta) \qquad (38)$$

The chemoelastic strain rate is given by:

$$\dot{\epsilon}^{ec} = \frac{\partial^2 \psi}{\partial \sigma \otimes \partial \sigma} : \dot{\sigma} + \frac{\partial^2 \psi}{\partial \sigma \partial c^*} \dot{c}^* \qquad (39)$$

The general form of $\psi$ is as follows:

$$\psi = \psi^e(\sigma) + \psi^{ch}(\sigma, c^*) + \psi^p(c^*, \vartheta) \qquad (40)$$

The three terms represent the energy contribution from mechanical elastic part, chemical elastic part, and plastic part respectively. The mechanical elastic part and plastic part are given as:

$$\psi^e(\sigma) = \frac{1}{4\Lambda} \left( \text{tr}(\sigma^2) - \frac{\lambda}{3\Lambda + 2\Lambda} \text{tr}^2(\sigma) \right) \qquad (41)$$
$$ \psi^p(c^*, \vartheta) = \frac{1}{2} H \zeta^2 + \sigma_y^0 \zeta \quad (42) $$

Linear isotropic hardening is assumed with $H$ being the hardening modulus, and $\zeta$ is a strain-like internal variable (normally it's equivalent to plastic strain). As a result, the chemoelastic part of the strain rate could be given as:

$$ \dot{\boldsymbol{\epsilon}}^{ec} = \frac{1}{2\Lambda} \left( \mathbb{I} - \frac{\lambda}{3\Lambda + 2\Lambda} \boldsymbol{I} \otimes \boldsymbol{I} \right) : \dot{\boldsymbol{\sigma}} + \frac{\partial \psi^{ch}}{\partial c^*} \dot{c}^* \quad (43) $$

For a homogenous and isotropic material, a linear expansion is assumed,

$$ \psi^{ch}(\boldsymbol{\sigma}, c^*) = \eta(c^*) \operatorname{tr}(\boldsymbol{\sigma}) \quad (44) $$

$\eta$ is given by:

$$ \eta = \alpha \Delta c^* = \alpha (c^* - c_0^*) \quad (45) $$

where $c_0^*$ is a reference concentration, usually set as 0. The consistency condition is $\dot{f} = 0$, i.e.,

$$ \dot{f} = \frac{\partial f}{\partial \boldsymbol{\sigma}} : \dot{\boldsymbol{\sigma}} + \frac{\partial f}{\partial \xi} : \dot{\xi} + \frac{\partial f}{\partial c^*} \dot{c}^* = 0 \quad (46) $$

The consistency parameter is given by:

$$ \dot{\gamma} = \frac{\left\langle \frac{\partial f}{\partial \boldsymbol{\sigma}} : \mathbb{C} : \dot{\boldsymbol{\epsilon}} - \frac{\partial f}{\partial \boldsymbol{\sigma}} : (3\lambda + 2\Lambda)\alpha \dot{c}^* \boldsymbol{I} + \frac{\partial f}{\partial c^*} \dot{c}^* \right\rangle}{H_0 + \frac{\partial f}{\partial \boldsymbol{\sigma}} : \mathbb{C} : \frac{\partial g}{\partial \boldsymbol{\sigma}}} \quad (47) $$

With $\langle x \rangle = \frac{1}{2}(x + |x|)$ denoting the positive part of $x$. For $J_2$ plasticity and associative flow rule, $\frac{\partial f}{\partial \boldsymbol{\sigma}} : \boldsymbol{I} = \frac{s}{\|s\|} : \boldsymbol{I} = \frac{\operatorname{tr}(s)}{\|s\|} = 0$, therefore,

$$ \dot{\gamma} = \frac{\left\langle \frac{\partial f}{\partial \boldsymbol{\sigma}} : \mathbb{C} : \dot{\boldsymbol{\epsilon}} + \frac{\partial f}{\partial c^*} \dot{c}^* \right\rangle}{H_0 + 2\Lambda} = \frac{\left\langle \frac{\partial f}{\partial \boldsymbol{\sigma}} : \mathbb{C} : \dot{\boldsymbol{\epsilon}} + \frac{\partial f}{\partial c^*} \dot{c}^* \right\rangle}{\frac{2}{3}H + 2\Lambda} \quad (48) $$

When plasticity is taken into consideration, the constitutive relations for a spherical particle can be given by:

$$ \left. \begin{aligned} \sigma_r &= \frac{E}{(1+\nu)(1-2\nu)} \left[ (1-\nu) \left( \frac{du}{dr} - \varepsilon_r^p - \alpha c^* \right) + 2\nu \left( \frac{u}{r} - \varepsilon_\theta^p - \alpha c^* \right) \right] \\ \sigma_\theta &= \frac{E}{(1+\nu)(1-2\nu)} \left[ \left( \frac{u}{r} - \varepsilon_\theta^p - \alpha c^* \right) + \nu \left( \frac{du}{dr} - \varepsilon_r^p - \alpha c^* \right) \right] \end{aligned} \right\} \quad (49) $$
And the deviatoric parts are:

$$ \sigma_r' = \frac{2}{3}(\sigma_r - \sigma_\theta), \sigma_\theta' = \sigma_\phi' = -\frac{1}{3}(\sigma_r - \sigma_\theta) \quad (50) $$

As a result, the effective stress/von Mises stress is:

$$ \sigma_e = \sqrt{\frac{3}{2} \mathbf{s} : \mathbf{s}} = |\sigma_r - \sigma_\theta| \quad (51) $$

When linear hardening is assumed, the yield function f and yield stress $\bar{\sigma}$ are given by:

$$ f(\sigma, \xi) = \sqrt{\mathbf{s} : \mathbf{s}} - \sqrt{\frac{2}{3}\bar{\sigma}(\xi)} \quad (52) $$

$$ \bar{\sigma}(\xi) = \bar{\sigma}(\bar{\varepsilon}_p) = \sigma_y^0 + H_0 \bar{\varepsilon}_p \quad (53) $$

where $\bar{\varepsilon}_p$ is the equivalent plastic strain, defined as below:

$$ \bar{\varepsilon}_p = \sqrt{\frac{2}{3} \varepsilon^p : \varepsilon^p} \quad (54) $$

The plastic strain components are given by:

$$ \dot{\varepsilon}_r^p = \dot{\gamma} \frac{\sigma_r'}{\|\mathbf{s}\|}, \quad \dot{\varepsilon}_\theta^p = \dot{\gamma} \frac{\sigma_\theta'}{\|\mathbf{s}\|} \quad (55) $$

For the associative flow rule,

$$ \frac{\partial g}{\partial \boldsymbol{\sigma}} = \frac{\partial f}{\partial \boldsymbol{\sigma}} = \operatorname{sign}(\sigma_r') \frac{1}{\sqrt{6}} \operatorname{diag}(2, -1, -1) \quad (56) $$

$$ \begin{aligned} \frac{\partial f}{\partial c^*} &= \frac{\partial f}{\partial \boldsymbol{\sigma}} : \frac{\partial \boldsymbol{\sigma}}{\partial c^*} + \frac{\partial f}{\partial \sigma_y^0} \frac{\partial \sigma_y^0}{\partial c^*} + \frac{\partial f}{\partial H_0} \frac{\partial H_0}{\partial c^*} \\ &= 2 \frac{\partial \Lambda}{\partial c^*} \frac{\partial f}{\partial \boldsymbol{\sigma}} : \boldsymbol{\varepsilon} - 2 \frac{\partial \Lambda}{\partial c^*} \frac{\partial f}{\partial \boldsymbol{\sigma}} : \boldsymbol{\varepsilon}^p - \sqrt{\frac{2}{3}} \left( \frac{\partial \sigma_y^0}{\partial c^*} + \frac{\partial H_0}{\partial c^*} \bar{\varepsilon}_p \right) \end{aligned} \quad (57) $$

When the yield function does not depend on the concentration field the consistency parameter can be determined by:

$$ \dot{\gamma} = \frac{\sqrt{6\Lambda}}{H + 3\Lambda} \operatorname{sign}(\sigma_r') \left( \frac{\partial \dot{u}}{\partial r} - \frac{\dot{u}}{r} \right) \Theta\left(\frac{\sigma_e}{\bar{\sigma}} - 1\right) \quad (58) $$

With $\Theta$ being the Heaviside unit step function, whose value is 1 for any positive variable and 0 otherwise. And $\sigma_e = |\sigma_r - \sigma_\theta|$ is the effective/ von Mises stress.
### 2.2.2. Numerical implementation of the phase field model

The problem becomes to seek solutions of $c_{n+1}^*$, $\mu_{n+1}$, $u_{n+1}$, and $\gamma_{n+1}$ at the time step $t_{n+1}$ given that $c_n^*$, $\mu_n$, $u_n$, and $\gamma_n$ at the time step $t_n$ are known. Since the weak form of the diffusion equation and the mechanical equilibrium equation are the same as in the elastic case, the same approach can be applied to the elasto-plastic case. The only difference comes from the plastic deformation. The current model does not take kinematic hardening and viscous effects into account, i.e. only the quasi-static problem is solved. When the time step is small enough, the increment of the consistency parameter can be given by:

$$ \Delta\gamma_n = \frac{\sqrt{6\Lambda}}{H + 3\Lambda}\mathrm{sign}((\sigma_r)_n)' \left(\frac{\partial(\Delta u_n)}{\partial r} - \frac{\Delta u_n}{r}\right) \Theta\left(\frac{(\sigma_e)_n}{\bar{\sigma}_n} - 1\right) \quad (59) $$

where

$$ \Delta u_n = u_{n+1} - u_n, $$

$$ (\sigma_r)'_n = \frac{2}{3}[(\sigma_r)_n - (\sigma_\theta)_n], \quad (\sigma_\theta)'_n = -\frac{1}{3}[(\sigma_r)_n - (\sigma_\theta)_n], $$

$$ (\sigma_e)_n = |(\sigma_r)_n - (\sigma_\theta)_n|, $$

$$ \begin{aligned} \bar{\sigma}_n &= \sigma_y^0 + H (\bar{\varepsilon}_p)_n, \\ (\bar{\varepsilon}_p)_n &= \sqrt{\frac{2}{3} (\varepsilon^p)_n : (\varepsilon^p)_n}, \\ \|\mathbf{s}\|_n &= \sqrt{\frac{2}{3} (\sigma_e)_n}, \\ (\Delta \varepsilon_r^p)_n &= \Delta \gamma_n \frac{(\sigma_r')'_n}{\|\mathbf{s}\|_n}, \quad (\Delta \varepsilon_\theta^p)_n = \Delta \gamma_n \frac{(\sigma_\theta')'_n}{\|\mathbf{s}\|_n}. \end{aligned} $$

As a result,

$$ \begin{aligned} \gamma_{n+1} &= \Delta\gamma_n + \gamma_n, \\ (\varepsilon_r^p)_{n+1} &= (\Delta\varepsilon_r^p)_n + (\varepsilon_r^p)_n, \\ (\varepsilon_\theta^p)_{n+1} &= (\Delta\varepsilon_\theta^p)_n + (\varepsilon_\theta^p)_n, \end{aligned} $$

and the other quantities at time step $t_{n+1}$ can also be computed. The mechanical equilibrium equation is essentially the same as for the elastic case, and the same solver is used. Once the internal iterations for solving the weak form of the mechanical equilibrium equation (Equation 26) converge, one increment of the mechanical field can be obtained by forward Euler integration method. Thus, the mechanical equilibrium equation can be solved, and the stress field can be obtained incrementally.
### 2.2.3. Numerical results

Here the stress evolution in a spherical Si electrode with radius 100 nm is studied. It has been shown that both the elastic modulus and the yield stress of Si are reduced with continuous Li-insertion Sethuraman et al. (2010), however, this concentration-dependence has not been accounted for theoretically. Furthermore, previous studies used either linear hardening or viscoplasticity to study the stress evolution. However, few looked at strain softening effects, which most likely take place based on the experimental observations Sethuraman et al. (2010); Soni et al. (2012). In the present study, both strain hardening and strain softening are considered. For the sake of simplicity, concentration-dependence of the modulus is not considered in this section (but in the next one), and therefore here the elastic modulus is kept the same for all cases. Normalized variables as in the elastic case are employed. The normalized radius is 1 and total number of the element mesh of the unit interval is 100. The other parameters are given by: $E_0=90$ GPa, $v=0.24$, $M^* = 1$, $i^* = 1$, $\kappa^* = 0.0004$, $\chi = 2.6$, $k = 0$ (same Young's modulus for both the core and the shell), $\sigma_y^0 = 0.05E_0$ Chen et al. (2014), $H = \pm0.01E_0$ (corresponds to strain hardening and softening respectively) Chen et al. (2014), $\Xi=3.1 \times 10^{-5}$ m$^3$/mol, $c_{max}=3.67 \times 10^5$mol/m$^3$, $T=300$ K, $\Delta t=0.0005$, $\alpha=\sqrt[3]{4.2}-1=0.613$ (corresponding to 420% volume expansion after full insertion).

At the early stage of lithiation (Figure 4), the differences in concentration, stresses and strain distributions are small. The concentration profiles (Figure 4(a)) almost coincide with each other. The radial stresses (Figure 4(b)) are almost constant in the core region and gradually decrease to zero as they approach the shell; the hoop stresses (Figure 4(c)) transit from tension in the core to compression in the shell; and the hydrostatic stresses (Figure 4(d)) show similar trends as the hoop stresses. The resulting von Mises stresses (Figure 4(e)) gradually increase from the core towards the shell, and as soon as they exceed the yield stress, plastic deformation begins to take place. If there's no plastic deformation, the stress evolution will be the same as in the elastic deformation.
<figure><img src="image_12.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(d)</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(e)</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(f)</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(g)</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(h)</figcaption></figure>

<figure><img src="image_19.png" /></figure>

Figure 4: Radial distribution of (a) Concentration profiles, (b) Normalized radial stress, (c) Normalized hoop stress, (d) Normalized hydrostatic stress, (e) Normalized von Mises stress, (f) Radial plastic strain, (g) Hoop plastic strain and (h) Equivalent plastic strain at time step 100 $\Delta t$ with/without chemo-mechanical coupling ('w cp' / 'w/o cp') effects under elasto-plastic deformation. Both strain hardening and strain softening cases are considered. All the stress components are normalized by Young's modulus of the delithiated electrode material.

case, i.e., the maximum von Mises stress will occur near the core-shell interface. However, once it exceeds the yield stress, there will be plastic deformation. As
a result, the von Mises stress (Figure 4(e)) is very small in the core and reaches a local maximum near the core-shell interface, and then there's a stress drop to a local minimum near the local maximum position but the von Mises stress increases rapidly towards the outer surface afterwards. The local stress changes are also reflected in the hoop stresses (Figure 4(c)) and the hydrostatic stresses (Figure 4(d)). As lithiation progresses, the phase interface (Figure 4(a), Figure 5(a), and Figure 6(a)) moves from the shell exterior towards the center. The coupling also predicts a smaller miscibility gap. But contrary to what happens in elastic electrodes, the phase boundary moves slower than the non-coupling case. This is in agreement with the in situ experiments for crystalline Si nanoparticles, where it is shown that a slower migration velocity occurs as lithiation progresses McDowell et al. (2012). In later lithiation stages, the radial stress (Figure 4(b), Figure 5(b), and Figure 6(b)), hoop stress (Figure 4(c), Figure 5(c), and Figure 6(c)) and hydrostatic stress (Figure 4(d), Figure 5(d), and Figure 6(d)) show similar trends for the three cases considered: strain hardening without CM coupling ('w/o cp'), strain hardening with CM coupling ('w cp') and strain softening without CM coupling. They all (stress/strain profiles) remain constant in the core regions with strain hardening without CM coupling being the largest and strain softening without coupling being the lowest in magnitude. The radial stress is compressive in both the core and shell, while the hoop and hydrostatic stresses are compressive in the core and transit from compression near the interface to tension in the shell exterior. This hydrostatic stress state hinders the diffusion of Li-ions, which is the reason for the slower migration of the reaction front. There are stress drops in the radial stress, hoop stress and hydrostatic stress profiles near the core-shell interface due to the occurrence of plastic deformation. The von Mises stress profiles (Figure 5(e), and Figure 6(e)) also show stress drops but the situation is a little bit different. The von Mises stresses (Figure 4(e), Figure 5(e), and Figure 6(e)) are almost zero in the core regions. As a result, the core remains elastic and no plastic deformation (Figure 4(f)-(h), Figure 5(f)-(h), and Figure 6(f)-(h)) occurs there. The von Mises stresses show abrupt changes near core/shell interfaces as they
increase first, then decrease and increase again quickly. The von Mises stress
towards the shell continues to increase for the strain hardening case while it
decreases for the strain softening case.

As for the plastic deformation, the plastic strain components are zero in
the core regions. The radial component of the plastic strain for the shell is
always positive while the hoop component is always negative. The equivalent
plastic strain continues to increase as lithiation continues, and is higher in the
strain softening case. It's interesting to notice that there's a relative plateau
in the equivalent plastic strain profiles (Figure 5(h), and Figure 6(h)). This
is due to the stress drops near the core-shell interface. The von Mises stress
near the phase boundary decreases quickly and then increases rapidly again.
However, according to the von Mises yield criteria, plastic deformation can
continue to occur only if the von Mises stress exceeds the effective yield stress.
Thus, there won't be plastic strain increment over the stress drop regions. The
radial and hoop stresses in the center (Figure 7(a)) are equal due to spherical
symmetry. For the initial stages (up to 1000 $\Delta t$), the stress profiles (Figure 7(a))
are almost the same and all in tension. However, they all undergo a transition
from tension to compression and the magnitude is the highest when it's strain
hardening without CM coupling and the magnitude is the smallest when it's
strain softening with CM coupling. The hoop stress is the driving force for crack
propagation and on the outer surface (Figure 7(b)) it is compressive during
the early lithiation stages and transits from compression to tension at later
stages. The hoop stresses on the outer surface are almost the same for the initial
stages when strain hardening/softening is considered either with or without CM
coupling. However, as soon as the electrode begins to yield, the hoop stress on
the outer surface keeps increasing for the strain hardening case while it decreases
for the strain softening case.

To have a better understanding of the stress evolution, the constitutive re-
lation (Equation 49) may be revisited. The stress state depends on both
lithiation-induced strain components and plastic strain components. In the
elastic case, the hoop stress in the shell is always compressive. When plastic
<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /><figcaption>Figure 5: Radial distribution of (a) Concentration profiles, (b) Normalized radial stress, (c) Normalized hoop stress, (d) Normalized hydrostatic stress, (e) Normalized von Mises stress, (f) Radial plastic strain, (g) Hoop plastic strain and (h) Equivalent plastic strain at time step 1500 Δt with/without chemo-mechanical coupling ('w cp' / 'w/o cp') effects under elasto-plastic deformation. Both strain hardening and strain softening cases are considered. All the stress components are normalized by Young's modulus of the delithiated electrode material.</figcaption></figure>

deformation is considered, the stress state depends on the competition between
the strain components. The lithiation-induced strain is isotropic and positive
during insertion for the elastic case. While for plasticity, the core remains elastic and yielding occurs in the shell region, with the radial plastic strain always being positive, (as can be seen from Figure 6(f), Figure 4(f), Figure 5(f)) the hoop plastic strain always being negative (as can be seen from Figure 6(g), Figure 4(g), Figure 5(g)), and $\epsilon_{\theta}^{p} = -\frac{1}{2}\epsilon_{r}^{p}$. This means that the plastic deformation leads to radial expansion and tangential shrinkage. As a result, the hoop stress can be given by:

$$ \sigma_{\theta} = \frac{E}{(1 + \nu)(1 - 2\nu)} \left[ \left( \frac{u}{r} - \alpha c^{*} \right) + \nu \left( \frac{du}{dr} - \alpha c^{*} \right) \right] + \frac{E}{2(1 + \nu)} \epsilon_{r}^{p} \quad (60) $$

The first part on the right-hand side of Equation 60 is the solution for the hoop stress under elastic deformation only, arising from mismatch strains due to the concentration gradient. The region close to the outer surface (Li-rich region/shell) contains significantly more Li-ions than the inner region (Li-poor region/core). As a result, the Li-rich regions expand to a larger degree. However, this expansion is constrained by the surrounding materials and the hoop stress in the shell is compressive initially. Due to the phase transformation, large stresses develop near the phase interface. When the effective stress (von Mises stress) is large enough, plastic deformation emerges near the phase interface. The second part of Equation 60 represents the contribution due to plastic deformation and it's always positive in the shell. As lithiation proceeds, the plastic region propagates and plastic deformation accumulates. When the stress component due to plastic deformation surpasses the component due to the concentration gradient, the hoop stress in the shell undergoes a transition from compression to tension. This is believed to be the cause of experimentally observed cracks during the lithiation process of spherical Si/Sn electrode particles. For example in paper Aifantis et al. (2012) it was shown that Sn and Sn nanoparticles severely fractured upon the first lithiation. Significant fracture was also observed in Si nanoparticles Liu et al. (2012); McDowell et al. (2012) during the first cycle upon Li-insertion. Particularly, in Liu et al. (2012); McDowell et al. (2012) the cracks were observed to propagate inward along the radial direction during the lithiation process. The transition of tangential com-
<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /></figure>

<figure><img src="image_31.png" /></figure>

<figure><img src="image_32.png" /></figure>

<figure><img src="image_33.png" /></figure>

<figure><img src="image_34.png" /></figure>

<figure><img src="image_35.png" /><figcaption>Figure 6: Radial distribution of (a) Concentration profiles, (b) Normalized radial stress, (c) Normalized hoop stress, (d) Normalized hydrostatic stress, (e) Normalized von Mises stress, (f) Radial plastic strain, (g) Hoop plastic strain and (h) Equivalent plastic strain at time step 3000 Δt with/without chemo-mechanical coupling ('w cp' / 'w/o cp') effects under elasto-plastic deformation. Both strain hardening and strain softening cases are considered. All the stress components are normalized by Young's modulus of the delithiated electrode material.</figcaption></figure>

pression to tangential tension explains the surface cracking of electrodes that
occurs experimentally upon maximum lithiation.
It should be pointed out that this model is based on the small deformation framework, so the numerical results cannot be compared quantitatively with experimental measurements. For quantitative comparison purposes, models with finite strain deformation are more appropriate. Our conclusions are in fact consistent with those obtained using finite strain Chen et al. (2014). Even though the predictions of our model are not quantitative, the main effects are captured and the model results can provide a complete overview of the different scenarios that can occur depending on the model parameters, and provide a comparison between the different cases, something which is more difficult to do in finite strain and has not been done before.

<figure><img src="image_36.png" /><figcaption>Figure 7: The evolution of (a) Normalized hoop stress in the particle center and (b) Normalized
hoop stress at the outer surface under elasto-plastic deformation. All the stress components
are normalized by Young's modulus of the delithiated electrode material. The unit for time
step is $\Delta t$.</figcaption></figure>

## 2.3. Comparison between elastic and plastic deformation

To consider the effects of plasticity on DISs, the stress evolution under elastic deformation and elasto-plastic deformation in a spherical Si electrode during lithiation is calculated. To make the comparison less complicated, CM coupling is not considered here. As a result, the same concentration distribution (Figure 8(a)) is obtained when the same set of parameters are employed. At early stages (100 $\Delta t$), phase equilibrium has not been reached. As a result, the Li concentration is higher near the outer surface and there are almost no
Li-ions in the interior region. As lithiation continuous (900 Δt, and 2000 Δt), a core-shell-like structure forms with the Li concentration being close to 1 and 0, respectively. The core and shell are separated by a sharp interface. As lithiation progresses, the interface also moves towards the center of the particle. And the Li concentration in the shell increases slightly. A reduction in the elastic modulus during Li-insertion is also taken into account (based on Equation 57). The parameters are the same as in the previous section except that $k = -\frac{5}{9}E_0$ and $H = 0.01E_0$. For plastic deformation, only strain hardening is considered. When the material properties of the core and shell are taken as being the same, the stress evolutions are as shown in Figure 8(b),Figure 8(d),and Figure 8(f). At the initial stage, the deformation of the electrode remains elastic. The radial and hoop stresses are the same in the core region. The radial stress gradually decreases to zero towards the outer surface, and there is an abrupt transition in the hoop stress from tension in the core to compression in the shell. The maximum von Mises stress occurs near the interface.

Under elasto-plastic deformation, yielding begins to occur near the outer surface during the initial lithiation stages, and the stress evolutions are almost the same as in the previous section. The hoop and radial stresses in the core transits from tension to compression at later stages, while the opposite trend is seen for the shell. The von Mises stress is almost zero in the core and increases to a local maximum first, then decreases to a local minimum and finally keeps increasing towards the shell.

When a reduction in the elastic modulus with increasing concentration is considered, similar distributions (Figure 8(c), Figure 8(e), and Figure 8(g)) are found for the stresses. Comparison between the DISs with and without concentration-dependent elastic moduli are also made under elastic (Figure 9(a)-(c)) and plastic (Figure 9(d)-(f)) deformation, respectively. It is found that when a reduction in the modulus with increasing concentration is considered, much smaller stresses are predicted in the elastic case, while for the plastic case, higher stresses are predicted. The stress evolution over time is also investigated, and it should be noted that there is a visible difference. The radial and hoop stresses
<figure><img src="image_37.png" /></figure>

Figure 8: (a) Concentration profiles at time step 100 $\Delta t$, 900 $\Delta t$, and 2000 $\Delta t$ respectively. Radial distribution of normalized radial stress ($\sigma_r$), normalized hoop stress ($\sigma_\theta$) and normalized von Mises stress ($\sigma_v$), with constant modulus, under elastic and elasto-plastic deformation at time steps of (b) 100 $\Delta t$, (d) 900 $\Delta t$, and (f) 2000 $\Delta t$ respectively. Radial distribution of normalized radial stress, normalized hoop stress and normalized von Mises stress, with reduction in modulus with increasing concentration, under elastic and elasto-plastic deformation at time steps of (c) 100 $\Delta t$, (e) 900 $\Delta t$, and (g) 2000 $\Delta t$ respectively. (h) The evolution of Normalized hoop stress on the outer surface under elastic ('Ela') and elasto-plastic ('Pla') deformation. All the stress components are normalized by the Young's modulus of the delithiated electrode material. The unit for time step is $\Delta t$. For the stress profiles ((b)-(g)), the solid lines are for the elastic deformation case and the solid lines with asterisks are for the elasto-plastic deformation case.
<figure><img src="image_38.png" /><figcaption>Figure 9: The effects of a reduced elastic modulus with lithiation on the radial distribution of (a)Normalized radial stress ($\sigma_r$), (b) Normalized hoop stress ($\sigma_\theta$), and (c) Normalized von Mises stress ($\sigma_v$) under elastic deformation. The effects of a reduced elastic modulus with lithiation on the radial distribution of (d) Normalized radial stress, (e) Normalized hoop stress, and (f) Normalized von Mises stress under elasto-plastic deformation. The evolution of (g) Normalized radial stress and (h) Normalized hoop stress in the center. All the stress components are normalized by the Young's modulus of the delithiated electrode material. The unit for time step is $\Delta t$. For the stress profiles ((a)-(f)), the solid lines are for the constant elastic modulus case and the dashed lines are for the case with reduction in modulus with increasing concentration.</figcaption></figure>
(Figure 9(g)-(h)) in the center are always positive for elastic deformation while they are only positive for the initial lithiation stages and then become negative at later stages under elasto-plastic deformation. The hoop stress on the outer surface is always compressive in the elastic case while for the plastic case it's compressive only for the initial lithiation stages and transits from compression to tension at later stages. The stresses predicted from the elastic model are much higher than those from plastic model, and the elastic model cannot explain the experimentally observed radial cracks during lithiation. This is due to the fact that the hoop stress near the outer surface is the crack driving force and compression cannot cause cracks to grow. As can be seen in Figure 8(h), the elastic model predicts a compressive hoop stress all the time during lithiation while the plastic model predicts compression only at early stages and it predicts tension at late stages, which allows for crack initiation and is consistent with experimental evidence showing crack formation right before Li de-insertion. As a result, the plastic model is more appropriate for high capacity electrodes coupled with large volume changes and phase transformations.

### 3. Conclusions

In this work, phase field models accounting for DISs in spherical active materials that undergo phase-transformations are developed. For electrodes with relatively small volume variations, elastic models can be employed while for electrodes with large volume changes, plastic models are preferred. The models accounts for the effect of phase. The sharp phase boundary is naturally captured by the phase field model, and the concentration field is obtained by a mixed formulation of the fourth-order Cahn-Hilliard equation. DISs are obtained by solving the variational form of the mechanical equilibrium equations. It is found that the DISs arise from the inhomogeneous volume expansions resulting from Li concentration gradients and the hydrostatic stress facilitates the diffusion of Li-ions under elastic deformation while it hinders diffusion in the plastic case. For elastic electrodes, higher radial, hoop, hydrostatic, and von
Mises stresses are predicted for the coupling cases when compared with the ones without CM coupling at the same time step. This is due to the fact that the coupling model predicts a smaller miscibility gap (concentration differences between the Li-rich phase and Li-poor phase), thus, less volume change and DISs are induced. For elasto-plastic solids, it's the opposite. The coupling model predicts a higher miscibility gap and the electrodes experience a slightly higher von Mises stress and equivalent plastic strain when it's compared with the ones without CM coupling at the same time step. When the lithiation results in a lower elastic modulus, as opposed to having the elastic properties independent of concentration, a decrease in DISs but an increase in strain occur under elastic deformation. The opposite occurs for the plastic case. Under elastic deformation, the radial stress is always tensile while the hoop stress is tensile in the core region and compressive in the shell. For the plastic case, the radial stress shows a transition from tension in initial stages to compression at later stages. The hoop stress in the core region also shows a similar trend while the hoop stress in the shell shows a transition from compression to tension, which is believed to be the cause of experimentally observed surface cracking during the Li insertion process. Furthermore, if strain softening due to plastic deformation is assumed, smaller stresses and higher plastic strains are predicted than the strain hardening case. To sum up, concentration-dependent material properties due to Li insertion and hardening behavior of the material due to plastic deformation plays a significant role on DISs in spherical phase transformation electrodes. By taking these factors into consideration, more accurate predictions of the DISs can be achieved, thus providing an improved theoretical basis and insight for designing next-generation mechanically stable phase transforming electrode materials.

## Acknowledgments

The authors are grateful to the National Science Foundation for supporting this work through the CMMI grant (CMMI-1762602).
References

(). https://fenics-dolfin.readthedocs.io/en/2017.2.0/apis/api_nls.html.

Aifantis, K., Dempsey, J., & Hackney, S. (2005). Cracking in si-based anodes for li-ion batteries. *Rev. Adv. Mater. Sci*, *10*, 403–408.

550

Aifantis, K., & Hackney, S. (2003). An ideal elasticity problem for li-batteries. *Journal of the Mechanical Behavior of Materials*, *14*, 413–427.

Aifantis, K. E., Dempsey, J., & Hackney, S. (2006). Fracture of nanostruc-
tured lithium batteries. In *Fracture of Nano and Engineering Materials and Structures* (pp. 63–64). Springer.

555

Aifantis, K. E., Huang, T., Hackney, S. A., Sarakonsri, T., & Yu, A. (2012). Capacity fade in sn-c nanopowder anodes due to fracture. *Journal of power sources*, *197*, 246–252.

Amatucci, G., Tarascon, J., & Klein, L. (1996). Cobalt dissolution in licoo2-
based non-aqueous rechargeable batteries. *Solid State Ionics*, *83*, 167–173.

560

Chang, L., Lu, Y., He, L., & Ni, Y. (2018). Phase field model for two-phase lithiation in an arbitrarily shaped elastoplastic electrode particle under galvanostatic and potentiostatic operations. *International Journal of Solids and Structures*, *143*, 73–83.

Chen, C.-H., Chason, E., & Guduru, P. R. (2017). Numerical solution of moving phase boundary and diffusion-induced stress of sn anode in the lithium-ion battery. *Journal of The Electrochemical Society*, 164, E3661–E3670.

565

Chen, L., Fan, F., Hong, L., Chen, J., Ji, Y., Zhang, S., Zhu, T., & Chen, L. (2014). A phase-field model coupled with large elasto-plastic deformation: application to lithiated silicon electrodes. *Journal of The Electrochemical Society*, 161, F3164–F3172.

570
Cheng, Y.-T., & Verbrugge, M. W. (2008). The influence of surface mechanics on diffusion induced stresses within spherical nanoparticles. *Journal of Applied Physics*, **104**, 083521.

Cheng, Y.-T., & Verbrugge, M. W. (2009). Evolution of stress within a spherical insertion electrode particle under potentiostatic and galvanostatic operation. *Journal of Power Sources*, **190**, 453-460.

Cheng, Y.-T., & Verbrugge, M. W. (2010). Diffusion-induced stress, interfacial charge transfer, and criteria for avoiding crack initiation of electrode particles. *Journal of the Electrochemical Society*, **157**, A508-A516.

Christensen, J., & Newman, J. (2006a). A mathematical model of stress generation and fracture in lithium manganese oxide. *Journal of The Electrochemical Society*, **153**, A1019-A1030.

Christensen, J., & Newman, J. (2006b). Stress generation and fracture in lithium insertion materials. *Journal of Solid State Electrochemistry*, **10**, 293-319.

Crank, J. (1979). *The mathematics of diffusion*. Oxford university press.

Deshpande, R., Cheng, Y.-T., & Verbrugge, M. W. (2010a). Modeling diffusion-induced stress in nanowire electrode structures. *Journal of Power Sources*, **195**, 5081-5088.

Deshpande, R., Cheng, Y.-T., Verbrugge, M. W., & Timmons, A. (2011). Diffusion induced stresses and strain energy in a phase-transforming spherical electrode particle. *Journal of the Electrochemical Society*, **158**, A718-A724.

Deshpande, R., Qi, Y., & Cheng, Y.-T. (2010b). Effects of concentration-dependent elastic modulus on diffusion-induced stresses for battery applications. *Journal of the Electrochemical Society*, **157**, A967-A971.

Esmizadeh, S., Haftbaradaran, H., & Mossaiby, F. (2019). An investigation of the critical conditions leading to deintercalation induced fracture in two-phase elastic electrode particles using a moving interphase core-shell model. *European Journal of Mechanics-A/Solids*, **74**, 96-111.
Golmon, S., Maute, K., Lee, S.-H., & Dunn, M. L. (2010). Stress generation in silicon particles during lithium insertion. *Applied Physics Letters*, **97**, 033111.

Hao, F., & Fang, D. (2013a). Diffusion-induced stresses of spherical core-shell electrodes in lithium-ion batteries: the effects of the shell and surface/interface stress. *Journal of The Electrochemical Society*, **160**, A595–A600.

Hao, F., & Fang, D. (2013b). Tailoring diffusion-induced stresses of core-shell nanotube electrodes in lithium-ion batteries. *Journal of Applied Physics*, **113**, 013507.

He, Y., Yu, X., Li, G., Wang, R., Li, H., Wang, Y., Gao, H., & Huang, X. (2012). Shape evolution of patterned amorphous and polycrystalline silicon microarray thin film electrodes caused by lithium insertion and extraction. *Journal of Power Sources*, **216**, 131–138.

Hertzberg, B., Benson, J., & Yushin, G. (2011). Ex-situ depth-sensing indentation measurements of electrochemically produced si–li alloy films. *Electrochemistry Communications*, **13**, 818–821.

Hu, P., Peng, W., Wang, B., Xiao, D., Ahuja, U., Réthoré, J., & Aifantis, K. E. (2019a). Concentration-gradient prussian blue cathodes for na-ion batteries. *ACS Energy Letters*, **5**, 100–108.

Hu, P., Wang, B., Xiao, D., & Aifantis, K. (2019b). Capturing the differences between lithiation and sodiation of nanostructured tis2 electrodes. *Nano Energy*, **63**, 103820.

Hu, Y., Zhao, X., & Suo, Z. (2010). Averting cracks caused by insertion reaction in lithium-ion batteries. *Journal of Materials Research*, **25**, 1007–1010.

Huang, S., Fan, F., Li, J., Zhang, S., & Zhu, T. (2013). Stress generation during lithiation of high-capacity electrode particles in lithium ion batteries. *Acta materialia*, **61**, 4354–4364.
Huttin, M., & Kamlah, M. (2012). Phase-field modeling of stress generation in electrode particles of lithium ion batteries. *Applied Physics Letters*, **101**, 133902.

630 Jaeger, J. C., & Carslaw, H. S. (1959). *Conduction of heat in solids*. Clarendon P.

Jagannathan, M., & Chandran, K. R. (2014). Analytical modeling and simulation of electrochemical charge/discharge behavior of si thin film negative electrodes in li-ion cells. *Journal of Power Sources*, **247**, 667–675.

635 Ko, S.-C., Lee, S., & Chou, Y. (2005). Chemical stresses in a square sandwich composite. *Materials Science and Engineering: A*, **409**, 145–152.

Larche, F., & Cahn, J. W. (1984). The interactions of composition and stress in crystalline solids. *J. Res. Natl. Bur. Stand.*, **89**, 467–500.

Larcht'e, F., & Cahn, J. (1982). The effect of self-stress on diffusion in solids. *Acta Metallurgica*, **30**, 1835–1845.

Lee, S., Wang, W., & Chen, J. (2000). Diffusion-induced stresses in a hollow cylinder. *Materials Science and Engineering: A*, **285**, 186–194.

Li, J. C.-M. (1978). Physical chemistry of some microstructural phenomena. *Metallurgical Transactions A*, **9**, 1353–1380.

Liu, H., Wu, Y., Rahm, E., Holze, R., & Wu, H. (2004). Cathode materials for lithium ion batteries prepared by sol-gel methods. *Journal of Solid State Electrochemistry*, **8**, 450–466.

Liu, X. H., Zhong, L., Huang, S., Mao, S. X., Zhu, T., & Huang, J. Y. (2012). Size-dependent fracture of silicon nanoparticles during lithiation. *ACS nano*, **6**, 1522–1531.

Liu, Y., Lv, P., Ma, J., Bai, R., & Duan, H. L. (2014). Stress fields in hollow core–shell spherical electrodes of lithium ion batteries. *Proceedings of*
the Royal Society A: Mathematical, Physical and Engineering Sciences, 470,
20140299.

Maxisch, T., & Ceder, G. (2006). Elastic properties of olivine li x fepo 4 from first principles. *Physical Review B*, 73, 174112.

McDowell, M. T., Ryu, I., Lee, S. W., Wang, C., Nix, W. D., & Cui, Y. (2012). Studying the kinetics of crystalline silicon nanoparticle lithiation with in situ transmission electron microscopy. *Advanced Materials*, 24, 6034-6041.

Natarajan, S., & Aifantis, K. (2020). Size effects in nanostructured li-ion battery cathode particles. *Journal of the Mechanical Behavior of Materials*, 29, 36-43.

Prussin, S. (1961). Generation and distribution of dislocations by solute diffusion. *Journal of Applied Physics*, 32, 1876-1881.

Qi, Y., Guo, H., Hector, L. G., & Timmons, A. (2010). Threefold increase in the young's modulus of graphite negative electrode during lithium intercalation. *Journal of The Electrochemical Society*, 157, A558-A566.

Qi, Y., Hector, L. G., James, C., & Kim, K. J. (2014). Lithium concentration dependent elastic properties of battery electrode materials from first principles calculations. *Journal of The Electrochemical Society*, 161, F3010-F3018.

Ryu, I., Choi, J. W., Cui, Y., & Nix, W. D. (2011). Size-dependent fracture of si nanowire battery anodes. *Journal of the Mechanics and Physics of Solids*, 59, 1717-1730.

Sethuraman, V. A., Chon, M. J., Shimshak, M., Van Winkle, N., & Guduru, P. R. (2010). In situ measurement of biaxial modulus of si anode for li-ion batteries. *Electrochemistry Communications*, 12, 1614-1617.

Shenoy, V. B., Johari, P., & Qi, Y. (2010). Elastic softening of amorphous and crystalline li-si phases with increasing li concentration: a first-principles study. *Journal of Power Sources*, 195, 6825-6830.
Shin, H.-C., & Pyun, S.-I. (1999a). An investigation of the electrochemical intercalation of lithium into a li1- δcoo2 electrode based upon numerical analysis of potentiostatic current transients. *Electrochimica acta*, **44**, 2235–2244.

Shin, H.-C., & Pyun, S.-I. (1999b). The kinetics of lithium transport through li1- δcoo2 by theoretical analysis of current transient. *Electrochimica acta*, **45**, 489–501.

Song, Y., Lu, B., Ji, X., & Zhang, J. (2012). Diffusion induced stresses in cylindrical lithium-ion batteries: analytical solutions and design insights. *Journal of The Electrochemical Society*, **159**, A2060–A2068.

Soni, S. K., Sheldon, B. W., Xiao, X., Bower, A. F., & Verbrugge, M. W. (2012). Diffusion mediated lithiation stresses in si thin film electrodes. *Journal of the Electrochemical Society*, **159**, A1520.

Srinivasan, V., & Newman, J. (2004). Discharge model for the lithium ironphosphate electrode. *Journal of the Electrochemical Society*, **151**, A1517–A1529.

Stournara, M. E., Guduru, P. R., & Shenoy, V. B. (2012). Elastic behavior of crystalline li-sn phases with increasing li concentration. *Journal of Power Sources*, **208**, 165–169.

Tsagrakis, I., & Aifantis, E. C. (2018). Gradient elasticity effects on the two-phase lithiation of lib anodes. In *Generalized Models and Non-classical Approaches in Complex Materials 2* (pp. 221–235). Springer.

Venkatraman, S., Subramanian, V., Kumar, S. G., Renganathan, N., & Muniyandi, N. (2000). Capacity of layered cathode materials for lithium-ion batteries—a theoretical study and experimental evaluation. *Electrochemistry communications*, **2**, 18–22.

Verbrugge, M. W., & Cheng, Y.-T. (2009). Stress and strain-energy distributions within diffusion-controlled insertion-electrode particles subjected to periodic
potential excitations. *Journal of The Electrochemical Society*, **156**, A927-A937.

Wang, D., Wu, X., Wang, Z., & Chen, L. (2005). Cracking causing cyclic instability of lifepo4 cathode material. *Journal of Power Sources*, **140**, 125-128.

Wang, J. W., He, Y., Fan, F., Liu, X. H., Xia, S., Liu, Y., Harris, C. T., Li, H., Huang, J. Y., Mao, S. X. et al. (2013). Two-phase electrochemical lithiation in amorphous silicon. *Nano letters*, **13**, 709-715.

Wang, W., Lee, S., & Chen, J. (2002). Effect of chemical stress on diffusion in a hollow cylinder. *Journal of applied physics*, **91**, 9584-9590.

Woodford, W. H., Chiang, Y.-M., & Carter, W. C. (2010). “electrochemical shock” of intercalation electrodes: a fracture mechanics analysis. *Journal of The Electrochemical Society*, **157**, A1052–A1059.

Yamada, A., Koizumi, H., Sonoyama, N., & Kanno, R. (2005). Phase change in li x fepo4. *Electrochemical and Solid-State Letters*, **8**, A409–A413.

Zhang, J., Lu, B., Song, Y., & Ji, X. (2012). Diffusion induced stress in layered li-ion battery electrode plates. *Journal of Power Sources*, **209**, 220-227.

Zhang, Q., & White, R. E. (2007). Moving boundary model for the discharge of a licoo2 electrode. *Journal of The Electrochemical Society*, **154**, A587–A596.

Zhang, T., & Kamlah, M. (2018). Sodium ion batteries particles: Phase-field modeling with coupling of cahn-hilliard equation and finite deformation elasticity. *Journal of The Electrochemical Society*, **165**, A1997–A2007.

Zhao, K., Pharr, M., Cai, S., Vlassak, J. J., & Suo, Z. (2011a). Large plastic deformation in high-capacity lithium-ion batteries caused by charge and discharge. *Journal of the American Ceramic Society*, **94**, s226–s235.

Zhao, K., Wang, W. L., Gregoire, J., Pharr, M., Suo, Z., Vlassak, J. J., & Kaxiras, E. (2011b). Lithium-assisted plastic deformation of silicon electrodes
in lithium-ion batteries: a first-principles theoretical study. *Nano letters*, **11**, 2962–2967.

Zuev, A. Y., & Tsvetkov, D. S. (2017). Conventional methods for measurements of chemo-mechanical coupling. In *Electro-chemo-mechanics of Solids* (pp. 5–33). Springer.