A multi-field model for charging and discharging of lithium-ion battery electrodes

An electrochemical–thermomechanical model for the description of charging and discharging processes in lithium electrodes is presented. Multi-physics coupling is achieved through the constitutive relations, obtained within a consistent thermodynamic framework based on the definition of the free energy density, sum of distinct contributions from different physics. The system is characterized by finite kinematics, under the assumption of locality of deformation, and the deformation gradient is decomposed into the product of elastic and inelastic parts. Specifically, a Taylor series expansion is used to approximate the inelastic deformation due to ion intercalation. The elastic part can be described alternatively by two finite kinematics models of neo-Hookean elasticity, and a Maxwell-type viscoelastic model accounts for time-dependent mechanical aspects. The model is implemented into a finite element code that uses B-spline basis functions. We illustrate the features of the model by means of selects examples, showing that chemo-mechanical interaction affects the equilibrium concentrations of the phases. The model captures the fundamental aspects of the anode charging and discharging processes.

Making use of the classical Boltzmann continuum assumptions [50], this study describes a multi-physics numerical model able to simulate the electrochemical process of charging and discharging of an electrode. The coupled constitutive relations are obtained through the thermodynamical approach of Coleman and Noll [12], based on the assumption that the local state of a system is entirely described by the deformation, the entropy, and a set of suitable internal variables defined within a small neighborhood of a material point. The proposed model extends previous models by [2,8,22,39,65] and others, by coupling the full electrochemical-thermomechanical process of charging and discharging with large reversible and irreversible mechanical deformations.
In general, the performance of the electrodes during charge-induced ionic flow is size, geometry, and structure dependent [7,19,29,58]. The bulk material may be subjected to structural changes, e.g., a transformation from crystalline to amorphous has been reported for silicon anodes [23,36]. Geometrical aspects coupled with the electrode design have been investigated experimentally, revealing that hollow micron-sized silicon structures are advantageous high-performance particles [64] and that solid amorphous silicon rods do not crack [10]. In common commercial electrode materials, the lithium ion content induces material softening or hardeningan effect which has been observed experimentally [41] and predicted in first principle density functional theory computations [37,42]. Moreover, the formulation of high-fidelity models for specific electrode compositions is an active area of research, see [8,17,31] among others.
Most approaches used to model batteries account for small deformations and linear elasticity, cf. [8,22,57]. The description is comparatively easy and is adequate for many applications, but new developments, such as silicon anodes, require to account for considerable changes of shape and volume. A simple coupled largedeformation model for anode charging, neglecting chemical decomposition and phase separation, was proposed in [13]. A continuum-level theory, coupling diffusion and finite deformations, was proposed in [2] and later solved in a simplified manner by standard finite elements (FEM), [26,38]. Typical cathode materials, such as lithium manganate, have been investigated for stress-related phase separation [22,66] also using a nonlocal diffusion model [62]. Advanced FEM simulations indicated that the formation of a core-shell structure able to inhibit further lithiation depends strongly on the choice of kinetic parameters and on the state of charge [67]. Cracking of the electrode material has also been investigated in [47], although it has not been observed experimentally [10].
In abstract terms, charging and discharging of a lithium-ion battery electrode result from particle exchange between the anode material A (e. g., silicon or graphite) and the electrolyte (e. g., LiPF 6 salt), The physical processes that affect the electrochemical reaction at the electrode are: ion flux and diffusion, electrical and thermal conduction, swelling and mechanical deformations, and electro-mechanical forces associated with the electrostatic field. The corresponding field equations are illustrated in Sect. 2. The thermodynamical consistent derivation of the general constitutive equations for the coupled electrochemical-thermomechanical process of charging and discharging is outlined in Sect. 3. In Sect. 4, we provide the constitutive model for lithium electrodes by defining the specific Helmholtz free energy density and using its partial derivatives to obtain the mechanical stresses, the thermal entropy, the electric induction, and the chemical potential of the Li + enriched electrode material. The properties of the coupled model are illustrated in Sect. 5, comparing the proposed approach with other models from literature. In Sect. 6, we provide FEM examples demonstrating the interactions of the different physical fields. The features of the proposed electrochemical-thermomechanical model are summarized in Sect. 7, while Appendix provides some fundamental relations.

Field equations
A point of a continuum in its material configuration B 0 ⊂ IR 3 is denoted by X. An arbitrary deformation in a time interval [0, t end ] moves the material point to a new position x(X, t) in the current configuration B t ⊂ IR 3 . The corresponding deformation mapping χ(X, t) : B 0 × [0, t end ] → B t is uniquely described by its material gradient, The local volume change at x is measured by the Jacobian determinant J = det F, which assigns an infinitesimal volume element dV of its material configuration to its spatial configuration dv, dv = J dV .
The field equations that govern the mass, concentration, elastic, electrostatic, and thermal fields are summarized in Table 1 in the common spatial form. Corresponding quantities of the material configuration will be written with capital letters. The complete notation is listed in Table 3 of Appendix, where also the transformation rules between the configurations are provided. In the following, we illustrate the equations of each of the physical processes in turn.

Partial and total mass balance
We consider a general multi-component system consisting of n species. The mass density ρ k (x, t) of the k-th species represents the mass m k (x, t) per current volume v, and the sum gives the total mass density ρ The electrode charging and discharging implies mass transfer, described by a rate per unit of current volume h. The evolution of the local density field follows from Eq. (2), which provides the conservation of the mass J = ρ 0 /ρ in the special case of h = 0.

