Comparison of Different Elastic Strain Definitions for Largely Deformed SEI of Chemo-Mechanically Coupled Silicon Battery Particles


Amorphous silicon is a highly promising anode material for next-generation lithium-ion batteries. Large volume changes of the silicon particle have a critical effect on the surrounding solid-electrolyte interphase (SEI) due to repeated fracture and healing during cycling. Based on a thermodynamically consistent chemo-elasto-plastic continuum model we investigate the stress development inside the particle and the SEI. Using the example of a particle with SEI, we apply a higher order finite element method together with a variable-step, variable-order time integration scheme on a nonlinear system of partial differential equations. Starting from a single silicon particle setting, the surrounding SEI is added in a first step with the typically used elastic Green–St-Venant (GSV) strain definition for a purely elastic deformation. For this type of deformation, the definition of the elastic strain is crucial to get reasonable simulation results. In case of the elastic GSV strain, the simulation aborts. We overcome the simulation failure by using the definition of the logarithmic Hencky strain. However, the particle remains unaffected by the elastic strain definitions in the particle domain. Compared to GSV, plastic deformation with the Hencky strain is straightforward to take into account. For the plastic SEI deformation, a rate-independent and a rate-dependent plastic deformation are newly introduced and numerically compared for three half cycles for the example of a radial symmetric particle.

finite elements , numerical simulation , solid-electrolyte interphase (SEI) , finite elastic strain , (visco-)plasticity 74C15 , 74C20 , 74S05 , 65M22 , 90C33

1 Introduction↩︎

Nowadays, the advancement of lithium-ion batteries is one key research in terms of climate change [1]. Alongside all solid-state batteries as promising candidates for future battery types [2], especially, amorphous silicon (aSi) is superior in terms of the nearly tenfold theoretical energy density compared to state of the art graphite anodes for classical battery types [3]. Unfortunately, this capacity increase is accompanied by a volume change up to 300% crucially effecting the surrounding solid-electrolyte interphase (SEI) during swelling and shrinking [3], [4]. The SEI features elastic and plastic deformation before it can break up and heal again during repeated cycling [4]. In addition, the voltage hysteresis of silicon is still an ongoing research topic because it influences battery lifetime and performance and makes further mechanical investigation of the SEI unavoidable [5]. Recent measurements of silicon anodes show large overpotentials with slow relaxation which could confirm the mechanical stress influence of the SEI and could explain the voltage hysteresis [6].

In this paper, we apply modern numerical techniques on the extended underlying system of equations, based on an efficient adaptive algorithm, compare [7] and [8]. Furthermore, we discuss the influence of the definition of the elastic strain tensor for the SEI domain on the numerical method, starting with a single particle setting [7]. Due to the large volume change, a finite deformation approach is needed for the particle as well as for the SEI domain. Firstly, a purely elastic case for the SEI is considered with two different definitions for the elastic strain tensor. We use the definition of the (Lagrangian) Green–St-Venant (GSV) strain tensor and the logarithmic (Hencky) strain tensor. In the situation of a single particle, both definitions lead to quantitative similar numerical results [9], [10]. The strain definition plays a significant role for the SEI domain whether the numerical simulation will abort or not. To the best of the authors’ knowledge, this has not yet been documented in the literature. Using the logarithmic strain approach, we can comfortably add a rate-independent and a rate-dependent plastic deformation for the SEI without increasing the system of equation for the plastic part of the deformation gradient [4], [9]. The rate-dependent case results is a typical stress-overrelaxation, which is also found in measurements of the electric field [6] and can confirm the hypothesis of [5].

The rest of the article is structured as follows: in the next section, we present and recap our used model approach and present the different definitions for the elastic strain approaches. In 3, we formulate our model equations and show important steps towards a numerical solution. 4 shows the failure of one of the used definitions of the elastic strains and successful numerical simulations for the other definition. In addition, the latter one is extended with plastic and viscoplastic deformation. In the end, we sum up our main findings and close with an outlook.

2 Theory↩︎

In this section, we recall and summarize our used model equations for the chemo-mechanically coupled particle-SEI approach based on the thermodynamically consistent theory by [4], [7], [9]. For the electrode particle, we use a pure elastic deformation model with a Green–St-Venant (GSV) strain tensor [4]. For the SEI, we introduce four different approaches for the deformation: 1.) a pure elastic with the GSV strain tensor, 2.) a pure elastic (Lagrangian) logarithmic Hencky strain as well as 3.) a plastic and 4.) a viscoplastic deformation approach, whereby last two use the logarithmic Hencky strain [9].

