Modeling and Simulation of Chemo-Elasto-Plastically Coupled Battery Active Particles

As an anode material for lithium-ion batteries, amorphous silicon offers a significantly higher energy density than the graphite anodes currently used. Alloying reactions of lithium and silicon, however, induce large deformation and lead to volume changes up to 300%. We formulate a thermodynamically consistent continuum model for the chemo-elasto-plastic diffusion-deformation behavior of amorphous silicon and it's alloy with lithium based on finite deformations. In this paper, two plasticity theories, i.e. a rate-independent theory with linear isotropic hardening and a rate-dependent one, are formulated to allow the evolution of plastic deformations and reduce occurring stresses. Using modern numerical techniques, such as higher order finite element methods as well as efficient space and time adaptive solution algorithms, the diffusion-deformation behavior resulting from both theories is compared. In order to further increase the computational efficiency, an automatic differentiation scheme is used, allowing for a significant speed up in assembling time as compared to an algorithmic linearization for the global finite element Newton scheme. Both plastic approaches lead to a more heterogeneous concentration distribution and to a change to tensile tangential Cauchy stresses at the particle surface at the end of one charging cycle. Different parameter studies show how an amplification of the plastic deformation is affected. Interestingly, an elliptical particle shows only plastic deformation at the smaller half axis. With the demonstrated efficiency of the applied methods, results after five charging cycles are also discussed and can provide indications for the performance of lithium-ion batteries in long term use.