Concentration and particle diffusion
The mass fraction of the k-th component ρ k is defined as and is further referred to as concentration. By definition, it holds n k=1 c k = 1. The global mass balance for each of the n individual species can be formulated by introducing the inward mass flux j k = − j k · n, where n is the unit vector normal to the boundary. The total mass change in Eq. (3) accounts for the single component mass transfer through the density rate h k , induced by the chemical reactions occurring inside the body (see below). According to the localization theorem (see Appendix, (66)), the total mass rate contributes to the local mass flux, Eq. (4), through the concentration c k (10).
The diffusive fluxes are assumed to depend on the spatial gradient of the molar chemical potential μ k through a second-order tensor of particle mobilities m ik , The molar chemical potential μ k comprises the free energy gain for adding one mole of k-th component to the system and does not depend on the configuration. For the special case of mass conservation, under the assumptions of isotropic mobility, absence of chemical reactions, and constituents of approximately equal molar mass, the temporal evolution of the concentration simplifies to where m ik = ρ D ik c i c k /θ (T ) is the degenerate mobility and m ik = m ik 1. The temperature dependent diffusivity D i j = D i j (T ) is typically modeled with the Arrhenius law. A straightforward calculation using the relations of the Appendix shows that Eq. (12) can be expressed in the same way with the material gradient ∇ X μ.

Chemical reactions
In addition to the pure diffusion of the species, we consider here a typical binary chemical reaction where the reactant consists of n A and n B moles of molecules of the species A and B, respectively, which react to produce n C and n D moles of molecules of species C and D, respectively. The reaction rate J of the chemical reaction gives information on the kinetics of the chemical reaction. Generally, we decompose the reaction into r stoichiometrically independent reaction pathways with rates J j , j ∈ {1, 2, . . . , r }. The driving forces of the process, i. e., the reaction affinities A j , are assumed to derive from a variational principle where the chemical species are subject to constant supply rates [18]. The chemical affinity A j of the j-th reaction is then defined as This contributes to the local mass evolution of the k-th component (last term in Eq. (4)). Since mass is conserved in each separate chemical reaction, it holds r k=1 ν k j = 0.

Linear and angular momentum balance
The deformation of a system follows Newton's laws of motion. For an overall Cauchy stress field σ (x, t), a mean velocity v(x, t) and a body force per unit volume b, Newton's law reduces to the local form of the linear momentum balance (5) and of the angular momentum balance (6). A surface traction t = σ · n may be assigned at the boundary. Here, the dynamics of n-component mixtures, i.e., convective effects and inertia of the different components, will be neglected, cf. [3]. For the material description, we employ the first Piola-Kirchhoff stress tensor P The angular momentum balance is then satisfied through the symmetry of the product P F = F P .

Electrostatics
The spatial electric field e is given by the spatial gradient of an electric potential φ e , e (x, t) = − ∇ x φ e (x, t), and the electric induction (or electric displacement) d is proportional to the electric field through the vacuum permittivity ε 0 . For a dielectric material, a common constitutive relation assumes d = ε 0 ε r e, where ε r is a concentration and temperature-dependent relative permittivity. For brevity, we write ε e (c, T ) = ε 0 ε r (c, T ). The governing equation for electrostatics is ∇ x ×e = 0 and, for a local charge density ρ e z, the electric induction must satisfy the equilibrium Eq. (7). Here, the mixture current fluxes between components are neglected, cf. [56].
For later reference, we formulate the electrostatic relations with the corresponding material fields, E = F e and D = J F −1 d. Then, the material electric induction follows as

Heat conduction
Heat transfer is assumed to be governed by Fourier's Eq. (8). The material thermal flux H T relates to its spatial counterpart h T by H T = J F −1 h T . Commonly the spatial thermal flux is assumed to be linearly dependent on the gradient of the temperature through a spatial second-order tensor of thermal conductivities k T . The constitutive relation in material form follows through a pull-back of the thermal conductivities, Remark: By accounting for a few reasonable simplificative assumptions, the previous field equations have been formulated in a fully decoupled form and independently one from the other. Coupling between species diffusion, mechanics, electrostatics, and heat conduction is introduced by means of the constitutive relations, which will be derived in the following section under suitable thermodynamical considerations.