Finite Deformation. Let \(\Omega_0\subset \mathbb{R}^3\) be a bounded domain which represents the particle-SEI domain in the reference (Lagrangian) configuration, divided into a particle domain \(\Omega_{0,\text{P}}\) and a SEI domain \(\Omega_{0,\text{S}}\). We denote the current (Eulerian) domain configuration without the subscript \(0\). With the displacement \(\boldsymbol{u}{}\) and the deformation \(\Phi\colon \mathbb{R}_{\geq 0}\times \Omega_0\rightarrow \Omega\), \(\Phi(t, \boldsymbol{X_0}{}) \mathrel{\vcenter{:}}= \boldsymbol{x}{} = \boldsymbol{X}{}_0 + \boldsymbol{u}{}(t, \boldsymbol{X}{}_0)\) from the reference configuration to the current configuration [11], [12], [13], [4], [7], [9], [14], [15], we can derive the deformation gradient \(\mathbf{F}{} (\boldsymbol{\nabla}_0\boldsymbol{u}{})\in \mathbb{R}^{3,3}\) as \(\mathbf{F}{} \mathrel{\vcenter{:}}= \partial\Phi / \partial\boldsymbol{X}{}_0 = \mathbf{Id}{} + \boldsymbol{\nabla}_0\boldsymbol{u}{}\) with the identity tensor \(\mathbf{Id}{}\). Following [16], [12] and [7], [15], we multiplicatively split the deformation gradient \(\mathbf{F}{} = \mathbf{F}{}_\text{ch}\mathbf{F}{}_\text{el}\mathbf{F}{}_\text{pl}\) into three parts: chemical, elastic and plastic part, respectively. Note that we have for each domain \(\Omega_{0,\text{P}}\) and \(\Omega_{0,\text{S}}\) own deformation gradients \(\mathbf{F}{}_\text{P}\) and \(\mathbf{F}{}_\text{S}\). We consider only chemical and elastic deformation for the particle domain (\(\mathbf{F}{}_{\text{pl}, \text{P}} = \mathbf{Id}{}\)), whereas only elastic and plastic deformations are considered for the SEI domain (no calculation of any chemical species, \(\mathbf{F}{}_{\text{ch}, \text{S}} = \mathbf{Id}{}\)). Furthermore, we identify two displacements: the particle displacement \(\boldsymbol{u}{}_\text{P} \colon \mkern 1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern 1.5mu_{0, \text{P}, t_\text{end}} \rightarrow \mathbb{R}^3\) and the SEI displacement \(\boldsymbol{u}{}_\text{S} \colon \mkern 1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern 1.5mu_{0, \text{S}, t_\text{end}} \rightarrow \mathbb{R}^3\) with the final simulation time \(t_\text{end}> 0\) and \(\mkern 1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern 1.5mu_{0, \square, t_\text{end}} \mathrel{\vcenter{:}}= [0, t_\text{end}] \times \mkern 1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern 1.5mu_{0, \square}\), \(\square \in \{ \text{P}, \text{S}\}\). The chemical part arises due to lithium insertion in the particle domain and is stated as \(\mathbf{F}{}_{\text{ch}, \text{P}} (\mkern 1.5mu\overline{\mkern-1.5muc\mkern-1.5mu}\mkern 1.5mu_\text{P}(t, \boldsymbol{X}{}_{0, \text{P}})) = \mathbf{F}{}_\text{ch}(\mkern 1.5mu\overline{\mkern-1.5muc\mkern-1.5mu}\mkern 1.5mu) = \lambda_\text{ch}(\mkern 1.5mu\overline{\mkern-1.5muc\mkern-1.5mu}\mkern 1.5mu) \mathbf{Id}{} = \sqrt[3]{1+ v_\text{pmv} c_{\text{max}}\mkern 1.5mu\overline{\mkern-1.5muc\mkern-1.5mu}\mkern 1.5mu} \mathbf{Id}{}\) with the partial molar volume \(v_\text{pmv}> 0\) of lithium inside aSi and the normalized lithium concentration \(\mkern 1.5mu\overline{\mkern-1.5muc\mkern-1.5mu}\mkern 1.5mu= c / c_{\text{max}}\in [0, 1]\) of the lithium concentration \(c \colon\) \(\mkern 1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern 1.5mu_{0, \text{P}, t_\text{end}}\) \(\rightarrow\) \([0, c_{\text{max}}]\) with respect to the maximal concentration of aSi, \(c_{\text{max}}\), in the reference configuration. Using the chemical part \(\mathbf{F}{}_\text{ch}\) and the plastic part \(\mathbf{F}{}_\text{pl}(t, \boldsymbol{X}{}_0)\), applied as internal variable [9], the elastic part can be computed as \(\mathbf{F}{}_{\text{el}, \text{P}} (\mkern 1.5mu\overline{\mkern-1.5muc\mkern-1.5mu}\mkern 1.5mu, \boldsymbol{\nabla}_0\boldsymbol{u}{}_\text{P}) = \mathbf{F}{}_{\text{el}} (\mkern 1.5mu\overline{\mkern-1.5muc\mkern-1.5mu}\mkern 1.5mu, \boldsymbol{\nabla}_0\boldsymbol{u}{}_\text{P}) = \lambda_\text{ch}^{{\scalebox{0.75}[1.0]{-}1}} \mathbf{F}{}\) and \(\mathbf{F}{}_{\text{el}, \text{S}} (\boldsymbol{\nabla}_0\boldsymbol{u}{}_\text{S}) = \mathbf{F}{}_{\text{el}} (\boldsymbol{\nabla}_0\boldsymbol{u}{}_\text{S}) = \mathbf{F}{} \mathbf{F}{}_\text{pl}^{{\scalebox{0.75}[1.0]{-}1}}\) with the gradient \(\boldsymbol{\nabla}_0\) in the respective domain. If it is clear from the context which domain is considered, the index \(\text{P}\) or \(\text{S}\) is omitted for reasons of better readability.

Free Energy. Following the thermodynamically consistent material approach [4], [9], we additively decompose the Helmholtz free energy \(\psi\) into a chemical part \(\psi_\text{ch}\) and a mechanical part \(\psi_\text{el}\) at constant temperature in the reference configuration: \[\begin{align} \psi(\mkern 1.5mu\overline{\mkern-1.5muc\mkern-1.5mu}\mkern 1.5mu, \boldsymbol{\nabla}_0\boldsymbol{u}{}, \mathbf{F}{}_\text{pl}) = \psi_\text{ch}(\mkern 1.5mu\overline{\mkern-1.5muc\mkern-1.5mu}\mkern 1.5mu) + \psi_\text{el}(\mkern 1.5mu\overline{\mkern-1.5muc\mkern-1.5mu}\mkern 1.5mu, \boldsymbol{\nabla}_0\boldsymbol{u}{}, \mathbf{F}{}_\text{pl}). \end{align}\] For the respective domains, we have \(\psi_\text{P}(\mkern 1.5mu\overline{\mkern-1.5muc\mkern-1.5mu}\mkern 1.5mu, \boldsymbol{\nabla}_0\boldsymbol{u}{}_\text{P}) = \psi_\text{ch}(\mkern 1.5mu\overline{\mkern-1.5muc\mkern-1.5mu}\mkern 1.5mu) + \psi_{\text{el}, \text{P}}(\mkern 1.5mu\overline{\mkern-1.5muc\mkern-1.5mu}\mkern 1.5mu, \boldsymbol{\nabla}_0\boldsymbol{u}{}_\text{P})\) and \(\psi_\text{S}(\boldsymbol{\nabla}_0\boldsymbol{u}{}_\text{S}, \mathbf{F}{}_\text{pl}) = \psi_{\text{el}, \text{S}}(\boldsymbol{\nabla}_0\boldsymbol{u}{}_\text{S}, \mathbf{F}{}_\text{pl})\). The chemical free energy density is defined by an experimental open-circuit voltage (OCV) curve \(U_\text{OCV}\) [4], [9], [14] \[\begin{align} \rho_0\psi_\text{ch}(\mkern 1.5mu\overline{\mkern-1.5muc\mkern-1.5mu}\mkern 1.5mu) = \scalebox{0.75}[1.0]{-} c_{\text{max}} \int_0^{\mkern 1.5mu\overline{\mkern-1.5muc\mkern-1.5mu}\mkern 1.5mu} \mathrm{Fa} \, U_\text{OCV}(z) \,\mathrm{d}z \end{align}\] with the mass density \(\rho_0\) of aSi in the reference configuration and the Faraday constant \(\mathrm{Fa}\). The mechanical free energy density is stated for both the particle and the SEI domain by a linear elastic approach [4], [7], [9] \[\begin{align} \rho_0 \psi_\text{el}(\mkern 1.5mu\overline{\mkern-1.5muc\mkern-1.5mu}\mkern 1.5mu, \boldsymbol{\nabla}_0\boldsymbol{u}{}, \mathbf{F}{}_\text{pl}) = \frac{1}{2}\mathbf{E}{}_\text{el}(\mkern 1.5mu\overline{\mkern-1.5muc\mkern-1.5mu}\mkern 1.5mu, \boldsymbol{\nabla}_0\boldsymbol{u}{}, \mathbf{F}{}_\text{pl}) \!:\! \mathbb{C}[\mathbf{E}{}_\text{el}] \end{align}\] with the elastic strain tensor \(\mathbf{E}{}_\text{el}\), the constant, isotropic stiffness tensor \(\mathbb{C}\) of aSi and \(\mathbb{C}[\mathbf{E}{}_\text{el}] = \lambda \mathop{\mathrm{tr \!}}\left(\mathbf{E}{}_\text{el}\right) \mathbf{Id}{} + 2G \mathbf{E}{}_\text{el}\). The first and second Lamé constants \(\lambda = 2G \nu / (1-2\nu)\) and \(G = E / \big( 2(1+\nu)\big)\) depend on the Young’s modulus \(E\) and the Poisson’s ratio \(\nu\). The parameters for the particle and the SEI material are specified in 4.

