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}\]

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 × 10³ degrees of freedom. The initial time step is chosen as  × 10⁻⁸, the maximal time step as  × 10⁻³, \(\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\).

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.

Acknowledgement↩︎

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.

ORCID↩︎

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

References↩︎