Thermodynamic derivation of constitutive relations
The constitutive equations characterizing the material properties must comply with the basic laws of thermodynamics, namely the conservation of mass, the energy balance, and the entropy principle. To this aim, we collect the energy contributions of all physical processes at first and define the body's internal energy. Next, we state the entropy principle, quantify the fluxes, and derive the corresponding form of the dissipation inequality.
As described in detail in [57], the formal derivation of the electro-chemo-thermo-mechanical model relies on the thermodynamics of irreversible processes and assumes local equilibrium conditions within a small volume of material at position X [14]. The internal energy U changes with contributions of heat supply Q, heat flux H, power of deformation P :Ḟ, mass diffusion, chemical reactions, and the presence of an electric field E ·Ḋ. Thus, the rateU in the material configuration can be expressed aṡ Both reaction affinities A j and chemical potentials μ k are referred to the material volume. The mass balance comes into play through the definition of the concentration fields (10), whereas gravity forces acting on the diffusive fluxes are neglected. Moreover, the electric term is included with a negative sign because it follows from the original term E ·Ḋ by a Legendre transform. The local rate of the entropy per unit material volume isṠ = − ∇ X · H S + Π S , with entropy flux H S and nonnegative entropy production Π S . Here we follow the classical Duhem's assumption and identify the entropy flux with the heat flux per unit temperature H T /T ; heat supply contributes to entropy production. Consequently, we obtain for the entropy production, where the last term expresses the dissipation by heat conduction, and the condition − 1 T H T · ∇ X T ≥ 0 determines the direction of heat flow. Further internal dissipation, D int ≥ 0, may arise as a consequence of material damage, viscous flow, and Lorentz forces, for example. We are aware that in the case of mixtures, the Duhem's assumption on the entropy flux is too restrictive. Nonetheless, for common continuum mechanical systems, more general formulations for H S result in the same thermodynamic restrictions for the constitutive equations, cf. [49].
To express the dissipation inequality (19) in terms of process energy contributions, we make use of an alternative thermodynamic potential, namely the Helmholtz free energy density Ψ = U − T S, obtained from (18) through a Legendre transform. Inserting its rate and (18) into (19) gives the dissipation inequality Note that the mechanical power includes the elastic power of deformation, P :Ḟ, and inelastic contributions that may result from internal dissipative processes, e. g., plasticity, and viscosity. The thermodynamic state of the material point is then entirely defined by F, T , c k , J j , E, and internal variables Z, which account for such dissipative phenomena. The Helmholtz free energy density is thus a function of the form Ψ = Ψ (F, T, c k , ∇c k , J j , E, Z). Taking the differential giveṡ where we introduced the thermodynamic forces conjugate to Z, The chemical affinities A j can also be derived from dissipation potentials subjected to the stoichiometric constraints, i.e., they are work conjugate to the reaction-path rates J j , cf. [18]. Additionally, we split the mechanical stress into an equilibrium part P eq , depending on the variables of state only, and a rate-dependent dissipative contribution P v , also depending on the rate of deformationḞ, Inserting these relations into the dissipation inequality (20) gives where we denote δ c • = ∂ c • − ∇ · (∂ ∇c •). Since inequality (24) must hold for any admissible process, the terms in brackets must vanish, leading to the equations of state (or constitutive relations) for the equilibrium stress, entropy, concentration, and electric induction: The last expression follows from mass conservation, n k=1ċ k = 0, and a restatement of Ψ . The surviving terms define the dissipation as which constrains the constitutive law for the heat flux H T and the set of internal variables Z. An admissible choice is Duhamel's law of heat conduction (17), which defines the heat to flow for a negative temperature gradient, Eq. (17). The general framework is completed with the definition of the viscous dissipation and the kinetic relations for Z.

Constitutive model for Lithium electrodes
Here we will present the constitutive model for the electrochemical-thermomechanical processes of anode charging in lithium batteries. In particular, we reduce the process to a binary diffusion system without additional chemical reactions and refer for the latter to [56]. The normalized concentration of Li + ions is denoted as c = c 1 , whereas the complementary concentration of the anode's material is c 2 = 1 − c 1 .

Kinematics
We start by specifying the kinematics of deformation and introduce the multiplicative decomposition of the deformation gradient (1) into elastic and (several) inelastic components, The multiplicative decomposition (27) is a convenient mathematical representation of configurational changes of a system undergoing multi-physics processes in large deformations, [16]. It was introduced in [24,43,48] for thermoelasticity, viscoelasticity, and phenomenological elastoplasticity, respectively. At the electrode, the inelastic deformation results from geometrical changes due to the Li + intercalation and swelling F s , viscous distortion F v , thermal expansion F th , and the active electric deformation F act . We remark that the inelastic parts of the deformation take place locally and result in a non-compatible intermediate configuration. The elastic part of the deformation gradient F e is related to the passive response of the material and relaxes the continuum from the intermediate to the current deformed configuration where equilibrium and compatibility conditions are fully satisfied, Fig. 1.
The swelling as a consequence of the Li + intercalation is unknown but assumed to be a purely volumetric deformation. Further inelastic mechanical deformations are identified with the isochoric viscous deformation components F v . Local deformations induced by the electric field are very small, cf. Sect. 5.3, which allows us to set F act to unity. Thermal expansions can be stated in the form F th = α th 1 + i a i N i ⊗ N i where α th , a i are temperature-dependent functions, cf. [48]. Because of the amorphous anode material, we assume an isotropic expansion with a i = 0 and an expansion coefficient which depends on T and c, α th (T, c) ≥ 0. With these assumptions at hand, the multiplicative decomposition of the deformation gradient (27) can be restated as The relation between the inelastic swelling J i and the concentration and phase segregation of intercalated lithium ions, i.e., the function J i (c) is unknown and must be determined. A straightforward ansatz is to assume a direct proportionality between volumetric expansion and c, [22,30], as commonly done in small strain elasticity models, [21,28,63]. The finite deformation model proposed by Anand [2] suggests that only known relation is ∂ J i /∂c ≥ 0. Using a virtual-power approach, this conditions leads directly toJ i = Ωċ, where the constant Ω is related to the atomic volume of Li + . Xu and coworkers extended this model in order to perform numerical calculations [66,67] and concluded that the inelastic swelling J i must be compensated by an elastic swelling J e to obtain a concentration-independent overall deformation as reported in [2]. This condition leads to a dependent elastic deformation J e (J i ) and Here we propose a different model and assume J i (c) to be an unknown function. Then, a Taylor expansion at c 0 gives where c 0 is a reference value. The elastic volumetric deformation will relax the continuum from the intermediate configuration to the fully compatible current configuration and thus we state: Next, we collect all the contributions to the free energy involved in the process. For ease of readability and implementation, we chose to formulate the free energy density as a function of independent variables F and J , cf. [11,20]. Moreover, we define the isochoric part of the elastic deformation gradient asF