Our main concern in this paper is the definition of the elastic strain tensor. In [9], two different definitions for the elastic strain tensor are stated: the GSV strain tensor, also called the Lagrangian strain tensor [12] \[\begin{align} \mathbf{E}{}_{\text{el}, \text{lag}} = \mathbf{E}{}_{\text{el}, \text{GSV}} = \frac{1}{2}\left(\mathbf{C}{}_\text{el}- \mathbf{Id}{}\right) = \frac{1}{2}\left(\mathbf{F}{}_\text{el}^\textsf{T}\mathbf{F}{}_\text{el}- \mathbf{Id}{}\right), \end{align}\] and the (Lagrangian) logarithmic Hencky strain tensor \[\begin{align} \mathbf{E}{}_{\text{el}, \text{log}} = \ln \left(\mathbf{U}{}_\text{el}\right) = \ln \left(\sqrt{\mathbf{C}{}_\text{el}} \right) = \sum_{\alpha = 1}^{3} \ln \left(\sqrt{\eta_{\text{el}, \alpha}}\right) \boldsymbol{r}{}_{\text{el}, \alpha} \otimes \boldsymbol{r}{}_{\text{el}, \alpha} \end{align}\] with the eigenvalues \(\eta_{\text{el}, \alpha}\) and eigenvectors \(\boldsymbol{r}{}_{\text{el}, \alpha}\) of \(\mathbf{U}{}_\text{el}\) being the unique, positive definite and symmetric right stretch part of the unique polar decomposition of \(\mathbf{F}{}_\text{el}= \mathbf{R}{}_\text{el}\mathbf{U}{}_\text{el}\) [11]. The right Cauchy–Green tensor \(\mathbf{C}{} = \mathbf{F}{}^\textsf{T}\mathbf{F}{}\) is also symmetric and positive definite [12]. For the situation of considering an elastic particle without SEI only, both elastic strain approach results have similar outcomes [9], [10]. In this paper, we use for the particle domain the GSV strain. For the SEI domain, we will apply both strain approaches, since especially the logarithmic strain is favorable to incorporate plastic deformation [9].

Chemistry. A generalized diffusivity equation [4], [17], [18] is used to describe the change of lithium concentration inside the reference particle domain \(\Omega_{0,\text{P}}\) \[\begin{align} \partial_t c = \scalebox{0.75}[1.0]{-}\boldsymbol{\nabla}_0 \!\cdot\!\boldsymbol{N}{} \end{align}\] with the lithium flux \(\boldsymbol{N}{}(\mkern 1.5mu\overline{\mkern-1.5muc\mkern-1.5mu}\mkern 1.5mu, \boldsymbol{\nabla}_0\boldsymbol{u}{}_\text{P}) = \scalebox{0.75}[1.0]{-}m (\mkern 1.5mu\overline{\mkern-1.5muc\mkern-1.5mu}\mkern 1.5mu, \boldsymbol{\nabla}_0\boldsymbol{u}{}_\text{P}) \boldsymbol{\nabla}_0\mu = \scalebox{0.75}[1.0]{-}D \left(\partial_c \mu \right)^{{\scalebox{0.75}[1.0]{-}1}} \boldsymbol{\nabla}_0\mu\), the scalar mobility \(m > 0\) for the applied isotropic case, the diffusion coefficient \(D > 0\) for lithium atoms in aSi. The chemical potential \(\mu \colon \mkern 1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern 1.5mu_{0, \text{P},t_\text{end}} \rightarrow \mathbb{R}\) is given as the partial derivative of the free energy density with respect to \(c\) [14], [19]: \[\begin{align} \mu & = \partial_c (\rho_0 \psi) = \scalebox{0.75}[1.0]{-}\mathrm{Fa} \, U_\text{OCV} - \frac{v_\text{pmv}}{3 \lambda_\text{ch}^5} \mathbf{F}{}^\textsf{T}\mathbf{F}{} \!:\! \mathbb{C}[\mathbf{E}{}_\text{el}]. \end{align}\]

At the surface of the particle domain, we apply a uniform and constant external lithium flux \(N_\text{ext}\). The sign of this flux is positive or negative, depending on lithium insertion and extraction, respectively, and is expressed in terms of the charging rate (C-rate). The state of charge (SOC) links the simulation time, the external lithium flux and the initial concentration via \(\text{SOC}(t) = \mkern 1.5mu\overline{\mkern-1.5muc\mkern-1.5mu}\mkern 1.5mu_0 + N_\text{ext}t\). For more information about \(N_\text{ext}\), the C-rate and the SOC, we refer to [4], [7], [9] and the references cited therein.