Introduction
Lithium (Li)-ion batteries gained an enormous amount of research interest in the past two decades [1], as a mean of storing electric energy and propelling electromobility [2,3].However, due to the complex electro-chemo-mechanically coupled processes occurring during charging and discharging of Li-ion batteries, ongoing research still aims at improving battery lifetime, reducing costs and increasing capacity by, e.g., varying the materials composing the battery [1,3].State of the art is the usage of graphite as anode material [1].A promising candidate to be used as anode material in Li-ion batteries is amorphous silicon (aSi), due to its large capacity and capability to form an amorphous Li x Si alloy [4] with the diffusing Li-ions, increasing battery capacity [2].A disadvantage is the large volume increase aSi particles undergo during alloying, which can reach up to 300 % [5].Numerous simulative studies have shown that these large deformations are accompanied by plastic deformations of aSi which are inherently linked to battery lifetime and capacity, see e.g.[3,[6][7][8][9].With the goal of using aSi as anode material, it is therefore imperative to study plastic deformation mechanisms at the particle level of aSi anodes and their interplay with battery performance during charging and discharging using physical models and computational investigations.
To this end, geometrically and physically nonlinear chemo-mechanically coupled continuum theories have proven to be a valuable tool, see e.g.[7,[10][11][12].For the mechanical part of the model, most works rely on a multiplicative split of the deformation gradient into a chemical, an elastic and a plastic part using finite deformation, see e.g.[13,Section 10.4] and [14, Section 8.2.2], [7,15].Discrepancies in modeling strategies occur in the nonlinear strain measure used, ranging from the Green-Lagrange strain tensor [11,12] to the Hencky strain tensor [7].In addition, several models consider plastic deformation to be rate-dependent [7,8], while others rely on a rate-independent plasticity theory [10,16].Unfortunately, neither the atomic-level structural evolution, nor the mechanical behavior of the aSi during lithiation and delithiation cycles is well understood [17].This also holds for the detailed mechanism of plastic deformation.However, several studies concluded that plasticity does occur during charging and discharging.In experimental studies, c.f. [18], a rate-dependent plastic behavior is considered to explain the observed behavior.In contrast, a numerical study conducted on a molecular level in [19] seems to indicate rate-independent plasticity.The chemical part of the models, describing diffusion of Li-ions during charging and discharging, is based on a diffusion equation relating changes in concentration to the gradient of the species' chemical potential and the species' mobility [7].Models differ in their approach to define the chemical contribution to the Helmholtz free energy, where approaches either rely on open-circuit voltage (OCV) curves [12] or assumptions for the entropy of mixing [7].In addition, the mobility is defined to be either derived as the change of the chemical part of the chemical potential with respect to the concentration [7] or the entire chemical potential [15].The coupling of deformations and diffusion arises due to the strains induced by Li-ions as well as the influence of mechanical stresses on the chemical potential.
Both finite difference [6,15] and finite element [7,11] schemes have been proposed to discretize the resulting equations, where the latter have been predominantly used lately, due to their superior applicability to complex geometries.Solving discretized non-linear coupled systems of equations is time consuming and expensive in terms of computational resources, due to small mesh sizes and small time step sizes required to resolve all mechanisms.In contrast to a finite difference discretization [15] or constant spatial discretization with uniform temporal discretization [20] space and time adaptive solution algorithms, such as the one proposed in [11], allow drastic reduction in computational resources.In addition, parallelization schemes [21] and automatic differentiation (AD) can reduce simulation times considerably.The AD concept is presented for the first time in combination with the variable-step, variable-order time integration scheme and the inelastic constitutive theory for a battery active particle.Introducing plastic deformation is another challenge, as the additional variables are either considered as degrees of freedom [15] or static condensation is used [22][23][24][25] to arrive at a primal formulation [7,26], where the variables are only computed at integration point level.A final remark concerns the choice of the chemical degree of freedom in chemo-mechanically coupled models.In non-phase separating models, commonly, only a single chemical degree of freedom is introduced, which is either the concentration or the chemical potential of the diffusing Li [7,27].However, when phase separation occurs, often modeled by a Cahn-Hilliard type theory, splitting methods [11] or micromorphic approaches [28] introduce a second chemical degree of freedom in order to circumvent the discretization of a partial differential equation of order four and stabilize the numerical scheme.
The goal of this work is to investigate predictions of a chemo-mechanically coupled model for large chemo-elasto-plastic deformation processes in aSi anode particles using modern numerical techniques including space and time adaptivity, parallelization schemes and AD.Regarding the model, we take into account plastic deformation of the aSi particles, where we consider the initial yield stress to be a function of lithium concentration [19].As no consensus exists in the experimental literature regarding the mechanisms of plastic deformation in aSi, we formulate both a rate-dependent viscoplasticity, as well as a rate-independent plasticity theory and discuss the different implications on particle behavior [17][18][19].We use static condensation to arrive at a primal formulation for the mechanical equations and consider plastic deformation at integration point level within each finite element [26].We explicitly derive a projector onto the admissible stresses, relying on the classical return mapping method [29,Chapter 3] and [26,30,31], which is rather straightforward for the Hencky strains used in our theory [32,33].The diffusion of Li into and out of the aSi anode particle follows classical diffusion theory.However, instead of relying on a standard mixing entropy ansatz for the chemical part of the free energy density, e.g., [7], we rely on a measured experimental OCV curve to model the chemical part [12,15,21], which has not received a lot of attention in literature.To the best of the authors' knowledge, a combination of a measured experimental OCV curve with a (visco-)plastic battery active particle has not been used before.In this work, we limit our model to show numerical results with a single-phase mechanism after the initial two-phase transformation of aSi in the first half cycle as described in [34].This finding of [34] is in contrast to [4], where a phase-separation is also measured after the first half cycle.As a boundary condition, various charging rates (C-rates) are applied for the lithium flux.For solving the coupled system of equations we extend the solution scheme proposed in [11,21,35], relying on a spatial and temporal adaptive algorithm.We consider a one-dimensional computational domain with radial symmetry and a two-dimensional elliptical computational domain and compare stress and plastic strain development as well as concentration distributions after a various number of half cycles for both rate-dependent and rate-independent plasticity models.In addition, we investigate the computational performance and numerical efficiency of our implementation scheme.
The remainder of this article is organized as follows: in Section 2 we introduce the theoretical basis for our work and derive the equations describing chemo-mechanically coupled diffusion processes in aSi anodes.Section 3 summarizes the numerical approach taken in this work to solve the derived system of equations.Subsequently, in Section 4, we present results for various investigated cases.We close with a conclusion and an outlook in Section 5.

Theory
In a first step we review and summarize our used constitutive theory adapted from [7,11,12,15] to couple chemical, elastic and plastic material behavior.We base our model on a thermodynamically consistent theory for the chemo-mechanical coupling during lithiation and delithiation.

Finite Deformation
We start by considering a mapping Φ : R ≥0 ×Ω 0 → Ω, Φ (t, X 0 ) := x = X 0 + u(t, X 0 ) from the Lagrangian domain Ω 0 to the Eulerian domain Ω, see for more information [36, Section 2], [14, Section 8.1], [37, Chapter VI] and [7,11,12,15,28].The deformation gradient F = ∂Φ/∂X 0 = I + ∇ 0 u with the identity tensor I and displacement u is multiplicatively decomposed into chemical, elastic and plastic parts [13,Section 10.4] and [14, Section 8.2.2], [7,15] with various expressions for the volume changes defined as respectively.The chemical and elastic deformations are reversible and summarized in F rev .The polar decomposition of the elastic part F el of the deformation gradient is given by its rotational and stretch part [36,Section 2.6]: with the right stretch tensor U el being unique, positive definite and symmetric.With the symmetric elastic right Cauchy-Green tensor the (Lagrangian) logarithmic Hencky strain can be defined as strain measure with a spectral decomposition where √ η el,α and r el,α are the eigenvalues and eigenvectors of U el , respectively.In literature, typically the Green-St-Venant (GSV) strain tensor, often called the Lagrangian strain tensor [14, Section 8.1], is used We will later compare results obtained for both strain measures.We consider an isotropic, volumetric swelling due to the Li concentration with F ch being defined as where λ ch = 3 1 + v pmv c max c, v pmv is the constant partial molar volume (pmv) of lithium inside the host material [11] and c ∈ [0, 1] is the normalized concentration in the reference configuration regarding the maximal concentration of the host material c max of the reference domain in mol m −3 .Thus, the total concentration in reference configuration is c = c max c in mol m −3 .With the chemical part F ch of the total deformation gradient and the plastic part F pl (t, X 0 ), introduced as an internal variable [38] in more detail in Subsection 2.5, the elastic part F el of the deformation gradient is given by

Free Energy
To obtain a thermodynamically consistent material model, which guarantees a strictly positive entropy production, we introduce a Helmholtz free energy ψ, being in general a function of the lithium concentration c and the displacement gradient ∇ 0 u [12,15,[39][40][41] as well as the internal variable F pl [38], while we assume the temperature T to be constant.Therefore, we use an additive decomposition into a chemical part ψ ch , being unaltered by elastic deformations, and a mechanical part ψ el , depending solely on elastic deformations, according to [15,[39][40][41] in the Lagrangian frame.
Following [12,15,21], we define the chemical part by incorporating an experimentally obtained open-circuit voltage (OCV) curve with the Faraday constant Fa and the mass density ρ 0 of aSi in the reference configuration.This measurement based model for the contribution to the chemical part of the free energy density is in contrast with most existing works, as often the chemical contribution is based on mixture entropy approaches for ideal gases, see, e.g., [6,7].The assumption of Li behaving like an ideal gas could be considered questionable, motivating the approach to base the chemical contribution part on a measurement.This leads to a chemical potential being non-symmetric to c = 0.5, c.f. Figure 11, which is not the case for the mixture entropy model.The mechanical part is given by a linear elastic approach being quadratic in the elastic Hencky strain [20,28,32]: Here,

Chemistry
Inside the host material, we use a continuity equation to describe the change in lithium concentration in the reference configuration via where is the scalar mobility, where m > 0, of Li in aSi for the considered isotropic case and D the diffusion coefficient for lithium atoms inside the active material [12,15].The chemical potential µ(t, X 0 ) is given as partial derivative of the free energy density with respect to the total concentration c using Equation ( 9)-( 11) Following [7,11], we apply a uniform and constant external flux N ext as constant current with either positive or negative sign for cycling the host particle in terms of the charging rate (C-rate).The C-rate characterizes the hours of complete lithiation of the particle: C-rate = 1/t cycle .So a charging rate of 1 C means that the battery is fully lithiated in 1 h, of 0.5 C in 2 h and so on [42].N ext is given in the dimension 1/h after scaling with c max /A V for the specific surface A V = surface/volume of Ω 0 in m 2 m −3 [11].The simulation time t and the state of charge (SOC) can be connected via with the volume V Ω0 of Ω 0 and a constant initial condition c 0 ∈ (0, c max ).

Mechanics
The deformation in the Lagrangian domain is considered by static balance of linear momentum [11,12,15] with the first Piola-Kirchhoff stress tensor The Cauchy stress σ in the Eulerian frame is linked to P via σ = PF T / det (F) [36,Section 3.1].Furthermore, we introduce the Mandel stress M(c, ∇ 0 u, F pl ) = C rev S rev = J rev F T rev σF −T rev = J el J ch F T el σF −T el with the second Piola-Kirchhoff stress tensor S rev = J rev F −1 rev σF −T rev , see for further information [15].Based on the derivations presented in [7,20,28], for C el and Ċel being coaxial and isotropic material behavior, a hyperelastic law relating the free energy density in the stress-free configuration and the Mandel stress M is retrieved.In the case considered in this work M is linear in E el and given by

Inelastic Constitutive Theory
Following [36, Section 2.7] and [7], the evolution equation for the plastic part F pl of the deformation gradient takes for elastically and plastically isotropic materials (where a plastic spin is negligible) the form As mentioned in the introduction, we consider two inelastic models and compare their influence on battery performance.We start with a rate-independent von Mises plasticity with isotropic hardening, which is formulated for the Mandel stress, see [29,Section 2.3] and [7,15,26,31,43].The yield function reads with the deviatoric stress tensor M dev = M − 1/3 tr M I, the Frobenius norm || • || and the yield stress σ F (c, ε eq pl ) := σ Y (c) + γ iso ε eq pl .The yield stress consists of two parts: a concentration dependent part σ Y (c), which will describe a softening behavior, and a linear isotropic hardening part with a scalar parameter γ iso > 0. ε eq pl (t, X 0 ) is the accumulated equivalent inelastic strain and will be introduced subsequently.Ideal plasticity is present for γ iso = 0.The softening behavior modeled by σ Y (c) is inspired by [7,19] and is incorporated via Plastic flow is only allowed if F Y = 0. To describe the plastic flow when this yield point is reached, we base on the maximum plastic dissipation principle [15,29,44].With this postulate from plasticity theory we can define the associated flow rule constraining the plastic flow to the normal direction of the yield surface ∂ M F Y : Here, D pl is the plastic strain rate measure [7,45].Now, we can define the scalar equivalent plastic strain which is used to describe an increase in yield stress σ F (c, ε eq pl ).Note that we always start with ε eq pl (0, •) = 0. To be consistent with a one-dimensional tensile test, we scale our concentration dependent yield stress σ Y (c) with the factor 2/3, compare [29,Section 2.3.1].
The second model studied is a viscoplastic material model which was proposed by Di Leo et al. [7], where a viscoplastic material behavior is considered without isotropic hardening, i.e. γ iso = 0.This results in a formulation for the equivalent plastic strain where σ Y * , ε0 and β are a positive-valued stress-dimensioned constant, a reference tensile plastic strain rate and a measure of the strain rate sensitivity of the material, respectively.
The classical loading and unloading conditions for rate-independent plasticity can be conveniently expressed via the Karush-Kuhn-Tucker (KKT) conditions [29, Section 1.2.1], [14, Section 3.2] and [15] Compared to classical notation of loading and unloading, the process is elastic if F Y < 0 requiring εeq pl ≡ 0 and no plastic deformation occurs.The consistency condition for the evolution of inelastic strains in the case of rate-independent plasticity reads when F Y = 0 : εeq pl ≥ 0, ḞY ≤ 0, εeq pl ḞY = 0, (26) so the plastic strain can increase during loading but not during unloading.In the viscoplastic case, the KKT conditions and the consistency condition are replaced by the evolution equation of the equivalent plastic strains as given by Equation ( 24), compare [29, Section 1.7].

Numerical Approach
In the following section we present the numerical treatment of our set of coupled partial differential equations of Section 2, i.e. the problem formulation, the normalization of the model parameters and the numerical solution procedure including the weak formulation, space and time discretization and our applied adaptive solution algorithm.

Problem Formulation
Before we state our problem formulation we introduce a nondimensionalization of the model to improve numerical stability.Our cycle time t cycle = 1/C-rate depends on the C-rate, i.e. the hours for charging or discharging of the particle.Further, the particle radius L 0 and the maximal concentration c max in the Lagrangian frame are used as reference parameters.For the yield stress σ Y (c) we use the same nondimensionalization as for the Young's modulus E. The resulting dimensionless numbers Ẽ and the Fourier number Fo relate the mechanical energy scale to the chemical energy scale and the diffusion time scale to the process time scale, respectively.All dimensionless variables are listed in Table 1 and will be used for model equations from now on, neglecting the accentuation for better readability.We state our general mathematical problem formulation [11,12] by solving our set of equations for the concentration c, the chemical potential µ and the displacements u, whereas the quantities F, F el , E el and P are calculated in dependency of the solution variables.The Cauchy stresses σ are computed in a postprocessing step.The handling of the internal variable F pl [38] is explained in more detail in Subsection 3.2.
The dimensionless initial boundary value problem with inequality boundary conditions is given as follows: let t end > 0 be the final simulation time and Ω 0 ⊂ R d a representative bounded electrode particle in reference configuration with dimension d = 1, 2, 3. Find the normalized concentration c : [0, with a boundary-consistent initial concentration c 0 and boundary conditions for the displacement excluding rigid body motions.Note that the original definition of the chemical part F ch of the deformation gradient is done in three dimensions, but all variables and equations are also mathematically valid in dimensions d = 1, 2. In this case, the deviatoric part is computed with the factor 1/d in contrast to the factor 1/3 after Equation (20).In case of d = 2, for example, we have no displacement and no stresses in the third direction, so it is purely a two-dimensional computational consideration.A final remark: Equation (27a) and Equation (27b) could be condensed to one equation.The second equation would be treated on integration point level.However, we choose the splitting method, detailed in [11], to easily extend the physical behavior to phase-separating materials with a Cahn-Hilliard approach and stabilize the numerical solution procedure.Using only the chemical potential as solution variable would lead to an additional calculation of the concentration with a scalar Newton method.
Table 1: Dimensionless variables of the used model equations.

Numerical Solution Procedure
This subsection describes the way to obtain a numerical solution and especially the handling of the KKT condition in Equation (27d): formulating a primal mixed variational inequality, using static condensation to obtain a primal formulation as well as space and time discretization, finally completed with an adaptive solution algorithm.

Weak Formulation
In a first step towards the numerical solution we state the weak formulation of Equation ( 27) as primal mixed variational inequality like in [26].However, it can also be derived from a minimization problem [29, Section 1.4.2] and [46,Section 7.3].We introduce the L 2 -inner product for two functions f , g ∈ L 2 (Ω 0 ) as (f, g) = Ω0 f g dX 0 , for two vector fields v, w ∈ L 2 Ω 0 ; R d as (v, w) = Ω0 v • w dX 0 , and for two tensor fields S, T ∈ L 2 Ω 0 ; R d,d as (S, T) = Ω0 S : T dX 0 and boundary integrals with the respective boundary as subscript.Now, we define the function spaces being sufficiently large in order to allow a mathematically well-posed formulation [38].V * includes displacement boundary constraints for the precise applications case from Section 4. Next, we multiply with test functions, integrate over Ω 0 and integrate by parts.Following [11,26,44,45], we finally have the weak primal mixed variational inequality formulation: find solutions {c, µ, u} and internal variables F pl and ε eq pl with c, µ ∈ V , for all test functions φ ∈ V , ξ ∈ V * and P * , ε eq * pl ∈ Y. Equation ( 28) becomes a saddle point problem requiring special techniques for solving the related linear system [26,37,47].However, we prefer a version of Equation (28) without Equation (28d) for easier handling of the numerical studies.Therefore, we eliminate Equation (28d) by using a primal formulation with a projector onto the set of admissible stresses [26,48], also known as static condensation [22][23][24][25].Following [26], we introduce the continuous projector onto the admissible Mandel stress for both considered inelastic constitutive theories.For the rate-independent model, it reads with κ(c, ε eq pl ) = 1 − 2G σY(c) ε eq pl and M tri follows from a purely elastic deformation, denoted as the trial part of M. The projector for ideal plasticity follows with γ iso = 0, while the projector for the rate-dependent viscoplastic approach is given in Subsection 3.2.3.
Finally, we can reformulate Equation ( 28) using the projector formulation for the specific plastic behavior and arrive at the primal formulation: for given F pl and ε eq pl This means that inequality Equation ( 27d) is condensed into Equation (30c) with Equation (29) or Equation (39), respectively.Note that by using the projector P Π we have transformed the plasticity inequality into a (non-smooth, actually Lipschitz continuous) nonlinear equation.

Space Discretization
Next, we introduce the spatial discretization and therefore choose a polytop approximation Ω h as computational domain for particle geometry Ω 0 .For the approximation of curved boundaries, an isoparametric Lagrangian finite element method is chosen [37, Chapter III §2] on an admissible mesh T n .For the spatial discrete solution we define the finite dimensional subspaces for the basis functions with bases On these finite dimensional subspaces we solve for the variables (30).Relating the discrete solution variables with the finite basis functions we gather all time-dependent coefficients in the vector valued function The internal variables F pl and ε eq pl are discretized in the same way as the solution variables and are denoted by F pl,h and ε eq pl,h .Note that both internal variables F pl and ε eq pl are not part of the solution vector y.However, they are handled separately as introduced in Subsection 3.2.3.Both projector formulas of Equation ( 29) and Equation ( 39) are also valid in the discrete case.This leads to our spatial discrete problem formulated as general nonlinear differential algebraic equation (DAE): for given F pl,h , ε eq pl,h The system matrix Mh is singular since it has only one nonzero-block entry given by M h = [(φ i , φ j )] i,j representing the mass matrix of the finite element space V h .The vector f consists of matrices and tensors, given as with the same indices as above: the mass matrix M h , the stiffness matrix

Time Discretization
Before we formulate the space and time discrete problem for y, we have to consider the time evolution of the internal variables F pl and ε eq pl and derive a time integration scheme for them.Therefore, we update the time integration separately from the time advancing of our Equation (34).Applying an implicit exponential map to Equation (19) to the continuous internal variable F pl leads to [38] from one time step t n to the next t n+1 = t n + τ n with time step size τ n > 0. For the rate-independent plasticity we use the well known return mapping algorithm [32,33] in an explicit form: where the trial Mandel stress is ] and the difference of the new time step and current time step of ε eq pl is △ n ε eq pl = ε eq,n+1 pl − ε eq,n pl .This is straightforward for isotropic linear hardening.With the solution for ε eq,n+1 pl from Equation (38) and the initial conditions F pl (0, •) = I and ε eq pl (0, •) = 0, all necessary quantities can be updated and F n+1 pl can be computed.For more details regarding the time integration scheme for F pl in the rate-independent case, see [20,Appendix C.5].
Following [7,28] for the rate-dependent viscoplastic approach, the projector is For this case, however, no explicit form can be retrieved for the accumulated plastic strain, due to the nonlinearity in Equation (24b).Therefore, we have to use a scalar Newton-Raphson method for the current time increment τ n for △ n ε eq pl of the implicit Eulerian scheme: solve Equation (24b) for M tri,dev > σ Y (c n+1 ) using the relation ||M n+1,dev || = ||M tri,dev || − 2G△ n ε eq pl from Equation (37) with the residual of Equation (24b) All computations from Equation (36) to Equation ( 40) are also valid for the spatial discrete variables but are written for the continuous case for sake of readability.
For the time evolution of the DAE (34), we apply the family of numerical differentiation formulas (NDFs) in a variable-step, variable-order algorithm, based on the Matlab's ode15s [49][50][51][52], since our DAE has similar properties as stiff ordinary differential equations [11].An error control handles the switch in time step sizes τ n and order.We have chosen the same time scale for this temporal discretization and the temporal discretization of the internal variables F pl and ε eq pl .We arrive at the space and time discrete problem to go on from one time step t n to the next t n+1 : for given F n pl,h and ε eq,n pl,h find the discrete solution y n+1 ≈ y(t n+1 ) satisfying and ε eq,n pl,h to the new time steps F n+1 pl,h and ε eq,n+1 pl,h as already introduced.

Adaptive Solution Algorithm
Since our DAE ( 34) is nonlinear, we apply the Newton-Raphson method and thus need to compute the Newton update of the Jacobian in each time step.For this, we have to linearize our equations, especially the projectors of Equation ( 29) and Equation (39).However, both projectors are not formally differentiable.Nevertheless, they fulfill the conditions to be slantly differentiable [26,53] and we can work with the formal linearization as an approximation.For the projector defined in Equation ( 29), we follow [20,Appendix C.6] and propose the consistent algorithmic modulus of P Π (c, ∇ 0 u, F pl , ε eq pl ) around E tri el [26,30] as with C G = 2G I − 1 3 I ⊗ I , C K = KI ⊗ I, the bulk modulus K and κ as in Subsection 3.2.1.The derivation for the approximation of the linearization of Equation ( 42) is given in Appendix C as well as the formulation of the linearization for the projector of Equation (39).Again, both linearizations are written for continuous functions but are also valid for the discrete variables.
Another possibility is the automatic differentiation (AD) framework, provided by [54], which offers the possibility to use an AD framework via Sacado (a component of Trilinos) [55].This allows us to compute the partial derivative, required for the Newton scheme, automatically instead assembling the Newton matrix by hand.We use the tapeless dynamic forward-mode Sacado number type (one time differentiable).To compute the Newton update, we use a direct LU-decomposition.The iteration number can be decreased if an appropriate initialization is chosen.Therefore, the starting condition for the first time step is stated in Subsection 4.1 while a predictor scheme is used during the further time integration [50].
Finally, we follow Algorithm 1 in [11] for the space and time adaptive solution algorithm.Both a temporal error estimator [49][50][51][52] and a spatial error estimator are considered for the respective adaptivity.To measure spatial regularity, a gradient recovery estimator is applied for the solution variables c, µ and u [56, Chapter 4].For marking the cells of the discrete triangulation for coarsening and refinement, two parameters θ c and θ r are applied with a maximum strategy [57].Altogether, we use a mixed error control with the parameters RelTol t , AbsTol t , RelTol x and AbsTol x .For further details we refer to [11].

Numerical Studies
This section deals with the investigation of the presented model from Section 2 with the numerical tools of Section 3.For this purpose we introduce the simulation setup in Subsection 4.1 and discuss the numerical results in Subsection 4.2 with 1D and 2D simulations, which are a 3D spherical symmetric particle reduced to the 1D unit interval and in addition a 2D quarter ellipse reduced from a 3D elliptical nanowire, respectively.The latter one is chosen to reveal the influence of asymmetric half-axis length on plastic deformation.

Simulation Setup
As mentioned in the introduction, amorphous silicon is worth investigating due to its larger energy density and is therefore chosen as host material.The used model parameters are listed in Table 2. Most parameters are taken from [15,21], however, for the plastic deformation we pick and adapt the parameters from [7] such that the yield stress is in the range of [58].In particular, the ratio between σ Y,max and σ Y,min is maintained.If not otherwise stated, we charge with 1 C and discharge with −1 C. Following [15], we charge between U max = 0.5 V and U min = 0.05 V, corresponding to an initial concentration of c 0 = 0.02 and a duration of 0.9 h for one half cycle which is one lithiation of the host particle.The OCV curve U OCV (c) for silicon is taken from [59] and defined as U OCV : (0, 1) → R >0 with The curve is depicted in Appendix G in Figure 11.

Geometrical Setup
We proceed by presenting two computational domains and present the boundary conditions due to new artificial boundaries.We choose a representative 3D spherical particle and reduce the computational domain to the 1D unit interval Ω 0 = (0, 1), with a new artificial boundary Γ 0 , compare Figure 1(a), in the particle center with a no flux condition and zero displacement: To ensure the radial symmetry, we adapt the quadrature weight to dX 0 = 4πr 2 dr in the discrete finite element formulation.In the 1D domain it is consistent to assume that the fields vary solely along the radius r.For the calculation of the displacement gradient, we refer to [  In Figure 1(b), the 2D simulation case is shown in terms of a quarter ellipse.We choose a factor of 0.6 for the short half-axis compared to the longer half-axis.We create this geometry by considering a 3D nanowire with no stresses and no changes in z-direction as well as symmetry around the x-and y-axes.Here, further artificial boundaries on Γ 0,x and on Γ 0,y with no flux conditions and only radial displacement: have to be introduced.We make use of an isoparametric mapping for the representation of the curved boundary on Γ ext .Again, we choose a constant initial concentration c 0 = 0.02 and a chemical potential µ 0 = ∂ c ρ 0 ψ ch (c 0 ) .However, we use for the displacement the condition u 0 = 0.

Implementation Details
All numerical simulations are executed with an isoparametric fourth-order Lagrangian finite element method and all integrals are evaluated through a Gauß-Legendre quadrature formula with six quadrature points in space direction.Our code implementation is based on the finite element library deal.II [54] , an initial refinement of seven and of four for the 1D case and the 2D case, respectively, an initial time step size τ 0 = 1 × 10 −6 and a maximal time step size For the marking parameters of local mesh coarsening and refinement, θ c = 0.05 and θ r = 0.5 are set.A minimal refinement level of five is applied for the 1D simulations, in order to achieve a parameterization as general as possible for all different 1D simulation.Due to limited regularity of the nonlinearity of the plastic deformation, we limit the maximal order of the adaptive temporal adaptivity by two.

Numerical Results
In this section we consider the numerical results of the 1D spherical symmetric particle and the 2D quarter ellipse computational domain.We analyze the computed fields such as stresses and concentrations as well as the computational performance of our presented model and implementation scheme.Further, we compare the computational times for using the derived linearization of the projector formulation and an automatic differentiation (AD) technique.

1D Spherical Symmetry
In a first step, we analyze the effect of plastic deformation on the chemo-physical behavior of the 1D domain depicted in Figure 1(a).Detailed studies for purely elastic behavior can be found in, e.g., [12,21] and are included in this study for comparison.Physical results after one half cycle.In Figure 2, we compare the numerical results for the concentration, the plastic strains as well as the tangential Cauchy stresses between the elastic, plastic and viscoplastic model after one lithiation of the host particle, that means one half cycle at SOC = 0.92 over the particle radius r.The changes of the concentration profiles due to plastic deformation are displayed in Figure 2(a).It is clearly visible that the concentration gradient increases in both the plastic and viscoplastic case in the vicinity of the particle surface, whereas lower concentration values occur in the particle center compared to the elastic case.This is due to the limited maximal stress in the plastic models.The lower stresses inside the particle lead to a lower mobility in the particle interior, c.f. Equation (13).This gradient in the Li mobility leads to the observed pile up of Li atoms at the particle surface.Comparing the plastic and viscoplastic case, a smoother transition from elastic to plastic is visible and therefore a little shift of the concentration values to the particle surface, see the lower magnifying glass in Figure 2(a).This is commonly referred to as viscoplastic regularization [29, Section 1.7.1.4.].The second magnifier shows an area, where the slope of the concentration profile changes close to the particle surface, revealing a second plastic deformation process during lithiation, which is also observable at the change in slope of the equivalent plastic strain, c.f. Figure 2(b).This becomes more apparent, when the tangential Cauchy stresses are investigated as a function of the SOC, c.f. Figure 2(d) and Figure 3(b).
Figure 2(b) shows the equivalent plastic strain ε eq pl , revealing plastic deformations near the particle surface.During lithiation, the particle deforms plastically twice.The first plastic deformation process occurs in the initial stages of charging at low SOC ≤ 0.13, c.f. Figure 2(d), and leads to equivalent plastic strains of 3.4 %.Upon further lithiation, this process repeats at SOC ≈ 0.8, c.f. Figure 2(d SOC / Tangential stress σ φ / GPa Fig. 2: Numerical results for the elastic (Ela.), plastic (Pla.) and viscoplastic (Vis.)approaches of the 1D radial symmetric case at SOC = 0.92 over the particle radius r: concentration c in (a), equivalent plastic strain ε eq pl in (b) and tangential Cauchy stress σ ϕ in (c) as well as tangential Cauchy stress σ ϕ at the particle surface r = 1.0 over SOC in (d).stress σ ϕ , c.f. Figure 2(c) and (d), there is a change of the stress direction from compressive stress to tensile stress at the particle surface for the plastic cases compared to the elastic case.This change in sign for the tangential stresses occurs in the area, where the particle undergoes plastic deformation beforehand.The heterogeneous plastic deformation thus leads to an eigenstrain, that results in tensile stresses near the particle surface, which cannot be observed in the elastic case.This means that for an almost fully lithiated particle, there is a significant shift in the stress development and the plastic deformation leads to tensile stresses close to the particle surface.This crucial change in the stress profile at a small spatial area of the particle is important to recognize for the battery life time.Figure 2(d) displays the tangential Cauchy stress at the particle surface versus the SOC over the complete half cycle.Directly after the start compressive Cauchy stresses occur where the elastic approach increases to larger values compared to the plastic approaches, which are limited due to the onset of plastic deformations.The viscoplastic model shows larger negative tangential stresses, i.e. an overstress above the yield stress, when compared to the rate-independent models, also allowing larger elastic strains, c.f. Figure 3(a).After reaching the maximal values the stresses reduce in all cases.However, the plastic approaches predict tensile stresses around SOC = 0.6.At SOC ≈ 0.8, the particle deforms plastically a second time, reducing the maximal tensile stresses in tangential direction.Striking is the difference between the elastic case revealing only compressive tangential Cauchy stresses compared to the plastic approaches featuring also tensile tangential Cauchy stresses.The radial Cauchy stress is not plotted since we have stress-free boundary condition at the particle surface.These findings are qualitatively comparable to numerical results from [61], Fig. 4(c), and [6], Fig. 5(d), compare also the numerical results at SOC = 0.1 in Figure 8 of Appendix D. We want to point out that both the plastic and the viscoplastic model lead to almost equivalent numerical results at the end of the first half cycle.This is, however, not the case for results after multiple half cycles, which is outlined below.Before we proceed by investigating the influence of several material parameters and multiple half cycles, we also compare the Green-St-Venant strain tensor or Lagrangian strain with our used logarithmic strain tensor.Both approaches predict almost identical results, see Appendix E, as also observed in the numerical results in [62].
Parameter studies.In a next step we take a closer look at the stress-strain curves of the different mechanical approaches and analyze the influence of the maximal yield stress σ Y,max on the tangential Cauchy stress σ ϕ at the particle surface r = 1.0 in Figure 3. Furthermore, we compare the dependency on the C-rate and the particle size in Figure 4.The AD concept is used for all parameter studies, see the next subsection numerical efficiency for more details.The comparison in Figure 3(a) shows the effects of the different mechanical approaches: elastic, plastic and viscoplastic deformation as well as ideal plastic deformation with γ iso = 0.For the ideal plastic case, we choose a uniform grid with ten refinements and a backward differentiation formula (BDF) time stepping of order two to increase numerical performance.Comparing the elastic and all plastic cases, the maximal compressive Cauchy stress of the elastic solution is not reached in the plastic cases, however, the stresses decrease rapidly after reaching the yield stress.Again, the viscoplastic overstress becomes apparent, c.f. Figure 2(d).Further, the influence of the concentration dependent yield stress is clearly visible, since with a constant yield stress, the curve would just move straight upwards (red dotted reference line) instead featuring a shift to the right.The tangential stress decreases rapidly after an initial plastic deformation, see also Figure 3(b), so that no further plastic deformation occurs.In contrast to the elastic model, tangential tensile stresses occur for large SOC values in all plasticity models, which lead to another onset of plastic deformation.For larger SOC values, only the plasticity approaches show tensile tangential Cauchy stresses with a second plastic deformation at the end of the first lithiation.As stated below, the viscoplastic model including the viscoplastic regularization leads to better numerical properties.In addition, the experiments of Pharr et al. [18] indeed indicate some sort of rate-dependent inelastic behavior, which, in the opinion of the authors, makes the usage of a viscoplastic model more plausible.Therefore, we continue our investigation with the viscoplastic model unless otherwise stated.Figure 3(b) compares the influence of a varying maximal yield stress σ Y,max on the tangential stress at the particle surface over the SOC.For smaller values of σ Y,max , the particle starts yielding at smaller tangential stresses, leading to an overall decrease in the observed minimal tangential stresses σ ϕ at small SOC.In contrast, the earlier the plastic deformation occurs, the larger the tangential tensile stresses at higher SOC.In addition, the decrease in yield stress with increasing concentrations can be observed, as the tangential tensile stresses lead to further plastic deformation at lower stress levels, indicated with the cyan arrow in Figure 3(b).However, no plastic deformation is visible for a maximal yield stress of 1.0 GPa, which shows a purely elastic response.
Next, we analyze the dependency of the plastic deformation on different C-rates and particle sizes in Figure 4 a comfortable user experience.However, Figure 4(a) shows that for higher C-rates higher tangential Cauchy stresses arise which lead especially for higher SOC values to a large area of plastic deformation.The smaller the C-rate, the lower the occurring stresses and the smaller the plastically deformed regions in the anode particle.For the parameter set considered in this study, decreasing the C-rate by 50 % leads to purely elastic deformations.Similar results are observed, when an increasing particle diameter is considered, as in Figure 4(b): the larger the particles, the greater the stresses and the area of plastic deformation.This results from larger concentration gradients for larger particles since the lithium needs longer to diffuse to the particle center.The simulation for particle size L 0 = 200 nm terminates at SOC ≈ 0.55, since at this simulation time a concentration of one is reached at the particle surface.For the largest particle radius, the heterogeneity in the concentration profile is even more pronounced.This also explains the apparent lower yield stress when the L 0 = 200 nm and L 0 = 100 nm curves are compared.In the larger particle, the concentration at the boundary is larger at smaller SOCs, due to the heterogeneity in the mobility, as outlined above.In conclusion, Figure 4 shows that small particles and low C-rates are preferred in order to avoid high stresses and irreversible plastic deformations.Numerical efficiency.To show the capabilities of the adaptive solution algorithm of Subsection 3.2.4,we consider in Figure 5(a) the time step size τ n and the used order of the NDF multi-step procedure over two half cycles, i.e., one lithiation and one delithiation step with t end = 1.8 h in total, and in Figure 5(b) the refinement level for the spatial refinement after one lithiation at SOC = 0.92, comparing the lithium concentration distribution over the particle radius r.Here, we use a maximal order for the time adaptivity of three for Figure 5(a) and a minimal refinement level of three in Figure 5(b), respectively.During the first lithiation, there are three changes in the time step size and used time order: after starting with order one and switching to order two and three, the first plastic deformation arises at around t = 0.04 h at the first vertical gray reference line.Here, the time step sizes decrease and the used order goes down to two.After recovering to larger step sizes and orders again, the plastic deformation has ended at the second gray line at t = 0.12 h and the particle deforms elastically again.Then, a large time range with the maximal time step size τ max and maximal order is passed through.Shortly before the end of the first half cycle, the step sizes and order decrease again (third gray line).In this instant, plastic deformation occurs again, due to the tensile stresses at the particle surface.Changing the external lithium flux direction from charging to discharging is not trivial for the adaptive algorithm: the order decreases to one and the time step sizes drops over five orders of magnitude at the red reference line.After recovering from this event, one further plastification occurs shortly after the red reference line before the maximal time step size and the maximal order are reached again.At the time of t = 1.7 h, indicated with the fourth gray reference line, the next plastic deformation happens accompanied with a reduction of time step sizes and lower order.In total, there are two further plastifications during delithiation.In Figure 5(b), the focus is on the spatial adaptivity.We see the concentration distribution c over the particle radius r and the refinement level of the cells of our triangulation.It is clearly visible, that areas with larger concentration gradients have a higher refinement level due to the used gradient recovery error estimator.
In Table 3 we consider the average number of Newton iterations per time step, the assembling time for the Newton matrix and the right hand side of the linear equation system for one half cycle.We compare the results for the derived linearization and the AD technique in the 1D and the 2D simulation setup.More numerical results for the latter case are given in Subsection 4.2.2.Additionally, we compare for the AD 2D simulation case also these simulations: • without MPI parallelization, but with space and time adaptivity • without MPI parallelization and constant spatial discretization, but with time adaptivity • without MPI parallelization and constant spatial and temporal discretization Table 3: Comparison of number of time steps, average number of Newton iterations per time step (Av.Newton), assembling time for the Newton matrix and the right hand side of the linear system, between the derived linearization and AD techniques for the 1D and the 2D simulation case as well as other modern numerical techniques for the 2D case.For the second case, we choose the spatial discretization for the whole duration of the simulation to be just below the largest number of degrees of freedom (DOFs) in the adaptive simulation case.For the last situation with uniform temporal discretization, we choose 2 × 10 −7 h and BDF with order one for the time integration scheme. 2 × 10 −7 h is close, but just above the lowest time step size of the space and time adaptive simulation case after the initial time step size.These values are necessary to correctly record all relevant physical phenomena.Note that we count all Newton steps, also the Newton steps when the rate for reduction of the new Newton update is not fast enough or the tolerances of the adaptive algorithm are not fulfilled.Considering Figure 5(a), the majority of time steps in one lithiation cycle (from t = 0 to t = 0.9) are purely elastic and result in larger time step sizes.Due to the variable-step, variable-order time integration scheme, the maximal number of Newton iterations is limited.If the maximal number of Newton iterations is reached, the rate for the Newton update is too small or the tolerances of the adaptive algorithm are not achieved, the time step size is reduced.A smaller time step size leads to a lower number of Newton iterations.These points explain the low average number of Newton steps for our particular simulation setting.
The number of time steps and the average number of Newton iterations per time step are similar for the 1D and 2D simulations, respectively.However, there are significant differences in the assembling times: the linearization is twice as slow as AD for the 1D simulation and even ten times slower for the 2D simulation.This is a remarkable acceleration in the assembling times when solving for the Newton update.Not using the MPI parallelization leads to assembling time roughly four times slower, as expected.Considering the constant grid, we see a deceleration of over a factor of two, but also a slight reduction of the number of time steps.For the last case, we computed 100 time steps and extrapolated the total assembling times.The average number of Newton steps for these 100 steps is lower than in the other cases, also due to the fact that within these first 100 time steps no plastic deformation occurs.Considering a simulation without MPI parallelization, but with constant spatial and temporal discretization, would lead to a around 10 5 higher assembling time.In addition, there would be another factor of ten without using the AD technique, so in total, we have an speed up by a factor of roughly estimated 10 6 .All in all, Table 3 shows the huge acceleration for the assembling during the simulation with modern numerical techniques.
Physical results after nine half cycles.Before proceeding to investigate the 2D computational domain, we use the capabilities of the spatially and temporally adaptive algorithm together with the AD technique, to study multiple charging and discharging cycles.To be more precise, we consider five charging and four discharging cycles.Figure 6 shows that the viscoplastic case with the derived linearization or use of AD deliver identical results.The differences between the plastic and viscoplastic deformation are remarkable after nine half cycles, especially for the concentration in Figure 6(a).The plastic approach is almost identical to the elastic results, whereas the viscoplastic results are more similar to Figure 2. The similarity between the elastic and plastic results can be explained by the isotropic hardening, which leads to a large increase in the yield stress.Thus, after nine half cycles, occurring stresses inside SOC / -Fig.6: Numerical results for the elastic (Ela.), plastic (Pla.),(Vis.) and viscoplastic with AD (AD) approaches of the 1D radial symmetric case at t = 8.12 h over the particle radius r: concentration c in (a), equivalent plastic strain ε eq pl in (b) and tangential Cauchy stress σ ϕ in (c) as well as tangential Cauchy stress σ ϕ over SOC in (d).
the particle do not reach the yield stress, preventing the particle to deform plastically any further.This in turn leads to a more homogeneous distribution of stresses and thus no pile up of Li atoms close to the particle surface.The lower magnifier reveals also for the viscoplastic approach some difference for higher cycling numbers: a small wave occurs due to the change in the sign of the lithium flux for charging or discharging the active particle.This is important to notice because the small wave is already more pronounced for the tangential Cauchy stress in Figure 6(c) and is a candidate for possible further inhomogeneities.It should also be noted the higher magnitude in Figure 6(b) compared to Figure 2(b), indicating large extent of plastic deformation which could eventually lead to damage of the particle.Figure 6(d) shows the development of the tangential Cauchy stress at the particle surface over the SOC.For the elastic case there is no difference recognizable after the first half cycle and all further half cycles are identical.Also for the viscoplastic case, here computed with AD, only a small difference is visible at the lower left corner, see the magnifier in the bottom center.The second and all further lithiation cycles feature a higher compressive stress, even slightly higher than the elastic case, than in the first half cycle.This can be explained by the constant initial data.The magnifying glass at the right top shows the effect of the change from lithiation to delithiation with a short increase and followed decrease indicating an additional plastic deformation.A larger difference between the elastic and viscoplastic approach at lower SOC is observable at the end of the discharging cycle at the left top corner with the magnifier at the left center.The plastic ansatz with isotropic hardening, however, shows the biggest difference in each cycling.As already discussed, the tangential Cauchy stresses tend from cycle to cycle to the values of the elastic case indicated with the cyan arrows.It should be noted that the small wave of Figure 6(a) and (c) is not visible in this representation.

2D Quarter Ellipse
We now evaluate numerical results for the 2D quarter ellipse simulation setup.Here, we adapt the initial time step size to τ 0 = 1 × 10 −8 and a minimal refinement level to three.Further, we set θ c = 0.005 and τ max = 5 × 10 −3 .With the capabilities of the adaptive space and time solution algorithm as well as the presented AD technique, we consider now the effects of half-axes of different lengths on plastic deformation and concentration distribution.
Physical results after one half cycle.Figure 7 shows the numerical results at the end of one half cycle with t end = 0.9 h and SOC = 0.92.The number of DOFs are in a range of [320,121472] and the total computation time is less than 15 minutes.In all subfigures, the adaptive mesh is visible in the background, indicating higher refinement levels at the external surface near the smaller half-axis.In Figure 7(a), the concentration profile is displayed in a range of [0.91, 0.94].It is noticeable that the gradient in the concentration is stronger at the particle surface on the short half-axis.Here, a larger area with lower concentration values in blue can be located with a small but steeply growing area near the particle surface.This property is comparable to the 1D simulation results with a higher increase near the particle surface in Figure 2(a) or in Figure 6(a).
Having a look on the scalar equivalent plastic strain ε eq pl in Figure 7(b) the area of plastic deformation can be plastic deformation occurs at the particle surface of the smaller half-axis.Large equivalent plastic strains of up to 22 % occur.Figure 7(c) shows the von Mises stress in the general plane state A significant inhomogeneity of the von Mises stress distribution in the particle is identifiable.Whereas the stress values at the larger half-axis are lower than 0.45 GPa, the von Mises stresses are more than 1.5 times higher at the smaller half-axis and feature two peaks near the particle surface: one due to tensile stresses because of the plastic deformation next to the surface and one due to compressive stresses a little further inside of the particle, both again like in the 1D case in Figure 2(c) or in Figure 6(c).The strong change of stress gradients near the surface of the smaller half-axis may appear unusual.However, as outlined, the change of the gradients results from the change in sign of the tangential stresses from tensile at the particle surface to compressive in the particle interior close to the surface, c.f. Figure 2(c).As the von Mises stress, c.f. Equation (46), is the absolute value of the stress deviator, this change results in the two maxima close together.This inhomogeneity of the stress distribution between compressive and tensile stresses is especially important for further investigations on particle damage or fracture.
However, the finding of the higher stresses and plastic deformation at the smaller half-axis are in contrast to the numerical outputs in Figure 4(b).There, larger particle size results in larger stresses, so it would be reasonable to feature higher stresses at the longer half-axis.Though, the longer half-axis is responsible for lower concentration values in the particle center.On the smaller half-axis, the diffusion would create higher concentrations in the particle center.Thus, the influence of the lower concentration values at the longer half-axis is responsible for higher concentration gradients at the smaller half-axis resulting in higher stress and plastic deformation.In this context it is important to recall that a constant external lithium flux is used at the particle surface, which favors this behavior.

Conclusion and Outlook
We conclude our work with a summary and an outlook.

Conclusion
Within this work, a large deformation, chemo-elasto-plastic model to predict the diffusion-deformation behavior of aSi anode particles within Li-Ion batteries is presented.The used plasticity model is in close analogy to the model presented in [7].As an extension, we specify our theory for both rate-dependent and rate-independent plastic deformations, see Subsection 2.5, and compare numerical results of both theories after multiple charging and discharging cycles.The chemical model relies on [12], where a measured experimental OCV curve is used to model the chemical contribution to the Helmholtz free energy.The derived model is implemented in an efficient finite element scheme, first presented in [11], relying on space and time adaptive algorithms and parallelization [21] as well as modern numerical techniques like AD.
For both plastic models, a return mapping algorithm [31] is used in the context of static condensation, allowing the evaluation of the plastic deformations at the finite element integration point level instead of treating them as additional degrees of freedom [15].We present in Subsection 3.2.1 a projection onto the set of admissible stresses inspired by [26], which can be given explicitly in the case of the rate-independent model with linear hardening, whereas in the viscoplastic case the projector is implicit, due to the non-linearity introduced by the evolution of the equivalent plastic strains.The linearization of these projectors, necessary for computing Newton updates in the nonlinear solution scheme, is approximated in Subsection 3.2.4inspired by [7].In addition, we apply AD techniques to circumvent the necessity of tedious analytic computations and assembling operations and compare the numeric performance of both approaches.We incorporate our DAE system into an existing efficient space and time adaptive solver, presented in Section 3.
The numerical results in Subsection 4.2 for a given set of parameters show a heterogeneous concentration distribution within the particle, when plastic models are considered.We attribute this to lower stresses inside the particle, which in turn affect the mobility of Li atoms and lead to a pile up at the particle's surface.In our studies, plastic deformations are limited to the outer parts of the particle and result in an eigenstrain that leads to tensile stresses at the particle surface during charging, which is not observed in the elastic model.Both the concentration pile up and the plastic deformation can be mitigated by decreasing the particle size and or decreasing the charging rate.
Comparing the plastic and viscoplastic model in Figure 6, differences become visible after multiple charging and discharging cycles.On the one hand, the plastic model hardens isotropic and after several cycles the occurring stresses do not reach the increased yield limit, preventing further plastification.On the other hand, the viscoplastic model does show an increased amount of plastic deformations, which results in a different stress and concentration distribution.As a result of Section 4, the differences are more pronounced the more cycles are considered, which seems not to be investigated before and affects battery performance and lifetime.This could be achieved by the use of several modern numerical techniques.
Investigating the performance of the approximation of the projector's linearization in comparison with the AD scheme in Table 3, we conclude that the AD scheme is way more efficient, due to the fact that the assembling of the Newton matrix is done simultaneously with the residual.This leads to decreased computation times.In addition, the use of AD circumvents the tedious analytical derivation as well as the fault prone implementation of the linearization.In total, with the application of efficient tools like parallelization, adaptive spatial and temporal discretization and AD on our higher order finite element approach, we have reduced our computational effort for the assembling significantly.Adding the plastic part of the deformation gradient as solution variable to the set of equations like in [15] would increase the computational time even more, especially for the purely two-dimensional computational consideration.
When studying a two-dimensional problem in Subsection 4.2.2, an asymmetry in the concentration and plastic strain distribution is observed, which we attribute to the ellipsoidal particle shape and the resulting asymmetric concentration distribution.
To conclude, our study indicates that small particle sizes with a spherical shape result in a smaller build up of stresses and is therefore desirable.Moreover, semiaxes with different lengths should be avoided to obtain homogeneously distributed mechanical stresses.Regarding material behavior, the considered isotropic hardening mechanism is favorable, as less plastic strains build up after multiple charging and discharging cycles.In addition, this would prevent the sharp concentration gradients observed in our study and thus less interference with battery performance is to be expected.A further combination with the viscoplastic approach could be a possible extension.

Outlook
To further extend the theory originating from this work, numerous paths can be taken.We consider the following ones to be especially interesting: • Due to the usage of the chemical potential as solution variable and the corresponding code structure, first introduced for the phase-separating material lithium iron phosphate in [11], we can add a Cahn-Hilliard approach for the phase-separation phenomena of aSi which is observed by [4].
• Previous research pointed towards the usage of nanowires as electrode geometry.The reason lies in the observation that lower electrode sizes decrease occurring stresses and plastic deformations.This is in agreement with the results of decreasing stresses for smaller particle radii.Still, a study of various nanowire shaped electrode geometries is a reasonable application for the derived model and the proposed implementation scheme.• Fracture and SEI formation typically occur in aSi anode particles [15].Both have been considered previously, see, e.g.[15,28,63], and should be included in future extensions of this work, to better capture the physical behavior of the anode particles.• A paper published recently by several of the authors [12], introduced an efficient scheme to compute particle deformation behavior when contact between particles or walls is introduced.A study of the elasto-plastic material behavior derived in this work, in combination with the obstacle, could provide further insight into the mechanics of anode materials during charging and discharging.• An integration of the Butler-Volmer condition into the external boundary condition for the lithium flux is of further interest to consider additional surface boundary reactions.• Finally, a more robust and scalable solver can be developed to provide additional speedup compared to the currently used LU-decomposition when using MPI parallelization [21].Further consideration can also be given to the matrix-free concept of deal.ii[54,64] and other solvers like PARDISO [65].

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

B Derivation of First Piola-Kirchhoff Tensor
The first Piola-Kirchhoff stress tensor is defined as the derivative of the free energy density by the deformation gradient F, compare [36, Section 6.1], as The chain rule for the second equality in Equation ( 47) can be computed as stated, since ρ 0 ψ el = Now, we have In total, the first Piola-Kirchhoff tensor P is computed via

C Linearization of Projector
We estimate the formulation of the exact Jacobian for the rate-independent plastic approach as in Di Leo [20,Appendix C.6] with the formal linearization of the projector P Π of Equation ( 29) around E tri el [26,30] I Π c n+1 , ∇ 0 u n+1 , F n pl , ε eq,n pl Then, the consistent algorithmic modulus is given as el ≤ σ Y (c n+1 ) + γ iso ε eq,n pl , C pl , C E tri,dev el > σ Y (c n+1 ) + γ iso ε eq,n pl . (66) With the expressions for M n+1 in the projector P Π of Equation ( 29 ], the linearization of P Π around E tri el can be expressed by Equation (42).A similar calculation leads for the rate-dependent case to Figure 8 shows the 1D radial symmetric results of the tangential Cauchy stress σ φ over the particle radius r for the elastic, plastic and viscoplastic case at SOC = 0.1.The area at the particle surface is plastically deformed for the non elastic approaches and is qualitatively comparable to numerical results in Fig. 4(c) of [61] and Fig. 5(d) of [6].With a reduction of the yield stress σ Y (c), the plastically deformed range would be larger.

E Comparison of Logarithmic Strain vs. Green-St-Venant Strain
Comparing the logarithmic strain approach and the Green-St-Venant strain approach as discussed in Subsection 2.1 and Subsection 4.2.1, both approaches result in equal outputs, displayed in Figure 9.This corresponds to the findings in [62].Pla.Vis.
Fig. 8: Numerical results for the elastic (Ela.), plastic (Pla.) and viscoplastic (Vis.)approaches of the 1D radial symmetric case at SOC = 0.1 for the elastic, plastic and viscoplastic approach of the tangential Cauchy stress σ φ over the particle radius r.

F Comparison of Radial Deformation Gradient
In this section, we take a closer look at the different components of the three different parts of the deformation gradient F pl , F el and F ch .At the start of the simulation, no displacement gradient is present and the total deformation gradient tensor is the identity tensor.In Figure 10, we consider the radial part of the each deformation gradient tensor, respectively.It can be clearly seen that the chemical part makes the largest contribution to the deformation, followed at a large distance by the plastic and the elastic parts which stay close to one.This finding justifies the used approach of a linear elastic theory.
el and the first and second Lamé constants are λ = 2Gν/ (1 − 2ν) and G = E/ 2 (1 + ν) , depending on the Young's modulus E and Poisson's ratio ν of the host material.

Fig. 1 :
Fig. 1: Computational domains in 1D in (a) and 2D in (b) used for the numerical simulations.
) and the magnitude of the equivalent plastic strain increases to 4 %.For the tangential Cauchy

Fig. 3 :
Fig.3: Influence of the different plasticity approaches in (a) and of the maximal yield stress of the viscoplasticity approach in (b) for the tangential Cauchy stress σ ϕ at the particle surface r = 1.0.

Fig. 4 :
Fig. 4: Tangential Cauchy stress σ ϕ over SOC at the particle surface r = 1.0 for varying C-rate in (a) and particle size in (b).

Fig. 5 :
Fig. 5: Advantages of the adaptive solution algorithm: time step size τ n over simulation time t for two half cycles (one lithiation and one delithiation) in (a) and concentration c as well as refinement level of the spatial discretization over the particle radius r at SOC = 0.92 in (b).

Fig. 7 :
Fig. 7: 2D quarter ellipse after one lithiation at SOC = 0.92 for concentration c in (a), equivalent plastic strain ε eq pl in (b) and von Mises stress σ vM in (c) using AD techniques.

Fig. 10 :
Fig. 10: Comparison of the radial component of the plastic, elastic and chemical part of the deformation gradient tensor F pl , F el and F ch over the particle radius r at SOC = 0.92.
composed of solutions on former time steps y n , . . ., y n−k and a constant α kn > 0 dependent on the chosen order k n at time t n [50, Section 2.3].The vector f depends explicitly on the time t due to the time-dependent Neumann boundary condition N ext due to the change of the sign for charging and discharging, respectively.After computing the new discrete solution y n+1 we can update F n
[60,emented in C++.Further, we use the interface to the Trilinos library[55, Version 12.8.1]and the UMF-PACK package[60, Version 5.7.8]for the LU-decomposition for solving the linear equation systems.A desktop computer with 64 GB RAM, Intel i5-9500 CPU, GCC compiler version 10.5 and the operating system Ubuntu 20.04.6 LTS is used as working machine.Furthermore, OpenMP Version 4.5 is applied for shared memory parallelization for assembling the Newton matrix, residuals and spatial estimates and message passing interface (MPI) parallelization with four MPI-jobs for the 2D simulations with Open MPI 4.0.3.Unless otherwise stated, we choose for the space and time adaptive algorithm tolerances of RelTol t reduction of the last two dimensions of the fourth order tensor A and the second order tensor B ∂ □ partial derivative with respect to □ △ n □ difference of next and current time step: □ n+1 − □ n or □ n+1 − □ n