Elastic energy
The sensitivity of a battery to mechanical deformation depends strongly on the material of the electrode. For example, the specific capacity of silicon is one order of magnitude larger than the one of graphite. The excellent capacity of a silicon anode is associated with a > 300% volumetric expansion during charging [9], relative to the 6-10% volume expansion observed for in graphite anode [1].
We introduce a hyperelastic material model characterized by an elastic free energy Ψ e independent of any dissipative internal process. Commonly, Ψ e is decomposed additively in volumetric and isochoric parts; thus, we formulate a neo-Hookean energy density, extended to the compressible range, of the form where K is bulk modulus and G a shear modulus. Both moduli may depend on Li + concentration and temperature. An alternative form of the hyperelastic material model, proven to be convenient in numerical applications [44][45][46], is the following The main advantage of the model (33) is its simplicity: It saves the decomposition of the deformation gradient and has a direct analogy to linear elasticity. Starting from the elastic parameters of the two species, the equivalent elastic parameters to be used in (32) and (33), we apply a Vergard's mixture rule

Thermal energy
As is any chemical system, the battery state depends on the temperature. The temperature rise during charging is a slow process. Because we are expecting a small effect of the anode's design, here we consider the temperature as driven by Joule heating but not affected by the deformation. Thus, the thermo-mechanical coupling can be modeled with the classical structural entropy potential whereas the potential for the thermal entropy is

Chemical energy
For the mixing energy of the diffusion system, we assume a constant pressure which allows us to formulate the free energy in equivalence to the Gibbs free energy. The configurational free energy density Ψ con for a binary mixture with two equilibrium phases at the concentration c α (α-phase) and c β (β-phase) is modeled with the common Margules form Here g 1 and g 2 are Gibbs free energy densities of both components. The third term accounts for energy contributions from the entropy of mixing. According to the classical Lewis-Randall ideal-solution model, the entropy contribution is typically weighted by a temperature-dependent material function θ(T ), e.g., θ(T ) = RT y ref , where y ref is reference molar concentration. For non-ideal mixtures, the energy must include an excess energy contribution that describes the decomposition of a non-ideal mixture into distinct phases. Here, the excess energy is chosen to be of Porter type [see [33]], with the Flory-Huggins interaction parameterχ . Forχ < 2, the mixing energy (38) is convex and corresponds to the energy of a homogeneous solution. Phase separation will be observed only forχ > 2, when the configurational energy turns into a double-well function, i.e., it has two relative minima and a concave (spinodal) region in between. The concentrations c α and c β are the binodal points where ∂ c Ψ con (c α ) = ∂ c Ψ con (c β ). The situation is illustrated for different values ofχ in Fig. 2.
Furthermore, the coexistence of phases causes additional energy contributions. The interfacial free energy density accounts for the nonlocal effect of surface tension with interface tension parameter κ > 0.

Electric energy contribution
For the effect of the electric induction on the material, we introduce the active part of the Helmholtz free energy. With E = F e and assuming a standard dielectric with permittivity ε e , it reads The thermo-electric coupling reduces to resistive heat generation; other thermoelectric effects like described in [40] are not expected here. The electric power determines Joule heating. For a given electric field E, it has the potential where R m is the material's resistivity and α ≤ 1 a conversion factor.

Inelasticity and viscosity
In the charging process of the electrode, some inelastic mechanical deformations may occur, cf. [10]. Within a variational framework, the existence of a free energy potential is assumed such that the inelastic driving forces are derived as work conjugate to the inelastic part of the deformation, e.g., and further forces are conjugate to the internal variables, Eq. (22). Additionally, a rate-dependent potential Φ may be introduced to characterize the rate-dependent dissipative processes. The rate of the inelastic deformations is subjected to the kinematic restrictions imposed by a flow rule, and for the evolution of the internal variables, suitable kinetic equations must be supplied. This technique leads to variational constitutive relations, as described in [6,34,59] for plasticity and viscoelasticity. From the typical cyclic loading of an electrode, we can interpret the inelastic mechanical deformation as viscous, i.e., the material state depends on the rate of deformationḞ. The equilibrium stress state is defined as P eq = P(Ḟ = 0). The viscous stress is derived from rate-dependent free energy contribution Ψ vis such that Here we choose viscosity to be volume preserving Newtonian. Specifically, we refer to the common generalized Maxwell model composed of N spring-dashpot elements in parallel, each with shear modulus G j , viscosity η j , and relaxation time τ j = η j /G j . The contribution to the body's free energy density is then where for each of the Maxwell elements, the rate-depended energy density is In (45), F v j denotes the viscous part of the deformation gradient of the j-th Maxwell element which follows the evolution equation The total F v is the weighted sum of all elements in time and leads to the solution We remark that our model corresponds to the classical hereditary integral theory of linear viscoelasticity [52] and can conveniently be solved with the help of a recursive scheme. For simplicity, we assume the moduli G j and η j to depend on c and T in the same way, i.e., τ j is independent of concentration and temperature.