Elastic and Inelastic Deformation. In both domains, the static balance of linear of momentum is used to consider the deformation [4], [7], [9], [14]: \[\begin{align} \boldsymbol{0}{} = \scalebox{0.75}[1.0]{-} \boldsymbol{\nabla}_0 \!\cdot\!\mathbf{P}{}_\text{P}(\mkern 1.5mu\overline{\mkern-1.5muc\mkern-1.5mu}\mkern 1.5mu, \boldsymbol{\nabla}_0\boldsymbol{u}{}_\text{P}), \qquad \boldsymbol{0}{} = \scalebox{0.75}[1.0]{-} \boldsymbol{\nabla}_0 \!\cdot\!\mathbf{P}{}_\text{S}(\boldsymbol{\nabla}_0\boldsymbol{u}{}_\text{S}, \mathbf{F}{}_\text{pl}) \end{align}\] with the first Piola–Kirchhoff tensor \(\mathbf{P}{} \in \mathbb{R}^{3,3}\), thermodynamically consistent derived as \(\mathbf{P}{} = \partial_\mathbf{F}{}(\rho_0 \psi)\). This results in \(\mathbf{P}{}_{\text{P}, \text{lag}} = \lambda_\text{ch}^{-2} \mathbf{F}{} \mathbb{C}_\text{P}[\mathbf{E}{}_{\text{el}, \text{lag}}]\) for the particle domain and in \(\mathbf{P}{}_{\text{S}, \text{lag}} = \mathbf{F}{} \left(\mathbf{F}{}_\text{pl}^{{\scalebox{0.75}[1.0]{-}1}} \right)^{\textsf{T}}\mathbf{F}{}_\text{pl}^{{\scalebox{0.75}[1.0]{-}1}} \mathbb{C}[\mathbf{E}{}_{\text{el}, \text{lag}}]\) or \(\mathbf{P}{}_{\text{S}, \text{log}} = \mathbf{F}{} \left(\mathbf{F}{}_\text{el}^\textsf{T}\mathbf{F}{}_\text{el}^{\vphantom{\textsf{T}}}\right)^{{\scalebox{0.75}[1.0]{-}1}} \left(\mathbf{F}{}_\text{pl}^{{\scalebox{0.75}[1.0]{-}1}} \right)^{\textsf{T}}\mathbf{F}{}_\text{pl}^{{\scalebox{0.75}[1.0]{-}1}} \mathbb{C}[\mathbf{E}{}_{\text{el}, \text{log}}]\) for the SEI domain depending on the definition of the GSV strain tensor or the logarithmic strain, respectively [4], [5], [9], [14]. The Cauchy stress \(\boldsymbol{\sigma} \in \mathbb{R}_{\text{sym}}^{3,3}\) in the current configuration is defined as \(\boldsymbol{\sigma} = \mathbf{P}{} \mathbf{F}{}^\textsf{T}/ \det \left(\mathbf{F}{}\right)\) [11] with \(\det \left(\mathbf{F}{}\right) > 0\) [11]. For the inelastic deformation approach, we base on the theory in [9]. There are a rate-independent plastic approach with isotropic hardening and a rate-dependent plastic approach developed and compared. The special feature is the usage of a projector formulation mapping the stresses onto the set of admissible stresses [20], a concept also known as static condensation [21], [22]. Following [9], we define for the SEI domain \(\Omega_{0,\text{S}}\) the Mandel stress \(\mathbf{M}{}_\text{S}(\boldsymbol{\nabla}_0\boldsymbol{u}{}_\text{S}, \mathbf{F}{}_\text{pl}) = \mathbf{M}{} =\partial_{\mathbf{E}{}_\text{el}} (\rho_{0, \text{S}} \psi_{\text{el}, \text{S}}) = \mathbb{C}_\text{S}[\mathbf{E}{}_{\text{el}}]\). For the rate-independent ideal plasticity (\(\gamma^\text{iso} = 0\) in [9], [20]), we express the classical loading and unloading conditions via the Karush–Kuhn–Tucker (KKT) [23], [12], [4] as \[\begin{align} F_\text{Y} \leq 0, \qquad \dot{\varepsilon}_\text{pl}^\text{eq} \geq 0, \qquad F_\text{Y} \dot{\varepsilon}_\text{pl}^\text{eq} = 0 \end{align}\] with the yield function \(F_\text{Y} (\boldsymbol{\nabla}_0\boldsymbol{u}{}_\text{S}, \mathbf{F}{}_\text{pl}, \varepsilon_\text{pl}^\text{eq}) = \|\mathbf{M}{}^\text{dev}\| - \sigma_\text{Y}\), the deviatoric Mandel stress \(\mathbf{M}{}^\text{dev} = \mathbf{M}{} - 1/3 \mathop{\mathrm{tr \!}}\left( \mathbf{M}{} \right) \mathbf{Id}{}\), the yield stress \(\sigma_\text{Y}\) and the accumulated equivalent inelastic strain \({\varepsilon}_\text{pl}^\text{eq}(t, \boldsymbol{X}{}_{0, \text{S}}) \geq 0\). To be consistent with the one dimensional tensile test, we rescale the yield stress with the factor \(\sqrt{2/3}\) [23]. In the viscoplastic approach, an evolution equation of the equivalent plastic strain, given by \[\tag{1} \begin{empheq}[ left={ \dot{\varepsilon}_\text{pl}^\text{eq} = \empheqlbrace} ]{alignat=3} &0,\tag{2} && \big|\!\big| \mathbf{M}{}^{\text{dev}} \big|\!\big|_{}^{} \leq \sigma_\text{Y}, \\ & \dot{\varepsilon}_0 \Bigg( \frac{\big|\!\big| \mathbf{M}{}^{\text{dev}} \big|\!\big|_{}^{} - \sigma_\text{Y}}{\sigma_{\text{Y}^{*}}}\Bigg)^\beta, \quad && \tag{3} \big|\!\big| \mathbf{M}{}^{\text{dev}} \big|\!\big|_{}^{} > \sigma_\text{Y}, \end{empheq}\] describes the plastic deformation instead of the KKT conditions [23], [9], [15]. The positive-valued stress-dimensioned constant \(\sigma_{\text{Y}^{*}}\), the reference tensile stress \(\dot{\varepsilon}_0\) and the measure of the strain rate sensitivity of the material \(\beta\) are defined in 4.

3 Numerical Approach↩︎

3.1 Problem Formulation↩︎

