Computers and Mathematics with Applications Formulation and analysis of fully-mixed methods for stress-assisted diffusion problems ✩

This paper is devoted to the mathematical and numerical analysis of a mixed-mixed PDE system describing the stress-assisted diffusion of a solute into an elastic material. The equations of elastostatics are written in mixed form using stress, rotation and displace- ments, whereas the diffusion equation is also set in a mixed three-field form, solving for the solute concentration, for its gradient, and for the diffusive flux. This setting simplifies the treatment of the nonlinearity in the stress-assisted diffusion term. The analysis of existenceanduniquenessofweaksolutionstothecoupledproblemfollowsascombination of Schauder and Banach fixed-point theorems together with the Babuška–Brezzi and Lax– Milgram theories. Concerning numerical discretization, we propose two families of finite element methods, based on either PEERS or Arnold–Falk–Winther elements for elasticity, andaRaviart–Thomasandpiecewisepolynomialtripletapproximatingthemixeddiffusion equation. We prove the well-posedness of the discrete problems, and derive optimal error bounds using a Strang inequality. We further confirm the accuracy and performance of our methods through computational tests.


Introduction
We are interested in the mathematical and numerical study of a stationary problem representing diffusion-deformation processes where the stress acts as a coupling variable. So-called stress-assisted diffusion models (derived from thermodynamic principles and phenomenological arguments in e.g. [1,2]) are relevant to numerous applications including diffusion of boron and arsenic in silicon [3], hydrogen diffusion in metals [4], voiding of aluminium conductor lines in integrated circuits [5], strain-aging measurements in iron [6], sorption in polymers [7], to name a few. Of special appeal to us is the study of microscopic electrode damage in lithium ion batteries [8][9][10][11][12]. When lithium diffuses into a secondary particle (an anode made of e.g. silicon), its core expands and its elastic response, also with that of neighbouring particles and the surrounding electrolyte, modify the diffusive properties inside the medium. If the process is confined inside the anode, then the electric field is practically constant and the system may be described solely in terms of diffusion and stress.
Regarding the mathematical and numerical analysis of related models, the literature is rather scarce. Some recent references include homogenization of concentration-electric potential systems [13], multiscale analysis of the deterioration ✩ Funding: This work was partially supported by CONICYT-Chile through the project AFB170001 of the PIA Program: Concurso Apoyo a Centros Científicos y Tecnológicos de Excelencia con Financiamiento Basal, and the Becas-Chile Programme for foreign students; by Centro de Investigación en Ingeniería Matemática (CI 2 MA), Universidad de Concepción, Chile; and by the EPSRC, United Kingdom through the research grant EP/R00207X/1. of binder in electrodes [14], and a general local-global well-posedness theory [15]. Differently from these approaches, in [16] we have recently proposed a mixed-primal formulation for stress-assisted diffusion. The model covers the linear elastic regime, it incorporates the rotation tensor as supplementary variable serving to impose stress symmetry in a weak manner; and this mixed problem is coupled with a primal formulation for diffusion. Here, in contrast, we consider an augmented mixed formulation for the diffusion equation. Similarly to [17], the concentration gradient and the diffusive flux are incorporated as auxiliary unknowns, which allows us to treat the stress-dependent diffusivity using a dual-mixed setting. In order to apply the regularity estimates from [16], we augment the formulation with redundant terms arising from a constitutive equation. Next, following the approach introduced in [18], we combine fixed-point arguments, regularity estimates, the Babuška-Brezzi theory, the Lax-Milgram lemma, the Sobolev embedding and Rellich-Kondrachov theorems, and small data assumptions to establish existence and uniqueness of solution of the continuous problem. The solvability of the Galerkin scheme follows from the Brouwer fixed-point theorem and properties of the finite element subspaces. Finally, the convergence analysis is conducted adapting Strang inequalities, Céa estimates, and using approximation properties of the finite element spaces.
The rest of the paper is organized as follows. In Section 2 we describe required notation and functional spaces to be employed along the paper. Then, we introduce the model problem and requirements on the specific constitutive functions. Next, in Section 3 we derive the augmented fully-mixed formulation and establish its well-posedness. The Galerkin scheme and the existence of discrete solution are then studied in Section 4. In addition, under similar assumptions we deduce error bounds in Section 5; and we close in Section 6 with a numerical example that confirms the theoretical rates of convergence, and a second test studying the applicability of the discrete formulation in the simulation of 3D microscopic lithiation processes.