Summary
The total Helmholtz free energy density of the anode material has the following expression:

Mechanical stresses
The first Piola-Kirchhoff stress tensor follows from the constitutive relation (25) 1 as the work conjugate to the deformations gradient. In equilibrium, it is and additional rate-dependent stress contributions may apply, Eq. (23). With decomposition (27), the isochoric expressionF e and involving the corresponding derivatives using the chain rule, we obtain for P e The resulting equilibrium stress of Eq. (49) is the pulled-back elastic stress. Defining a conjugate stress P = ∂Ψ iso /∂F e and a pressure term p = ∂Ψ vol /∂ J e we re-state it as which illustrates that any ansatz for the hyperelastic energy density just enters the expressions forP and/or p. For the specific free energy model (32), we evaluate the first Piola-Kirchhoff stress tensor For model (33), the first term of (52) changes, and we obtain the mechanical stress Clearly, the coupling of diffusion, temperature, electric field, and elasticity manifest itself in the different stress contributions (and also in the field-dependent material moduli). The first lines of (52) and (53) include the thermo-mechanical stresses; next follows the viscous stress contribution P v , and the last terms account for the presence of an electric field, which is usually referred to as the Maxwell stress and its heating contribution. The Cauchy stress σ in the current configuration can be obtained by a push-forward operation, according to Eq. (15).

Thermal field
The entropy field (25) 2 is conjugate to the temperature T with the result The fact that S = 0 for T → 0 needs to be considered in the constants α th (T ), G(T ), . . . but is of no relevance here.

Chemical potential
The chemical potential of a species follows from (25) 3 . It expresses the free energy change for adding or removing one mole of a component to or from the system; thus, it is the same for any configuration. With energy (48), we obtain where we use (52) (or alternatively (53)) for the derivative of the elastic energy. Moreover, ∂ K /∂c is the derivative of the overall elastic bulk modulus function (34) with respect to the corresponding concentration; the same holds for G, α th , . . .. These terms vanish for constant elastic parameters. Obviously, the first two lines in (55) include the coupling of mass diffusion and the thermoelastic field; the third line comes from the evolution of the chemical field and the thermo-electro-chemical coupling gives the fourth line.

Electric field contribution
The constitutive Eq. (25) 4 provides the electric induction as and corresponds to (16) for isothermal conditions.

Properties of the coupled model
This subsection presents characteristic features of our coupled electrochemical-thermomechanical model. For clarity, univariate parameter studies are performed and, because a thermo-mechanical material behavior is well investigated, the focus is set here on the chemo-mechanical and electro-mechanical interaction. Additionally, we provide a comparison to the finite deformation model of Anand [2] and Xu [66].

Chemical energy, phase decomposition and mechanical properties
We consider two material states with bulk moduli K 0 = K (c = 0) and K 1 = K (c = 1). Depending on the electrode material, these values may be rather similar, e.g., K 1 ≈ 0.9 K 0 for a typical cathode material like Li x Mn 2 O 4 [25], or quite different, e.g., K 1 ≈ 0.05 K 0 for the Li-Si mixture in a silicon anode [27,42]. For both situations, the influence on the equilibrium concentrations c eq α and c eq β is studied. The Helmholtz free energy density (48) follows here from Eqs. (32), (34) and (38) with g 1 = g 2 = 0 and a material-dependent normalization by θ in the dimensionless form For further investigation, we set θ = 10 kN/mm 2 and the Flory-Huggins parameter toχ = 2.5, cf. Fig. 2. The chemical potential results from the total derivative wrt. c, Eq. (25). A mechanically frozen state is considered, F = 1, and following Cahn-Hilliard's theory of phase decomposition, the common tangent construction is used to determine the equilibrium concentrations. Specifically, the concentrations c α and c β are calculated with a bisection algorithm for the coupled system of equations where Ψ α = Ψ c = c eq α , Ψ β = Ψ c = c eq β , μ α = μ c = c eq α and μ β = μ c = c eq β . Clearly, both μ and Ψ depend on the inelastic volume swelling J i in a nonlinear way.
We recall that we approximate J i = 1 + Ω(c − c 0 ) according to Eq. (30), i.e., there is no inelastic deformation when c = c 0 . Here, we set Ω = 0.8. The Taylor series expansion at reference point c 0 determines the reference configuration and thus affects the shape of the potentials Ψ (c). At first, we expand J i (c) at c 0 = 0 (as done in the models of [2,26,66]). The resulting double well potentials Ψ (c) are displayed in the left plot of Fig. 3. We see that a higher value ΔK = (K 0 − K 1 )/θ straightens and shifts the Ψ (c)-curve which means that a high difference of K 0 and K 1 works against phase decomposition. In other words, phase segregation can be suppressed by a large enough elastic modulus.
Also a different reference value c 0 shifts the positions of the equilibrium concentrations c eq α,β and changes the shape of the Ψ (c) curve. In the central plot of Fig. 3, we have c 0 = 0.5 and see already a convex energy Ψ for ΔK = 0.4; only small values and differences in K lead to a double well potential. The situation is even more obvious for the limit case c 0 = 1 where both curves are convex, see left plot of Fig. 3.    We conclude that concentration-dependent bulk moduli and also the reference concentration for Taylor series expansion may qualitatively change the chemical energy. To illustrate the effect over the full range of c 0 ∈ [0, 1] the equilibrium concentrations are plotted in Fig. 4 for two values of ΔK . The outer blue lines correspond to K = 0, i.e., no elastic effects on the equilibrium concentrations. There are always two phases, here with values c eq α,β = 0.5 ± 0.355. The dark red lines correspond to a high bulk modulus and indicate two equilibrium phases only for small values of c 0 ; for larger concentration, the elastic compression suppresses phase decomposition. The white stripe in the middle separates the equilibrium states. The situation is different for the strongly softening material in the right plot of Fig. 4. When K 1 K 0 , two equilibrium phases exits also under large compression.
Similar results can be obtained for an increasing relative atomic volume Ω of the intercalating ions, see Fig. 5. Here we set K 0 and K 1 as in Table 2 and observe that a high bulk modulus, rather insensitive to concentration, suppresses phase decomposition (left), whereas strongly softening material decomposes even for 'large' intercalating ions (right).
In summary, it can be said that the equilibrium concentration strongly depends on the material parameters and also on the choice of reference value c 0 . The proper choice of c 0 would be the same value employed to measure the relative atomic volume Ω. For practical application, it should be set to the initial concentration. We remark that, for other choices, numerical problems may arise, because a system with reference configuration Fig. 6 Comparison of chemical potentials μ for proposed and existing model for uniaxial stretch (left) and spherical expansion (right); setting: λ = 1.5, K 0/θ = 1, K 1/θ = 0.9 and c 0 = 0 J i = 1, initial concentration c ini and given relative atomic volume Ω is over-constrained and inconsistent for c ini = c 0 .