Using Table 1 in [9], we arrive at the dimensionless initial boundary value problem with inequality conditions: find the normalized concentration \(c \colon \mkern 1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern 1.5mu_{0, \text{P}, t_\text{end}} \rightarrow [0,1]\), the chemical potential \(\mu \colon \mkern 1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern 1.5mu_{0, \text{P}, t_\text{end}} \rightarrow \mathbb{R}\), the displacements \(\boldsymbol{u}{}_\text{P} \colon \mkern 1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern 1.5mu_{0, \text{P}, t_\text{end}} \rightarrow \mathbb{R}^3\) and \(\boldsymbol{u}{}_\text{S} \colon \mkern 1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern 1.5mu_{0, \text{S}, t_\text{end}} \rightarrow \mathbb{R}^3\) satisfying \[\tag{4} \begin{empheq}[left=\empheqlbrace]{alignat=4} \partial_t c &= \scalebox{0.75}[1.0]{-}\boldsymbol{\nabla}_0 \!\cdot\!\boldsymbol{N}{} (c, \boldsymbol{\nabla}_0\boldsymbol{u}{}_\text{P}) && \qquad && \text{in }\Omega_{0, \text{P}, t_\text{end}}, \tag{5} \\ \mu &= \partial_c (\rho_0 \psi (c, \boldsymbol{\nabla}_0\boldsymbol{u}{}_\text{P})) && && \text{in }\Omega_{0, \text{P}, t_\text{end}}, \tag{6} \\ \mathbf{0}{} &= \boldsymbol{\nabla}_0 \!\cdot\!\mathbf{P}{}_\text{P} (c, \boldsymbol{\nabla}_0\boldsymbol{u}{}_\text{P}) && && \text{in }\Omega_{0, \text{P}, t_\text{end}}, \tag{7} \\ \mathbf{0}{} &= \boldsymbol{\nabla}_0 \!\cdot\!\mathbf{P}{}_\text{S} (\boldsymbol{\nabla}_0\boldsymbol{u}{}_\text{S}, \mathbf{F}{}_\text{pl}) && && \text{in }\Omega_{0, \text{S}, t_\text{end}}, \tag{8} \\ F_\text{Y} \leq 0,\quad& \dot{\varepsilon}_\text{pl}^\text{eq} \geq 0, \quad F_\text{Y} \dot{\varepsilon}_\text{pl}^\text{eq} = 0 && && \text{in }\Omega_{0, \text{S}, t_\text{end}}, \tag{9} \\ \boldsymbol{N}{} \cdot \boldsymbol{n}{}_\text{P} &= N_\text{ext} && && \text{on }\partial\Omega_{0, \text{P}, t_\text{end}} \tag{10} \\ \boldsymbol{u}{}_\text{P} &= \boldsymbol{u}{}_\text{S} && && \text{on }\Gamma_{\text{inter}, t_\text{end}} \tag{11} \\ \mathbf{P}{}_\text{P}\cdot \boldsymbol{n}{}_\text{P} &= \mathbf{P}{}_\text{S}\cdot \boldsymbol{n}{}_\text{P} && && \text{on }\Gamma_{\text{inter}, t_\text{end}} \tag{12} \\ \mathbf{P}{}_\text{S}\cdot \boldsymbol{n}{}_\text{S} &= \mathbf{0}{} && && \text{on }\partial\Omega_{0, \text{S}, t_\text{end}} \tag{13} \\ c(0, \boldsymbol{\cdot}) &= c_0 && && \text{in }\Omega_{0, \text{P}}, \tag{14} \\ \mathbf{F}{}_\text{pl}(0, \boldsymbol{\cdot}) &= \mathbf{Id}{} && && \text{in }\Omega_{0, \text{S}}, \tag{15} \\ \varepsilon_\text{pl}^{\text{eq}}(0, \boldsymbol{\cdot}) &= 0 && && \text{in }\Omega_{0, \text{S}} \tag{16} \end{empheq}\] with the interface of particle and SEI domain \(\Gamma_{\text{inter}, t_\text{end}} \mathrel{\vcenter{:}}= [0, t_\text{end}] \times \Gamma_{\text{inter}} = [0, t_\text{end}] \times (\partial\Omega_{0, \text{P}} \cap \partial\Omega_{0, \text{S}})\), the outward unit normal vector \(\boldsymbol{n}{}_\text{P}= \scalebox{0.75}[1.0]{-}\boldsymbol{n}{}_\text{S}\) and the external boundary \(\partial\Omega_{0, \text{S}, t_\text{end}} \mathrel{\vcenter{:}}= [0, t_\text{end}] \times \partial\Omega_{0, \text{S}}\). We assume a constant lithium flux with changing sign due to cycling at the particle boundary as well as equal displacements and stresses at the interface boundary \(\Gamma_{\text{inter}, t_\text{end}}\) as well as no stresses at the SEI boundary \(\partial\Omega_{0, \text{S}}\). The initial concentration \(c_0\) is chosen to be constant and rigid body motions are excluded with appropriate displacement boundary conditions. For the detailed treatment of the plastic part \(\mathbf{F}{}_\text{pl}\) of the deformation gradient and the equivalent plastic strain \(\varepsilon_\text{pl}^\text{eq}\) as internal variables, we refer to [9]. Finally, the desired quantities, the Cauchy stresses \(\boldsymbol{\sigma}_\text{P}\) and \(\boldsymbol{\sigma}_\text{P}\), are computed in a postprocessing step using the solution variables, \(\mathbf{F}{}\), \(\mathbf{F}{}_\text{el}\), \(\mathbf{E}{}_\text{el}\) and \(\mathbf{P}{}\).

3.2 Numerical Solution Procedure↩︎