The model problem
Let Ω ⊆ R n , n ∈ {2, 3}, be a given bounded domain with polyhedral boundary Γ = ∂Ω and denote by ν the outward unit normal vector on the whole boundary ∂Ω. Standard notation will be adopted for Lebesgue spaces L p (Ω) and Sobolev spaces H s (Ω) with norm ∥·∥ s,Ω and seminorm |·| s,Ω . In particular, H 1/2 (Γ ) is the space of traces of functions of H 1 (Ω), while H −1/2 (Γ ) denotes its dual. By |·| we will denote both the Euclidean norm in R n and the Frobenius norm in R n×n . Let div τ be the divergence operator acting along the rows of the tensor τ. We recall that the tensorial H(div) space equipped with the norm ∥τ∥ 2 div;Ω := ∥τ∥ 2 0,Ω + ∥div (τ)∥ 2 0,Ω ∀ τ ∈ H(div; Ω), is a Hilbert space. Let I stand for the identity tensor in R n×n . For any tensors τ = (τ ij ) i,j=1,n , and ζ = (ζ ij ) i,j=1,n , we denote the transpose, trace, tensor product, and deviatoric tensor, respectively, as Finally, we will denote by ∥·∥ ∞,Ω the norm of the Banach space L ∞ (Ω) as well as of its vectorial version L ∞ (Ω).
Let us consider the following system of PDEs, governing the diffusion of a solute interacting with the motion of an elastic solid occupying the domain Ω: The involved quantities and model parameters are the Cauchy solid stress σ, the displacement field u, the infinitesimal strain tensor ε(u) := 1 2 ( ∇u + ∇u t ) , the Lamé constants λ, µ > 0 characterizing the material, the diffusive fluxσ, the solute concentration φ, the tensorial diffusivity ϑ : R n×n → R n×n , the vector of body loads f : R → R n , and a displacementdependent source term g : R n → R. For the load, source, and diffusivity functions we will require uniform boundedness and Lipschitz continuity, that is there exist positive constants f 1 , f 2 , L f , g 1 , g 2 , L g , and ϑ 1 , ϑ 2 , L ϑ , such that Additionally, ϑ is of class C 1 and uniformly positive definite, the latter meaning that there exists ϑ 0 > 0 such that Finally, we assume that f (φ) ∈ H 1 (Ω) for each φ ∈ H 1 (Ω), and that for each γ ∈ (0, 1) there exists a constant C γ > 0 such that g(w) ∈ H γ (Ω) for each w ∈ H γ (Ω) and ∥g(w)∥ γ ,Ω ≤ C γ ∥w∥ γ ,Ω . (2.8) Examples of stress-dependent diffusivity functions and concentration-dependent body loads may include exponential functions of the volumetric stress for lithiation of batteries [10], simple polynomial relationships for biological materials [19], or Carreau-type laws for ϑ, that is respectively, where C 0 , C 1 , C 2 are constants, whereas for f linear dependences modelling isotropic swelling in composite materials [20], saturation-based descriptions for viscous layers [21], or concentration gradient modulations for single-cell mechanics [22] are considered, that is respectively, where C ∈ R n and m > 1. Nevertheless, not all of these fulfil each one of the above described hypotheses. For instance, the constitutive equation for diffusivity ϑ(σ) = C 0 I + C 1 σ + C 2 σ 2 violates (2.6) and it is not necessarily satisfying (2.7) for certain stress configurations. Illustrations of the latter issue are discussed in [19].