Comparison to other chemo-mechanical models
Anand [2] suggests a constitutive model that couples ion diffusion with large elastic and inelastic deformations of the material. This model is derived from a microforce balance theory, introduces a Mandel stress in the chemical potential, and uses phenomenological plasticity to account for inelastic effects. In the corresponding numerical implementations [26,38], several simplifications are made, e.g., the material parameters are treated as concentration independent quantities. In order to conduct comprehensive numerical simulations, Xu et al. [66] extended the model to include configurational changes, concentration-dependent elastic moduli, and a Lagrangian description of the constitutive equations. In the following, we compare Xu et al. [66] approach to our model. Anand's model [2] relies on the multiplicative decomposition of F = F e F i , where F accounts for elastic and inelastic mechanical deformations, but disregards temperature or electric effects. Macroscopic and microscopic force balances are formulated in the current configuration, and the principle of virtual power is used to derive evolution equations for the unknown quantitiesχ ,Ḟ e ,ċ, and others. For the inelastic Jacobian, this procedure givesJ i = Ωċ, supporting the conclusion that J i has the form of Eq. (29) [67], i.e., the inelastic Jacobian is set to J i (c) = 1 + Ωc (c 0 = 0 in our model). Unlike in our model, in [67] the evolution of J i is assumed to be completely compensated by the elastic swelling J e . This, in turn, follows from the underlying assumption that the right-Cauchy Green tensor C = F T F is independent of the concentration, namely C = C e (c)C i (c). As a consequence, the elastic free energy refers to the intermediate configuration. Expressed in terms of the neo-Hookean energy density (32), this leads to the form Ψ e = J i (Ψ iso + Ψ vol ). The derivation of the Piola-Kirchhoff stresses requires the application of the chain rule; thus, an additional term J i Ψ e arises from the configurational change. The remaining part of the stress is nearly the same as the first term in Eq. (52), because J i Ψ e' compares with modified elastic constants in Ψ e . Clearly, our model adopts a different kinematics. To show the differences between the models, Fig. 6 plots the chemical potential for different mechanical deformations. The curves on the left plot refer to uniaxial loading at different stretches, revealing a minor difference in the response. Our numerical experience allows to state that the results of the microforce model are very similar to the ones of our model. The observed difference in quantitative terms (i.e., magnitude of the stress) cannot be assessed because of the lack of reference data.

Electro-mechanical coupling
Using a standard finite element model, we want to assess the relevance of electro-mechanical coupling under realistic working conditions in a Li-ion anode. Typical lithiation potentials have been reported to be in the range of few mV [15,32]. Experimentally, silicon-based half cells are subjected to open-circuit voltages of 10 mV and lithiate for 24 h [60]. Exceeding a threshold potential leads to a significant reduction in cycling efficiency and to the failure of the solid-electrolyte interface [61]. We consider a Li-ion anode subject to a electric potential growing with 1 mV steps up to a maximum value of 10 mV. The specimen is cylindrical with a radius of 0.5 µm and a length of 7 µm, clamped at the two extremities, Fig. 7. The linear momentum balance (5) and the electrical induction (7) are coupled in a weak form and solved with a staggered approach [35,54]. Both bulk modulus and shear modulus are set to 10 GPa; the relative permittivity is ε r = 0.15.
The growing potential causes a growing active (inelastic) deformation, which is constrained in the longitudinal direction and free in the radial direction. The cylinder expands the radius increasing quadratically with increasing electric potential, see Fig. 7. Nevertheless, the radial displacement remains very small, with a maximum strain of ≈ 10 −6 . We conclude that the local deformations induced by the electric field are negligible, which justifies the choice to set F act to unity, Eq. 28.