Weak Formulation. Following [7], [9], [20], we multiply with test functions, integrate over the respective domain, integrate by parts and formulate a weak primal mixed variational inequality for 4 . Using a projector \(\mathbf{P}{}_\Pi\) onto the set of admissible stresses [20], [24], we can reformulate the occurring saddle point problem as a primal formulation: for given \(\mathbf{F}{}_\text{pl}\) and \(\varepsilon_\text{pl}^\text{eq}\) find solutions \(\lbrace c, \mu, \boldsymbol{u}{}_\text{P}, \boldsymbol{u}{}_\text{S}\rbrace\) with \(c, \mu \in V\), \(\partial_t c \in L^2(\Omega_{0}, \mathbb{R})\), \(\boldsymbol{u}{}_\text{P}\in \boldsymbol{V}{}_\text{P}^{*}\) and \(\boldsymbol{u}{}_\text{S}\in \boldsymbol{V}{}_\text{S}^{*}\) such that \[\tag{17} \begin{empheq}[left=\empheqlbrace]{alignat=2} \big( \varphi,\partial_t c \big) &= \scalebox{0.75}[1.0]{-} \big( m (c, \boldsymbol{\nabla}_0 \boldsymbol{u}{}_\text{P}) \boldsymbol{\nabla}_0\varphi, \boldsymbol{\nabla}_0\mu \big) - \big(\varphi, N_\text{ext}\big)_{\Gamma_{\text{inter}}},\tag{18} \\ 0 &= \scalebox{0.75}[1.0]{-} \big( \varphi, \mu \big) + \big( \varphi, \partial_c ( \rho_0 \psi_\text{ch}(c)) \nonumber \\ & \quad \quad \quad \quad \quad + \partial_c ( \rho_0 \psi_\text{el} (c, \boldsymbol{\nabla}_0 \boldsymbol{u}{})) \big),\tag{19} \\ \mathbf{0}{} & = \scalebox{0.75}[1.0]{-} \big( \boldsymbol{\nabla}_0\boldsymbol{\xi}{}, \mathbf{P}{}_\text{P} ( c, \boldsymbol{\nabla}_0\boldsymbol{u}{}_\text{P}) \big) + \big(\boldsymbol{\xi}{}, \mathbf{P}{}_\text{S} \cdot \boldsymbol{n}{}_\text{P} \big)_{\Gamma_{\text{inter}}} \tag{20} \\ \mathbf{0}{} & = \scalebox{0.75}[1.0]{-} \big( \boldsymbol{\nabla}_0\boldsymbol{\chi}{}, \mathbf{P}{}_\text{S} (\boldsymbol{\nabla}_0\boldsymbol{u}{}_\text{S}, \mathbf{F}{}_\text{pl}, \mathbf{P}{}_\Pi (\boldsymbol{\nabla}_0 \boldsymbol{u}{}_\text{S}, \mathbf{F}{}_\text{pl}, \varepsilon_\text{pl}^\text{eq}) ) \big) \nonumber \\ & \quad - \big(\boldsymbol{\chi}{}, \mathbf{P}{}_\text{P} \cdot \boldsymbol{n}{}_\text{P} \big)_{\Gamma_{\text{inter}}} \tag{21} \end{empheq}\] for all test functions \(\varphi \in V\), \(\boldsymbol{\xi}{} \in \boldsymbol{V}{}_\text{P}^{*}\), \(\boldsymbol{\chi}{} \in \boldsymbol{V}{}_\text{S}^{*}\) and \(\mathbf{P}{}_\Pi (\boldsymbol{\nabla}_0 \boldsymbol{u}{}_\text{S}, \mathbf{F}{}_\text{pl}, \varepsilon_\text{pl}^\text{eq}) = \mathbb{C}[\mathbf{E}{}_{\text{el}, \text{log}}]\) with appropriate function spaces \(V, \boldsymbol{V}{}_\text{P}^*\) and \(\boldsymbol{V}{}_\text{S}^*\) to guarantee a well-posed formulation [25]. The last two spaces include constraints to avoid rigid body motions. In 17 , \((\boldsymbol{\cdot}, \boldsymbol{\cdot})\) indicates the \(L^2\)-inner product for two scalar functions, vectors and tensors, respectively. The notation with boundary index \((\boldsymbol{\cdot}, \boldsymbol{\cdot})_{\Gamma_{\text{inter}}}\) denotes the interface integral at the interface \({\Gamma_{\text{inter}}}\) regarding \(\boldsymbol{n}{}_\text{P}\). 9 results in a saddle point problem, which is condensed into 21 with the respective projector formulation for the rate-independent or rate-dependent plastic approach [9].

Discretization. Considering an admissible mesh \(\mathcal{T}_h\) as discretization of the computational domain \(\Omega_\text{h}\), we apply the isoparametric Lagrangian finite element method [13] and insert spatial approximations \(c_h, \mu_h, \boldsymbol{u}{}_{\text{P}, h}\) and \(\boldsymbol{u}{}_{\text{S}, h}\) with the finite dimensional subspaces of the respective function space: \[\begin{alignat}{2} &V_h &&= {\text{span}} \{ \varphi_i \,:\,i=1,\dots, N_\text{P}\} \subset {V}, \\ &\boldsymbol{V}{}_{\text{P}, h}^{*} &&= {\text{span}} \{ \boldsymbol{\xi}{}_j \,:\,j=1,\dots, 3 N_\text{P}\} \subset \boldsymbol{V}{}_{\text{P}}^{*}, \\ &\boldsymbol{V}{}_{\text{S}, h}^{*} &&= {\text{span}} \{ \boldsymbol{\chi}{}_k \,:\,k=1,\dots, 3 N_\text{S}\} \subset \boldsymbol{V}{}_{\text{S}}^{*}. \end{alignat}\] For the temporal discretization, we collect all time-dependent coefficients in the vector valued function \[\begin{align} \label{eq:discrete95solution95vector} \boldsymbol{y}{} \colon[0, t_\text{end}] \rightarrow \mathbb{R}^{(2+3)N_\text{P}+ 3 N_\text{S}}, \quad t \mapsto \boldsymbol{y}{}(t) = \begin{pmatrix} \boldsymbol{c}{}_h(t) \\ \boldsymbol{\mu}{}_h(t) \\ \boldsymbol{u}{}_{\text{P},h}(t) \\ \boldsymbol{u}{}_{\text{S},h}(t) \end{pmatrix} \end{align}\tag{22}\] satisfying \(\mathbf{M}{} \partial_t \boldsymbol{y}{} - \boldsymbol{f}{}(t, \boldsymbol{y}{}, \mathbf{F}{}_{\text{pl}, h}, \varepsilon_{\text{pl}, h}^\text{eq}) = \mathbf{0}{}\) for \(t \in (0, t_\text{end}]\) with \(\boldsymbol{y}{}(0) = \boldsymbol{y}{}^0\) and the discrete version of the internal variables \(\mathbf{F}{}_{\text{pl}, h}, \varepsilon_{\text{pl}, h}^\text{eq}\) [9]. For the temporal discretization of the internal variables, we use an implicit exponential map [9]. In total, we get the space and time discrete problem to go on from one time \(t_n\) to the next \(t_{n+1}= t_n + \tau_n\) with \(\tau_n> 0\) applying the numerical differential formulation (NDF) of linear multistep methods [26][28]. Then, for given \(\mathbf{F}{}_{\text{pl},h}^n\) and \(\varepsilon_{\text{pl},h}^{\text{eq},n}\) find the discrete solution \(\boldsymbol{y}{}^{{n+1}} \approx \boldsymbol{y}{}(t_{{n+1}})\) satisfying \[\begin{align} \label{eq:space95and95time95discretization} \alpha_{k_n} \mathbf{M}{} \left(\boldsymbol{y}{}^{n+1}- \boldsymbol{\Phi}{}^n\right) - \tau_n \boldsymbol{f}{} \left(t_{n+1}, \boldsymbol{y}{}^{n+1}, \mathbf{F}{}_{\text{pl}, h}^{n}, \varepsilon_{\text{pl},h}^{\text{eq}, n} \right) = \mathbf{0}{} \end{align}\tag{23}\] with \(\boldsymbol{\Phi}{}^n\), consisting of solutions on former time steps \(\boldsymbol{y}{}^n, \dots, \boldsymbol{y}{}^{n-k}\), and a constant \(\alpha_{k_n}>0\), which dependents on the chosen time integration order \(k_n\) at time \(t_n\) [27]. With the computed discrete solution \(\boldsymbol{y}{}^{n+1}\), we can then update \(\mathbf{F}{}_{\text{pl},h}^{n+1}\) and \(\varepsilon_{\text{pl},h}^{\text{eq},{n+1}}\) for the new time steps as explained in [9]. Finally, we solve the nonlinear system in each time step with the Newton–Raphson method and follow the adaptive algorithm in [7].

4 Numerical Experiments↩︎