Weak formulation and solvability analysis
In this section we derive an augmented fully-mixed variational formulation for (2.1)-(2.3) and propose a fixed-point strategy for its analysis. We show that the fixed-point operator is well-defined and apply the Schauder's theorem to prove existence of solution, whereas Banach fixed-point theorem will lead to uniqueness of solution under small data assumptions.
In turn, defining the concentration gradient t := ∇φ, we can recast the diffusion equation as We then test the three-field problem (3.3) against s ∈ L 2 (Ω),τ ∈ H(div; Ω) and ψ ∈ L 2 (Ω). Integrating by parts the expression ∫ Ω ∇φ ·τ and using the Dirichlet boundary condition for φ (second equation in (2.3)), we arrive at the weak formulation: find (t,σ, φ) ∈ L 2 (Ω) × H(div; Ω) × L 2 (Ω) such that ∫ In view of modifying the regularity properties of the coupled problem, we proceed to enrich the foregoing equations with the following residual terms: where κ 1 , κ 2 , κ 3 and κ 4 are positive parameters to be specified later on. We remark that the identities required in (3.5) are nothing but the constitutive and the equilibrium equations concerningσ, along with the relation defining t, and the Dirichlet boundary condition for φ; all of them tested differently from (3.4). Instead of (3.4), we will now focus on the following augmented formulation for the diffusion problem: and Consequently, we arrive at the following augmented fully-mixed formulation for (2.1)-(2.3): find (3.9)
We also collect the following two technical lemmas, whose proofs can be found in [ In what follows we show that T has at least one fixed point. Firstly we will prove that the uncoupled problems defined by S andS are well-posed, where we emphasize that S is defined similarly as in [16], and therefore we omit parts of the proofs whenever necessary. Our analysis will focus on the uncoupled problem (3.6) and its repercussion on T. Let us start by recalling the continuity of a and b. For a proof we refer to e.g. [24].
Furthermore, it is not difficult to see that a is strongly elliptic in the kernel of b. In fact, we denote the operator induced by the bilinear form b as B, and note that where c 1 is the constant provided by Lemma 3.1. Additionally, as a slight modification of the proof of [24, Section 2.4.3], we find that B is surjective. Finally, we observe that G and F φ are bounded with This analysis confirms the well-posedness of (3.1), which is abridged in the following lemma. (3.14) Proof. It follows from estimates (3.11)-(3.13) and a direct application of the Babuška-Brezzi theory (see, e.g. [ In turn, we prove the well-posedness of problem (3.6) with the next result.

Solvability analysis of the fixed-point equation
We now verify the hypotheses of the Schauder fixed-point theorem (see, e.g. [31, Theorem 9.12-1]). Before starting the result to be proved, we restrict T to a ball and show that this operator maps into itself.

Lemma 3.5. Let W be the closed and convex subset of H 1 (Ω) defined by
wherec S is the constant given by (3.14). Then T(W ) ⊆ W .
We are in a position to establish the announced property of the operator T.
Analogously to the proof of [16, Lemma 2.7], we apply the ellipticity of A σ (cf. (3.18)) and then, by adding and subtracting appropriate terms, we find that . (3.26) Using the definition of A σ , G u , Cauchy-Schwarz's inequality, and (2.5), (2.6), we can assert that (3.29) Next, according to the definitions given when starting the proof, we can rewrite (3.29) as It is important to note here that, when needed, ∥S 1 (S 1 (ϕ), S 2 (ϕ))∥ ∞,Ω can be bounded by (3.23), for each ϕ ∈ H 1 (Ω). Finally, applying estimates (3.24) and (3.30), we find that which proves continuity of T. In turn, let {φ k } k∈N be a sequence of W , which is clearly bounded. Then, there exists a subsequence {φ (1) k } k∈N ⊆ {φ k } k∈N and φ ∈ H 1 (Ω) such that φ (1) k w − → φ ∈ H 1 (Ω). In this way, thanks to the compactness of i c , we deduce that φ (1) k → φ ∈ L 2 (Ω), which, combined with (3.25), implies that T(φ (1) k ) → T(φ) ∈ H 1 (Ω), and proves the compactness of T(W ). □ We are ready now to prove that (3.10) is well-posed. From Lemmas 3.5 and 3.8, the existence of solution is merely an application of Schauder's theorem. Furthermore, assuming that the data is small enough, we can prove uniqueness of solution. This is indeed possible thanks to the regularity estimates established at the end of Section 3.2.
Details of the proof are similar to those available in [16, Thm. 2.9].

Theorem 3.9.
Let W be as in Lemma 3.5. Then, the augmented fully-mixed problem (3.9) has at least one solution Moreover, if the data satisfy then the solution φ is unique in W.

The Galerkin scheme and well-posedness of the discrete problem
In this section we introduce and analyse a Galerkin scheme for (3.9). We adopt the discrete analogue of the fixed-point strategy introduced in Section 3.2 and apply the Brouwer fixed-point theorem to prove existence of discrete solution. We start by considering generic finite dimensional subspaces which will be specified later on. Hereafter, h denotes the size of a regular partition T h of Ω into triangles K of diameter h K , i.e. h := max {h K : K ∈ T h } . A Galerkin scheme for (3.9) reads: In order to address the well-posedness of (4.3), we proceed analogously as in Section 3.2 and apply a fixed-point strategy.
In fact, we define the operator with A σ h and G u h being defined by (3.7) with σ = σ h and (3.8) with u = u h , respectively. In this way, by introducing the operator 3) can be rewritten as the fixed-point problem: Analogously to the continuous case, we first study the well-posedness of S h andS h , and hence the well-definiteness of T h .
To this end we proceed as in [ and assume the following discrete inf-sup conditions: (H.0) There exists a constant α 1 > 0, independent of h, such that (H.1) There exists a constant β 1 > 0, independent of h, such that Deriving well-posedness of the discrete problem (4.4) results as a straightforward application of the discrete Babuška-Brezzi theory. Firstly, the operators related to a and b, and the functionals G and F φ h are all bounded on subspaces of the corresponding continuous spaces. Next, the inf-sup conditions are given by (H.0) and (H.1). The unique solvability of (4.4) is abridged in the following lemma.

Lemma 4.1. For each
In regard to the problem defined byS h we state next the discrete analogue of Lemma 3.4.
,δ ∈ (0, 2), and κ 2 , κ 4 > 0. Then, Proof. We first observe that for each We notice in advance that, instead of the regularity estimates employed in the continuous case (not applicable in the present discrete case), we simply utilize properties of the discrete subspaces chosen. In what follows, we verify the hypotheses of the Brouwer fixed-point theorem (see, e.g. [31,) to prove that T h has at least one fixed point.

Lemma 4.5. For each
(4.10) Proof. Proceeding as in [16,Lemma 3.5], Then, analogously to the proof of Lemma 3.7, we get At this point, we are able to state the main result of this section.

Then (4.3) has at least one solution
Proof. After using Lemmas 4.3 and 4.6, the result is a straightforward consequence of Brouwer's fixed-point theorem. In turn, bounds (4.11) and (4.12) follow from Lemmas 4.1 and 4.2, respectively. □

Error analysis for the proposed Galerkin method
In this section we advocate the derivation of error estimates for (4.3). For this purpose, we consider in what follows be the solutions of (3.9) and (4.3), respectively. We seek an upper bound for for which, we suggest to estimate Next, we recall from [32, Thm. 11.2 and 11.1] two instrumental results. First, a Strang inequality for saddle point problems where continuous and discrete formulations differ only in the functional. This will be applied to (5.1). Second, the standard Strang Lemma for elliptic problems, which fits (5.2). We will not write them explicitly here but will refer to these lemmas as Saddle-point Strang Lemma, and Elliptic Strang Lemma, respectively. From now on, we denote as usual Next, a straightforward application of the Saddle-point Strang Lemma yields the following result concerning a priori estimates for ) In turn, an estimate for ∥(t,σ, φ) − (t h ,σ h , φ h )∥ reads as follows.
Proof. Applying the elliptic Strang Lemma in the context of (5.2), gives Then, proceeding analogously as in the proof of Lemma 3.7, we deduce that In turn, in much the same way as [33, Lemma 5.2], we add and subtract suitable terms to write thus, the estimates for the first and third terms follow by applying the boundedness (3.16), whereas for the second one, we find that whence, we deduce that Finally, by replacing (5.6)-(5.7) into (5.5), we get (5.4), which ends the proof. □ Now, to derive the Céa estimate for the total error we combine Lemmas 5.1 and 5.2. To this end, and for notational convenience, we introduce the following constants Next we replace the bounds for ∥u − u h ∥ 0,Ω and ∥σ − σ h ∥ 0,Ω into (5.4), and apply from (3.23) that We then perform algebraic manipulations to find that ) Consequently, we can establish the following result which provides the complete Céa estimate.

Theorem 5.3. Suppose that the data satisfy
Then, there exist positive constants C 4 and C 5 independent of h, such that . integer k ≥ 0, for each K ∈ T h we let P k (K ) be the space of polynomial functions on K of degree ≤ k and define the local Raviart-Thomas space of order k as 2 , and x is the generic vector in R 2 . Let b K be the element bubble function defined as the unique polynomial in P k+1 (K ) vanishing on ∂K with ∫ K b K = 1. Then, for each K ∈ T h we consider the bubble space of order k, defined by

Numerical results
In this section we present some examples illustrating the performance of our augmented fully-mixed scheme (4.3), and confirming the rates of convergence provided by Theorem 5.4. These numerical results also include examples in which some of the data do not necessarily satisfy all the hypotheses required, thus confirming the potentiality of the method proposed, and also evidencing that only technical limitations are preventing us from extending our theoretical analysis to more general cases. Our implementation is based on the FEniCS library [36]. In turn, a Picard algorithm with tolerance of 10 −6 on the ℓ ∞ -norm of the residual has been employed for the fixed-point problem (4.6). The boundary conditions employed in Examples 2 and 3 were motivated by the specific application of stress-assisted diffusion problems in lithium batteries, and they correspond to mixed boundary conditions, which are currently not supported by our theoretical analysis. Nevertheless the obtained results still show stable and robust computations, which insinuates that only technical difficulties prevent us of extending our analysis to the case of mixed boundary conditions. For the diffusion sub-problem in Example 2 we have utilized the variational formulation (3.4) applying a fixed-point on σ and t; and σ and φ, respectively. Example 1. In our first numerical test we take the unit square as computational domain Ω = (0, 1) 2 and choose the following manufactured exact solutions and coupling terms to (3.9): , g(u) = 2 + 1 1 + |u| 2 . 2(D 0 + D 1 ). According to (3.20), the stabilization parameters are taken as κ 1 = ϑ 0 /ϑ 2 2 , κ 2 = ϑ 0 /2 ϑ 2 2 , κ 3 = ϑ 0 /2 and κ 4 = ϑ 0 /4. The convergence of the approximate solutions is assessed by computing errors in the respective norms and experimental rates, that we define as usual where e,ê denote errors computed on two consecutive meshes of sizes h,ĥ, respectively. We choose the finite element spaces (5.11) and (5.12), that is BDM k+1 − P k − P k − RT k − P k − P k+1 approximations with k = 0 and k = 1. Errors and decay rates are summarized in Table 6.1, where we observe that optimal convergence O(h k+1 ) is attained for all fields in their relevant norms. These findings are in agreement with the bounds given by Theorem 5.4. In all cases, three Picard steps were required to reach the desired tolerance. Sample solutions are displayed in Fig. 6.1.

Example 2.
Next we concentrate on the simulation of microscopic lithiation of an anode. Details on model derivation and physical considerations can be found, for instance, in [9,12,13], whereas the specific settings that motivate this example are summarized in [38]. The domain consists of a truncated sphere of radius 10 µm, representing the silicon core of a secondary particle (see Fig. 6.3(a)), which we discretize using an unstructured mesh of 104913 tetrahedral elements. For this test we consider lowest order Raviart-Thomas elements for the flux and the concentration gradient, and piecewise polynomials for the concentration (see the method developed in [39]). We assume that the face of the truncated sphere which is closest to the plane x 1 = 0 (denoted Γ D ) is in contact with a region of electrolyte, that is, the zone between the sphere and the surrounding cube. On Γ D we set zero-flux of lithium and also consider that the anode has an external layer that does not permit displacement of the body, so there we set u = 0. On the remainder of the boundary, Γ N = Γ \ Γ D , we prescribe a maximum lithium concentration φ = φ max with φ max = 26390 mol m −3 , as well as σν = βφIν, where β is a parameter to be specified later on. We assume that the source term is zero, and the diffusivity is specified as ϑ(σ) = D 0 I + D 1 σ with D 0 = 1.2 e−21 m 2 s −1 , D 1 = 3.9e−14 m 2 s −1 , and the elastic material properties of silicon are E = 60 GPa and ν = 0.25. Following the referenced models, here the total stress contains a contribution due to lithium concentration. More specifically, we consider σ tot = σ − βφI, with β =Ω(3λ + 2µ)/3, whereΩ = 4.926e−6 m 3 mol −1 is the partial molar volume. The balance of momentum is then −div σ = −β∇φ, or equivalently −div σ tot = 0 and the zero traction boundary condition can be recast as σ tot ν = 0 on Γ N . Example 1: Convergence history and Picard iteration count for the augmented BDM k+1 − P k − P k − RT k − P k − P k+1 approximations with k = 0, 1. Here N stands for the number of degrees of freedom associated to the each triangulation T h . Augmented In order to have a model with fewer chemical and physical parameters, and also to accommodate a model with adimensional units, we proceed to rescale the strong form of the governing equations and testing different deformation regimes to match the expected values found in the literature. We introduce the following parameters: the intrinsic size of the domain L = 1.6e−5 m 2 , ∇ * = ∇/L, div * = div/L, φ * = φ/φ max , u * = u/L and σ * = L 2 σ. Thus, taking D 0 = 1.0e−2D 1 L 2 , we reduce the parameters D 0 , D 1 , β,Ω given above to only β * = βφ max /L 2 and α = 1.0e−2D 1 L 2 φ max . Making abuse of notation, we rename β * by β, u * by u, and so on. The proper scaling of the parameters implies that the baseline case corresponds to β = 5.0e1 and α = 1.0e−3. Fig. 6.3(i) illustrates the sharp transition between high and low concentrations as lithium diffuses from Γ N into the secondary particle. In addition, Fig. 6.3(c) shows more pronounced displacements near Γ N (which is precisely the region where the silicon is fully lithiated), and the particle swelling is indeed influenced by the lithium gradient distribution. The stress-assisted diffusion mechanism together with the dilation-dependent source term, also contributes to maintain maximum lithium concentration near Γ N . This two-way coupling effect implies in turn that the lithium concentration is less important in regions where the secondary particle is clamped.
In Fig. 6.2 we show three different constitutive relations defining ϑ as a function of the first component of the Cauchy stress tensor. The first and third specifications correspond to the functions used in this test and in the accuracy example, respectively, whereas the second relationship has been used in [19] in the context of biological materials. Depending on the values attained by the stress, one could then easily derive the values of the augmentation constants. On the other hand, in Figs. 6.2(b) and 6.2(c) we report a study on the influence of different values of the coupling constants β and α into the norms of selected solution fields for the elasticity and diffusion problems. We remark that the ℓ ∞ -norm of u h is practically invariant to moderate values of β, but it increases abruptly when this parameter approaches 70. Furthermore, the L 2 -norm of the stress increases linearly with β. As this constant drives the intensity of the deformation as well as the coupling strength, we also observe an increase in the Picard iteration count (where we stress that all fields are normalized). We also observe an increase of the L 2 -norm of the concentration gradient with respect to α, while for smaller values of α the method produces higher values of the L 2 -norm ofσ. Example 3. In our last example we test a similar model defined on a perforated cylindrical particle (see a sketch in Fig. 6.4(a)).
The problem setup has been adapted from [13]. The outer and inner radii of the bases are 5 µm and 1 µm, respectively, and the height of the cylinder is 25 µm. We discretize the domain using an unstructured mesh of 101907 tetrahedral  elements, and employ the method that uses the lowest-order spaces defined in (5.10), (5.12). In this test we consider that the particle is clamped on the inner wall Γ I , while zero lithium fluxes are prescribed on Γ B ∪ Γ I . Also, we fix a maximum lithium concentration on Γ O , whereas zero traction will be imposed on Γ B ∪ Γ O . We let E = 10 GPa, ν = 0.3 andΩ = 3.497e−6 m 3 mol −1 . The diffusive source is zero and the diffusivity tensor and body load source are given by ϑ(σ) = αI+α 2 σ+α 3 σ 2 and f (φ) = βrφ, respectively, where r is the radial vector r = (x, y, 0) t and α, β are the adimensional parameters given in Example 2, assuming the values α = 5.0e−3 and β = 75. Fig. 6.4 shows the approximate solutions, indicating that the cylindrical particle deforms on the faces and outer radius and having a more important displacement on the faces. Finally, as in Example 2, we observe that the lithium concentration induces the swelling of the cylindrical particle, however as on the faces Γ B we now have zero-traction and zero concentration flux conditions coexisting, the lithium concentration is no longer maximal on the outer radius.