Numerical examples
The multi-field problem described in the previous section has been solved by means of a finite element approach. The method requires to formulate the governing equations in weak form and to introduce a spatial discretization by means of approximation functions. In the following, we state the weak form of the balance equations for the coupled problem, and present a selection of illustrative numerical examples concerning elastic and inelastic phenomena in a cylindrical structure. The material parameters used in the simulations are summarized in Table 2.

Finite element approximation
The weak form of a differential equation requires the multiplication of the equation by a suitable test function (or variation) having the same mathematical features of the unknown field. The resulting scalar function is integrated over the physical domain, making use of the Green's formula and including natural boundary conditions.
The weak form of the mechanical problem is obtained from the material form of the linear momentum, Eq. (5), using the variation δχ The deformation χ(x, t) is dependent on concentration, temperature, and electric fields through the stress, Eq. (52). The weak form of the heat equation is derived through the variation δT as where entropy (54) and heat flux boundary ∂B H 0 are accounted for. The diffusion problem, using (11) with chemical potential (55), is described by a fourth-order partial differential Eq. (4). The corresponding weak form is where m c = m i j for all i, j; no-flux boundary conditions apply. Finally, the weak form of the electrostatic problem becomes where the electric induction is specified in (56). The unknown fields of the problems, χ, T , c, and D, are coupled through the constitutive relations. The resulting multi-physics problem is highly nonlinear and, to solve the corresponding numerical equations, it will be necessary to proceed with a suitable linearization of the constitutive relations. Linearization implies the definition of the hessian of the free energy that will include the fourth-order material tensor, the conductivity, the second-order electrical conductivity tensor, and others, as well as the corresponding coupling terms.
Because of the second-order derivatives in (61), the finite element approximation of the coupled problem requires the introduction of C 1 (at least) continuous piecewise-smooth approximation functions. Here, we use non-rational B-splines, which offer a high-order accuracy and robustness. For details of the implementation and demonstration of the efficiency, we refer to previous works [4,5]. The spline-based construction of the geometry requires to address some design issues and will be described in another paper.

Purely elastic rod
We start by analyzing the phase decomposition in an elastic cylindric specimen (rod) characterized by an initially uniform distribution of Li ions. The overall concentration is set initially to a constant value, c(X, 0) = 0.4, and then let to evolve within a time interval of 2520 s. Results are presented in terms of normalized time, i.e.,t = t/τ where τ = 2 /D, see Table 2.
The specimen has a length of l = 7 μm and a radius of r = 0.5 μm, Fig. 8. The adopted mesh is composed by 189 elements with uniform size. The specimen is clamped at both ends, i. e., longitudinal displacements are constrained on the green planes in Fig. 8. The central points of the bases, indicated with yellow dots, are constrained in all directions. A no-flux condition is imposed on all boundaries through Lagrange multipliers [4]. Figure 8 shows a sequence of snapshots of the specimen at different times. The system immediately decomposes, exclusively by diffusion, into α and β-phases, whereby the specimen deforms as a consequence of swelling due to the intercalation Li ions. With time some regions grow at the expense of other ones and the phases coarsen. The material parameters chosen for the simulation lead to a symmetric equilibrium state characterized by two phases with concentration c α = 0.855 and c β = 0.145, respectively. The application of an additional mechanical stress moves these values towards the mean concentration, as shown in Sect. 5. The expansion point for J i (c) is c 0 = 0.25. Figure 9 shows the concentration distribution and the corresponding von Mises stress distribution at the timet = 100. High-concentration regions expand and are characterized by a low stress on the cylinder surface, while low-concentration regions are characterized by a high compressive stress. The mutual dependence is less evident in the internal parts of the specimen, which reveal a more homogeneous stress.
Next, the phase decomposition problem is solved for a cylinder with unconstrained and load free bases, constrained only at the central section to exclude rigid body motion, see Fig. 10. For this free case, the interaction between diffusion and mechanical deformation is reduced to the presence of stresses related to intercalation. The specimen still shows phase coarsening and swelling, but these structures develop in much longer times. For example, at the timet = 180, the clamped specimen appears fully decomposed, whereas the free specimen is still in the initial configuration characterized by uniform concentration.

Parametric analysis
In order to assess the robustness of the proposed approach, we perform a parametric analysis considering only the geometry of the clamped specimen.  Figure 11 compares the concentration fields at the timet = 100 for three specimens of length l = 3, 6, 12 µm, respectively. The shortest specimen results decomposed entirely att = 100, whereas the concentration field in larger specimens is still evolving. The size of the evolving phases in the three specimens is similar but not equal. In particular, by measuring a single concentration phase as the distance between two half hills, the longest specimen develops five phases, whereas the medium specimen develops two phases, and the shortest develops only one phase. We can define the intrinsic wavelength of the system as approximately 2.4 µm long. Note that, because of the driving surface energy, the c α -phases arrange initially at the system boundaries for c α > c β and c = 0.5; thus, all the specimens show the presence of high concentration phases at both clamped ends. For an additional 30 µm long specimen, Fig. 12 shows an intermediate state of concentration decomposition. Figure 11 (right) compares the concentration fields obtained with two different levels of mesh refinement, revealing a good agreement in terms of concentration phase location. Interestingly, further mesh refinements do not improve the results of the fine mesh in Fig. 11 (right). Figure 13 describes the dependence of the results on the coefficient κ, a gradient energy factor that enters the interface term |∇c| 2 in Eq. (39). Using the 12 μm long clamped elastic specimen, results have been obtained for the reference value of κ reported in Table 2 and for a 10 times larger value, i. e., κ 1 = 10κ. In dimensionless form, the interface term can be written as |∇c| 2 = −2 |∇c| 2 , where∇ is a dimensionless operator and is a characteristic length. Consequently, the scaling of κ → κ 1 can be seen as a scaling of → 1 It follows the transition from αto β-phase requires a longer time and lead to the formation of a smaller number of phases. We remark that the characteristic length may enter the normalized time, but the scaling is not one-to-one because of the presence of chemo-elastic interaction. It has been shown that besides the concentration gradient coefficients, strain gradient coefficients or related gradient expressions also play a stabilizing role favoring uniform concentrations [51,55].