Simulation Setup. We consider a 3D spherical symmetric aSi particle with surrounding spherical symmetric SEI. Therefore, we reduce the computational domain due to symmetry assumptions to the 1D unit interval for the particle \(\Omega_{\text{com}, \text{P}}\) with additional computational SEI domain \(\Omega_{\text{com}, \text{S}}\) as displayed in 1. This means that we have only changes along the radius \(r\), compare also [8] and [4], [7], [9]. For the computational setting, we introduce additional necessary boundary conditions with no flux and no displacement \[\begin{align} \boldsymbol{N}{} \cdot \boldsymbol{n}{}_{\text{P}} = 0, \quad \boldsymbol{u}{}_\text{P}= 0 \quad \text{on }\Gamma_{\text{in}, t_\text{end}}. \end{align}\]

Figure 1: Dimension reduction of a three-dimensional unit sphere with surrounded SEI to the one-dimensional interval, based on [8].

We charge our particle with 1 C and discharge with \(\scalebox{0.75}[1.0]{-}1~\si{C}\) applying the OCV curve stated in [4], [9] with \(c_0 = 0.02\) and a cycling duration of \(0.9~\si{\hour}\). Further initial conditions are selected as follows: \(\mu_0 = \partial_c (\rho_0 \psi_\text{ch}(c_0))\), \(u_\text{P}= r (\lambda_\text{ch}(c_0) - 1)\) and \(u_\text{S}= \lambda_\text{ch}(c_0) - 1\) to decrease the number of required Newton steps at the beginning [7]. All numerical simulations are performed with an isoparametric fourth-order Lagrangian finite element method using an integral evaluation through a Gauß-Legendre quadrature formula with six quadrature points in spatial direction and the finite element library deal.II [29]. Further hardware specifications are given in [9]. Shared memory parallelization is used for assembling of the Newton method. We solve the linear system with a LU-decomposition. Due to implementation reasons of the coupled domain, we use a constant and uniform distributed mesh. In our experiments, the mesh has around 15 × 103 degrees of freedom. The initial time step is chosen as  × 10−8, the maximal time step as  × 10−3, \(\mathop{\mathrm{RelTol}}_t = \num{e-5}\) and \(\mathop{\mathrm{AbsTol}}_t = \num{e-8}\).

Numerical Results. The parameter setup for the particle can be found in Table 2 in [9] with further parameters for the SEI taken from [4] with \(L_{0, \text{S}} = 0.1 L_{0, \text{P}}\), \(\nu_\text{S}= 0.25\), \(E_\text{S}= 900\si{\mega \pascal}\), \(\sigma_\text{Y} = 49.5 \si{\mega \pascal}\). The parameter for viscoplastic plasticity are chosen as \(\dot{\varepsilon}_0 = \num{e-3}~\si{\second}\) or \(\dot{\varepsilon}_0 = \num{e-4}~\si{\second}\), \(\sigma_{\text{Y}^*} = \sigma_\text{Y}\) and \(\beta = 2.94\).

Figure 2: Radial (a) and tangential (b) Cauchy stress over the radius (\(0 \leq r \leq 1.0\) for \(\Omega_{\text{h}, \text{P}}\), \(1.0 \leq r \leq 1.1\) for \(\Omega_{\text{h}, \text{S}}\)) at the end of the simulation; for the 1.) purely elastic (ela.) case with the GSV strain approach (Lag. ela.), the simulation stopped at \(\text{SOC}= 0.34\), see 3 (a), whereas all other cases (2.) elastic, 3.) plastic (pla.) and 4a.), 4b.) viscoplastic (vis.)) for the logarithmic (Log.) Hencky strain approach reached the final simulation time.

Figure 3: Electrical voltage \(U\) over SOC for three half cycles for different elastic strain approaches for elastic and plastic cases with aborted GSV Lagrangian approach for the elastic case in orange (a) and tangential Cauchy SEI stress over SOC at the particle SEI interface with increasing stress-overrelaxation for smaller \(\dot{\varepsilon}_0\) (blue arrow) in the viscoplastic case (b).

We consider in our numerical simulation three half cycles which means two lithiations and one delithiation.

During the model extension from a single particle setting to the coupled particle-SEI approach, we started with a purely elastic SEI with the GSV elastic strain. However, the numerical simulation stopped around \(t = \num{0.32}~\si{\hour}\) and \(\text{SOC} = \num{0.34}\), respectively, failing to find an appropriate Newton update. 2 (a)–(b) and 3 (a)–(b) shows the numerical results in solid orange. It is clearly visible in 2 (b) that the tangential Cauchy stress of the SEI shows a large increase at the interface between the particle and the SEI leading to unnatural behavior and the stop of the numerical simulation. On the other hand and also in the purely elastic setting, the use of the logarithmic strain definition results in a stabilization of the numerical simulation and suppresses the large increase at the interface (solid green). 2 (a) shows the required displacement interface condition between the particle and the SEI. Due to the application of the exponential map of the plastic time integration, the logarithmic strain approach is favorable compared to the GSV strain approach and is used in all plastic simulations. 2 (a)–(b) and 3 (a) show almost no change between the rate-independent and rate-dependent plasticity. However, 3 (b) shows the typical stress-overrelaxation for the tangential Cauchy stress of the SEI at the particle-SEI interface, consisting of an overshooting at the beginning of the first plastification and followed by some relaxation towards the rate-independent plastic results. A decrease of \(\dot{\varepsilon}_0\) from \(\num{e-3}~\si{\second}\) to \(\num{e-4}~\si{\second}\) results in a higher overshooting (blue arrow). An electrical overrelaxation can also be seen in measurements in [6] explaining the voltage relaxation after low current GITT-pulses [5]. The results fit qualitatively to that ones in [4] The purely elastic case has no hysteresis meaning that lithiation and delithiation result in the same stresses, plotted in dotted green in 3 (b). Finally, 3 (a) shows the evaluation of the electric voltage with the Butler–Volmer condition in a postprocessing step leading to similar results compared to the single particle setting [9].

5 Summary and Conclusion↩︎

In this work, we consider a thermodynamically consistent chemo-elasto-plastically coupled model for spherical symmetric aSi particles with surrounded spherical symmetric SEI using a finite deformation approach. We base our theory for the coupling of particle and SEI on [4]. However, we use the rate-independent and rate-dependent plasticity ansatz by [9], which is new for a plastic SEI approach. We apply straightforwardly, for our model extension of the particle-SEI setup with plastic SEI approach, the efficient adaptive temporal algorithm [7] combined with higher order finite element methods on a uniform mesh. During the development of purely elastic effects, we discovered that the Green–St-Venant (Lagrangian) strain approach leads to an abortion during the numerical simulation, whereas a logarithmic strain approach can overcome this failure and results in a successful numerical simulation. The logarithmic approach is also favorable to add plastic deformation due to the used exponential map of the plastic part of the deformation gradient. For the rate-dependent plasticity, the typical stress-overrelaxation can be observed.

In future, we want to consider a spatial adaptive algorithm for the coupled particle-SEI setting with a combined discretization and iteration error estimation [30] or a residual based error estimator [19].

Declaration of competing interest↩︎

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work in this paper.

CRediT authorship contribution statement↩︎

R. Schoof: Methodology, Software, Validation, Formal analysis, Investigation, Data Curation, Writing – original draft, Visualization. G. F. Castelli: Methodology, Writing – review & editing. W. Dörfler: Conceptualization, Resources, Writing – review & editing, Supervision, Project administration, Funding acquisition.


The authors thank L. von Kolzenberg and L. Köbbing for intensive and constructive discussions about modeling silicon particles. R.S. acknowledges financial support by the German Research Foundation (DFG) through the Research Training Group 2218 SiMET – Simulation of Mechano-Electro-Thermal processes in Lithium-ion Batteries, project number 281041241.


R. Schoof: 0000-0001-6848-3844, G. F. Castelli: 0000-0001-5484-6093, W. Dörfler: 0000-0003-1558-9236


S. Sarmah, Lakhanlal, B. K. Kakati, D. Deka, Recent advancement in rechargeable battery technologies, WIREs Energy and Environment 12 (2) (2023) e461.
M. Landstorfer, S. Funken, T. Jacob, An advanced model framework for solid electrolyte intercalation batteries, Phys. Chem. Chem. Phys., 13 (28) (2011) 12817.
Y. Zhao, P. Stein, Y. Bai, M. Al-Siraj, Y. Yang, B.-X. Xu, A review on modeling of electro-chemo-mechanics in lithium-ion batteries, J. Power Sources 413 (2019) 259–283.
L. von Kolzenberg, A. Latz, B. Horstmann, Chemo-mechanical model of sei growth on silicon electrode particles, Batter. Supercaps 5 (2) (2022) e202100216.
L. Köbbing, A. Latz, B. Horstmann, Voltage hysteresis of silicon nanoparticles: Chemo-mechanical particle-sei model, Adv. Funct. Mater.
D. Wycisk, G. K. Mertin, M. Oldenburger, O. von Kessel, A. Latz, Challenges of open-circuit voltage measurements for silicon-containing li-ion cells, Available at SSRN:
G. F. Castelli, L. von Kolzenberg, B. Horstmann, A. Latz, W. Dörfler, Efficient simulation of chemical-mechanical coupling in battery active particles, Energy Technol. 9 (6) (2021) 2000835.
G. F. Castelli, Numerical investigation of Cahn–Hilliard-type phase-field models for battery active particles, Ph.D. thesis, Karlsruhe Institute of Technology (KIT) (2021).
R. Schoof, J. Niermann, A. Dyck, T. Böhlke, W. Dörfler, Efficient modeling and simulation of chemo-elasto-plastically coupled battery active particles, Comput. Mech.
T. Zhang, M. Kamlah, Sodium ion batteries particles: Phase-field modeling with coupling of Cahn–Hilliard equation and finite deformation elasticity, J. Electrochem. Soc. 165 (10) (2018) A1997–A2007.
G. A. Holzapfel, Nonlinear Solid Mechanics, John Wiley & Sons, Ltd., Chichester, 2010.
J. Lubliner, Plasticity theory, Pearson Education, Inc., New York, 2006.
D. Braess, Finite Elements, 3rd Edition, Cambridge University Press, Cambridge, 2007.
R. Schoof, G. F. Castelli, W. Dörfler, Simulation of the deformation for cycling chemo-mechanically coupled battery active particles with mechanical constraints, Comput. Math. Appl. 149 (2023) 135–149.
C. V. Di Leo, E. Rejovitzky, L. Anand, Diffusion-deformation theory for amorphous silicon anodes: The role of plastic deformation on electrochemical performance, Int. J. Solids Struct. 67-68 (2015) 283–296.
A. Bertram, Elasticity and Plasticity of Large Deformations: Including Gradient Materials, 4th Edition, Springer Nature Switzerland AG, Cham,, 2021.
L. Anand, A Cahn-Hilliard-type theory for species diffusion coupled with large elastic-plastic deformations, J. Mech. Phys. Solids 60 (12) (2012) 1983–2002.
C. V. Di Leo,, Ph.D. thesis, Massachusetts Institute of Technology (MIT) (2015). ://
R. Schoof, L. Für, F. Tuschner, W. Dörfler, Residual based error estimator for chemical-mechanically coupled battery active particles, Springer Proc. Math. Stat.
J. Frohne, T. Heister, W. Bangerth, Efficient numerical methods for the large-scale, parallel solution of elastoplastic contact problems, Internat. J. Numer. Methods Engrg. 105 (6) (2016) 416–439.
E. L. Wilson, The static condensation algorithm, International Journal for Numerical Methods in Engineering 8 (1) (1974) 198–203.
D. A. Di Pietro, A. Ern, A hybrid high-order locking-free method for linear elasticity on general meshes, Comput. Methods Appl. Mech. Engrg. 283 (2015) 1–21.
J. C. Simo, T. J. R. Hughes, Computational inelasticity, Interdisciplinary applied mathematics, Springer, New York, 1998.
F.-T. Suttmeier, On plasticity with hardening: an adaptive finite element discretisation, Int. Math. Forum 5 (49-52) (2010) 2591–2601.
P. Neff, C. Wieners, Comparison of models for finite plasticity: a numerical study, Comput. Vis. Sci. 6 (1) (2003) 23–35.
M. W. Reichelt, L. F. Shampine, J. Kierzenka, Matlab ode15s, copyright 1984–2020 The MathWorks, Inc.(1997).
L. F. Shampine, M. W. Reichelt, The MATLABODE suite, SIAM J. Sci. Comput. 18 (1) (1997) 1–22.
L. F. Shampine, M. W. Reichelt, J. A. Kierzenka, Solving index-\(1\)DAEs in MATLAB and Simulink, SIAM Rev. 41 (3) (1999) 538–552.
D. Arndt, W. Bangerth, M. Bergbauer, M. Feder, M. Fehling, J. Heinz, T. Heister, L. Heltai, M. Kronbichler, M. Maier, P. Munch, J.-P. Pelteret, B. Turcksin, D. Wells, S. Zampini, The deal. II library, version 9.5, J. Numer. Math. 31 (3) (2023) 231–246.
T. Carraro, S. Dörsam, S. Frei, D. Schwarz, An adaptive Newton algorithm for optimal control problems with application to optimal electrode design, J. Optim. Theory Appl. 177 (2) (2018) 498–534.