Elasticity models
The comparison is extended to the two elasticity models of Sect. 4.2, namely the classical neo-Hookean model with energy density (32) and the simplified finite elasticity model with energy density (33). The corresponding stresses are given by (52) and (53), respectively.
The geometry of the specimens (l = 7 μm, r = 0.5 μm, c ini = 0.4, c 0 = 0.25) and the discretization (189 elements) are similar to the ones adopted in the previous section, but the time necessary for phase decomposition is longer, since here the elastic parameters K and G are kept constant. Figure 14 compares the time evolution of the radii of the specimen as a function of the position along the longitudinal axis for the two elasticity models. At both clamped ends, the radius shows an initial jump from 0.5 to ≈ 0.53 µ, describing at the interior a waving profile characterized by three phases. Subsequently, the three

Viscoelastic rod
The effects of viscoelasticity are studied only for the 7 μm long rod. To account for the viscous behavior, we adopt a one-parameter Maxwell model with τ 1 = 100 2 /D, G 1 = G 0 , keeping the other values according to the list in Table 2. Figure 15 illustrates the time evolution of the elastic, chemical, and interface energy for the viscous model (red curves) and for the elastic model (blue curves) described in Sect. 6.2. The chemical energy plot reveals that phase decomposition is delayed in the viscous case; the minimum for the viscous line is attained at about 10 time steps later than for the elastic line. In later times, the kinetics of the reaction is accelerated as a consequence of the reduction of the elastic stresses. Figure 16 shows snapshots the deformed configurations of the viscoelastic rod.

Charging of a viscoelastic rod
We conclude with an investigation on the charging process of a 7 μm rod with initial concentration c ini = 0.25, and undergoing a constant normalized incoming charging Li flux of 0.064 at one of the basis, while the second basis is set as a no flux boundary. In order to display the results, we select ten points P i , i ∈ {1, . . . , 10}, located at the center of external finite elements, see Fig. 17. The point P 1 is close to the Li flux boundary, while the point P 10 is close to the no-flux boundary at the opposite end of the rod. The dynamics of charging is shown in the plots of Fig. 17. At the point P 1 , the concentration increases until it reaches the equilibrium concentration c α . In contrast, at the point P 10 , energy minimization leads to a local concentration c β . The highest concentration values are observed at the charging side and decrease progressively with the distance from it, i. e., c(P i ) > c(P j ), j > i. Due to the constant incoming flux, the charging process is characterized by an almost constant velocity ∂ t c(P i ) at intermediate times between the initial and final equilibrium concentrations. Plots show that, within the first 2300 time steps before full charging, concentrations remain in the interval [0, 1]. It bears emphasis that in real charging processes, the incoming Li flux is not constant but reduces with raising concentration. An analogous effect is induced by viscoelasticity.
Clearly, because of dissipative phenomena, the relevance of viscosity or inelasticity comes heavily into play when charging and discharging cycles are considered at the anode, calling for more sophisticated simulations that will be conducted in future studies.

Conclusions
We derived a multi-field model for lithium electrodes during charging and discharging. Ion flux, diffusion and phase decomposition, temperature, electricity, ion intercalation, and swelling contribute to the material internal energy and to the development of large elastic and inelastic mechanical deformations. Coupling between concentration, electricity, temperature, and displacement fields arises from the constitutive relations and is accounted for in the kinematics through the multiplicative decomposition of the deformation gradient, which splits into elastic and inelastic parts. A Taylor series expansion is used to approximate the unknown inelastic deformation due to ion intercalation. We show that the reference concentration c 0 strongly influences the decomposition behavior.
Our model is consistently derived from the basic laws of thermodynamics and reflects the fundamental aspects of the anode charging process. We have illustrated the properties of the proposed model through several examples and have pointed out that chemomechanical interaction affects the equilibrium concentrations of the phases. For typical charging processes, the electromechanical coupling remains marginal, and the model compares-for special cases-to previous electrode charging models available in the literature.
We observed that two neo-Hookean models of finite elasticity lead to similar evolutions of the concentration distribution, i. e., they deliver the same intrinsic wavelength during the coupled decomposition. A Maxwelltype viscoelastic material model reveals that a viscoelastic specimens delays the onset of phase decomposition, but the decomposition accelerates at later stages.
The finite element simulations here described demonstrate the robustness of the proposed electrochemicalthermomechanical model.