Enriched continuum for multi-scale transient diffusion coupled to mechanics

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website. • The final author version and the galley proof are versions of the publication after peer review. • The final published version features the final layout of the paper including the volume, issue and page numbers.


Introduction
Coupled diffusion-mechanics problems arise in many application areas, when the diffusion of solute particles causes volumetric swelling of a host material, inducing chemical stresses which in turn affect the mass flux [1,2]. It has a broad range of applications ranging from biological tissues to microelectromechanical systems devices. For example, the swelling of brain tissues, known as edema, due to water diffusion [3] or the bending of thin plates due to chemical saturation [4]. Another typical example is the swelling of the active material due to the lithiation process in lithium-ion batteries [5], which is also governed by coupled diffusion-mechanics phenomena.
The fundamental physics behind coupled diffusion-mechanics takes place at the atomic scale where the atomic or ionic diffusion occurs [1,6]. The diffusion rate of the solute particles and the swelling of the host material depends on the atomic size of the materials involved [7] and on the activation energy which causes the jump of atoms inside the crystal lattice [8]. The jump directions and the frequencies are affected by the stresses inside the material, which in turn alter the activation energy and hence the mass flux. At the continuum level, the diffusion of species are described as driven by the gradient of the chemical potential [9]. The induced chemical stresses affect the chemical potential, which in turn influences the mass flux in the material [9][10][11][12]; this is known as the Gorsky effect [13].
This article addresses the application of coupled diffusion-mechanics described by the simulation of swelling in lithium ion batteries. A lithium ion battery consists of four components: two electrodes-a cathode and an anode, an electrolyte and the separator. Through an electro-chemical reaction, the chemical energy is converted to electrical energy in a discharge cycle; the reverse reaction takes place during the charging cycle. During charging the chemical potential across the cell forces the lithium ions to diffuse towards the anode compartment via the electrolyte while passing through the separator [14]. At the anode, the lithium ions are deposited in the active particles during an intercalation process which increases the volume of the active particle. Upon discharging, a similar reaction occurs in which the lithium in the anode is oxidized into lithium ions and electrons. The electrons flow through the external circuit to the cathode and lithium ions diffuse towards the cathode where they intercalate into the active particles.
The amount of swelling of the active particles depends on the cathode and anode materials. For example, swelling in the cathode of up to 6.5% is reported in different lithium-metal-oxides and up to 10% in lithium-cobalt-oxide [15], Silicon-based anode active particles can swell up to 300% [16]. Even when the deformation of cathode materials is small, e.g. LiCoO 2 , LiMn 2 O 3 and LiFePO 4 , the cyclic lithiation and delithiation of active particles leads to cracks and loss of contact with the matrix, which gradually results in a capacity loss and eventually failure of the battery [17,18]. Hence, to design a longer lifetime and higher energy density batteries, simulation of coupled diffusion-mechanics is of primary importance [19,20].
Most of the work done in the literature on the simulation of coupled diffusionmechanics in batteries is based on the pioneering works of Larché and Cahn [21], in which a framework for solid-state diffusion was developed for compositional changes in the solid state [5]. In general, due to its multiphysics and multiscale nature, the simulation of lithium ion batteries is a challenging task [17]. Analytical methods for the solution of coupled diffusion-mechanics problems are limited to simple geometrical shapes [22], therefore approximate solutions using numerical techniques such as finite elements are often required [17]. However, with a complex microstructures [23] and transient phenomena [24], the direct numerical simulations (DNS) become prohibitively expensive.
Computational homogenization is a well known technique to reduce the computational costs associated with the modeling of physical phenomena in complex microstructures [25,26]. It replaces a highly heterogeneous medium with an equivalent homogeneous one by decomposing the problem into smooth macroscale and highly oscillatory microscale problems. The effective behavior is computed from a representative microscopic element (RVE) [27] and transferred to the macroscale. The computational homogenization of transient phenomena, as associated with lithium diffusion in batteries, has been the focus of research recently [24,28]. Effective responses have to be computed at each macroscopic material point at each time step, making homogenization of transient phenomena computationally demanding. For a general overview of multi-scale computational modeling of lithium ion batteries, see [29,30].
In this work, we propose a computationally efficient method for the homogenization of coupled diffusion-mechanics for the cathode material of a lithium ion battery. The homogenization of the underlying diffusion and mechanical problems is performed separately by using the method proposed in [24]. For the diffusion problem, equivalence of virtual power (extended Hill-Mandel condition) is considered, while for the mechanical problem equivalence of virtual work is used (standard Hill-Mandel). Assuming linear material properties and small strains, the relaxed separation of scales allows the decomposition of the microscopic fields into their steady-state and transient parts. The mechanical response relies on the assumption of full scale separation since the characteristic time of the elastic deformation for the considered problem is very small compared to the characteristic diffusion time [2,9]. Moreover, the characteristic diffusion time in the active material particles is several orders of magnitude larger, than the one of electrolyte (considered here as a matrix in which active particles are embedded) [28]. Therefore, the lithium ions travel instantly in the electrolyte as compared to their diffusion in the active material. This allows for a so-called relaxed separation of scales, in which diffusive species migrate instantly in the matrix and very slowly in the inclusions [31]. Next, a model reduction technique, inpired by [32] for elasto-dynamic problems and applied in [31] to heat conduction problems, can be performed to extract the reduced bases for the steady-state and transient parts of the microscopic response. Although mechanical inertia effects can be neglected, the mechanical deformation is coupled to the transient diffusion, i.e. it evolves in time with the concentration field. Hence, a decomposition of the microscopic displacement field into a steady-state and a transient part is also required.
Through model reduction, the fine scale coupled-diffusion equations are replaced by a set of ordinary differential equations for the emergent macroscopic field variables, giving rise to an enriched continuum at the macroscale. These equations are to be solved at the macroscale together with the macroscopic mass and linear momentum conservation and the constitutive effective mass flux, rate of change of concentration and stress, obtained through the reduced order homogenization. The resulting enriched continuum macroscopic problem is computationally significantly less expensive than the original fully resolved problem or the direct transient computational homogenization.

Outline
The general framework of the coupled diffusion-mechanics framework is presented in "Coupled diffusion-mechanics formulation" section, where the classical formulation in terms of concentration and strain is summarized. Next, a computationally more convenient formulation expressed in the terms of the chemical potential and strain is derived. "Computational homogenization" section presents the homogenization framework, in which the relaxed separation of scales is defined. The downscaling is performed and the macroscopic effective constitutive responses are obtained through an upscaling procedure. The model reduction is carried out in "Model reduction leading to an enriched continuum" section. First, a finite element discretization is introduced and the partitioned equations are shown. The reduced bases are identified, the macroscopic quantities are written in terms of the coefficients of the reduced bases and finally mode selection criteria are discussed. Numerical examples for the cathode material of a lithium ion battery are presented in "Numerical examples" section.

Symbols and notation
Macroscopic quantities are represented with a bar on top: for example scalar, vector and second-order tensor macroscopic quantities are written as a, a, and A, respectively. Microscopic quantities are represented without a bar; microscopic scalar, vector and second-order tensorial quantities are written as a, a and A, respectively. The same Cartesian basis is adopted at the macro and micro scales. The dot products between two vectors, and between a second-order tensor and a vector are represented as a · b := a i b i and A · a := A ij a j e i , respectively. A tensorial dyadic product is denoted as a ⊗ b := a i b j e i ⊗ e j and A ⊗ a := A ij a k e i ⊗ e j ⊗ e k . The gradient of a scalar and a vector is defined as ∇a := ∂a ∂x i e i and ∇a := ∂a i ∂x j e i ⊗ e j . Similarly, the divergence operates as ∇ · a := ∂a i For linear algebra operations, columns are represented with a tilde underneath a lowercase letter, e.g. a , and matrices are represented with a bar underneath an uppercase letter e.g. A . The matrices and columns of vectors and tensor quantities are written with bold symbol, for example a matrix of a vector or a tensor quantity is written as A . A tensorial product between two column arrays of vectors is defined as a T ⊗ b , where The microscopic domain and its boundary are represented by and ∂ , respectively. The volume average of a microscopic quantity • is defined as where V = d is the volume of the microscopic domain .

Coupled diffusion-mechanics formulation
Coupled diffusion-mechanics equations describing the fully resolved (heterogeneous) problem are presented in this section. The conservation laws and the boundary conditions are written for the chemical and mechanical problems, followed by the derivation of the form of the constitutive equations [9]. First, the formulation considering the concentration and the displacement (strain) as the primary field variables is presented, which requires C 1continuity and is therefore cumbersome to implement numerically. Next, using a Legendre transform, the primal field variables are transformed to the chemical potential and strain [24]. This formulation requires only C 0 -continuity and standard finite elements can be used for the implementation. Finally, the material model to be used for the microscale constituents is presented.

Conservation laws
To take into account the large time scales associated with the mass diffusion problem a transient mass conservation equation is considered (without the volumetric source/sink term) which states that the divergence of the mass flux j in a domain 1 is opposite to the time rate of the concentration fieldċ. Equation (3) is supplemented with Dirichlet and Neumann boundary conditions, plus an initial condition where c is the prescribed value of the concentration field on the Dirichlet part of the boundary ∂ c , and j n is the prescribed normal outward mass flux on the Neumann part of the boundary ∂ j n such that ∂ c ∪ ∂ j n = ∂ and ∂ c ∩ ∂ j n = ∅. The initial value of the concentration at time t = 0 is denoted by c 0 . Considering the short characteristic times of phenomena associated with the mechanical problem, it is justified to assume a conservation of linear momentum neglecting inertial terms, which without volumetric forces reads requiring the divergence of stress field σ in a body vanish. Conservation of linear momentum (5) is also supplemented with the Dirichlet and Neumann boundary condi- where u is the displacement field, u is the prescribed displacement value on the Dirichlet part of the boundary ∂ u , and t n is the traction force applied on the Neumann part of the boundary ∂ t n such that ∂ u ∪ ∂ t n = ∂ and ∂ u ∩ ∂ t n = ∅. Constitutive equations for the mass flux j, the concentration c and the stress σ are required to close the problem (3)- (6).

(c, ε) formulation
Following [9], the dissipation inequality for a coupled diffusion-mechanics problem can be written as where μ is the chemical potential, ϕ is the dissipation density at a material point x anḋ ψ is the time derivative of the Helmholtz's free energy density. For coupled diffusionmechanics problems, the Helmholtz's free energy density ψ depends on the concentration field c and the strain ε, related to the displacement field u by ε = sym(∇u) (assuming linear kinematics). Using the chain-rule, its material time derivative can be written aṡ substituting the expression ofψ from Eq. (8) into the dissipation inequality (7) and rearranging terms yields In the inequality (9), the restriction on the dissipation density to be positive is partially fulfilled by setting which provides, for a given expression for the Helmholtz potential ψ, the constitutive equations for the stress and the chemical potential, respectively. Considering a quadratic Helmholtz free energy density [9] results in linear constitutive expressions for the stress and the chemical potential given by and where C is the elastic stiffness tensor, S the chemical strain modulus tensor, c 0 the initial concentration and is the chemical modulus. The constitutive model based on energy density function given in Eq. (11) is an alternative approach to the approach in which the microscopic strain field is decomposed in an elastic and volumetric swelling part, for more details see [33] . The remaining dissipation term in (9) − j · ∇μ ≥ 0, asserts a restriction on the constitutive form of the mass flux j. Here, we use Fick's second law which states that the mass flux j depends linearly on the gradient of the chemical potential ∇μ i.e.
where M is the second-order mobility tensor which has to be positive definite to satisfy (14). Next, the constitutive Eqs. (12), (13) and (15) can be introduced in the mass conservation Eq. (3), and in the conservation of linear momentum (5) The mass and the linear momentum conservation Eqs. (16) and (17) can be solved together for the concentration and displacement fields (c, u). Equation (16), however, involves the third-order derivative of u and its numerical solution therefore requires a C 1 -continuous finite element formulation. Various other solution techniques have also been proposed in the literature for this type of problems, see for example [24,34]. In the current work, following [24], a Legendre transform is performed on the Helmholtz's free energy density function ψ(c, ε) to obtain a dual energy density function ω, for which the primary field variables are the chemical potential μ and the strain ε.

(μ, ε) formulation
Now, we derive the constitutive equations for stress σ, concentration c and mass flux j considering (μ, ε) as the primary field variables. A Legendre transform can be performed on the Helmholtz's free energy density function (11) to obtain the dual energy density function ω which is now a function of the chemical potential μ and the strain ε. The constitutive equations for the concentration and the stress fields (c, σ) can be obtained by the standard Coleman-Noll procedure. Substituting ψ = ω + μc from (18) into (7) provides the dissipation inequality Using the chain-rule, the time derivative of the dual energy densityω can be written aṡ Substituting the expression forω, from Eq. (20), into the dissipation inequality (19) and rearranging terms yields From here, the constitutive forms for the stress σ and the concentration field c are found as Using constitutive Eq. (22) in conjunction with (18) and (11) provides the constitutive equations for the stress and for the concentration field For the remaining dissipation term −j · ∇μ ≥ 0, again Fick's second law (15) can be used. Introducing the stress (23), the concentration field (24) and the mass flux (15) into the mass conservation (3) gives while the conservation of linear momentum (5) reads Equations (25) and (26) are solved for the chemical potential μ and the displacement u.
The requirement of C 1 -continuity on u is now relaxed by using the (μ, ε) formulation, as can be seen from (26), for which a standard C 0 -continuous finite element formulation can be used.  Fig. 1 Steps involved in the development of the enriched continuum formulation: first, performing computational homogenization, the solution of a highly heterogeneous problem also known as the direct numerical simulation (a) is replaced by a two-scale macro and microscale coupled problem (b), next model reduction is performed at the microscale yielding an enriched continuum formulation (c)

Linear isotropic constitutive model
A isotropic material model is considered for both mass diffusion and mechanical problems. The isotropic mobility tensor is given by where I is the second order identity tensor and M is the scalar mobility coefficient. The chemical strain modulus S is assumed to have the following form [2] where K = 3λ+2G 3 is the bulk modulus, λ, G are Lamé's constants and γ is the partial molar volume of the material, which is the volumetric increase of a material by the introduction of one mole of other substance. The linear elastic stiffness tensor C is expressed in terms of Lamé's constants as where I is the fourth order identity tensor. Next, the computational homogenization framework for the two-scale coupled diffusion-mechanics problem will be presented.

Computational homogenization
In this section, the computational homogenization of a two-scale coupled diffusionmechanics problem is presented. First, the separation of scales regimes are defined for the mass diffusion and mechanics problems. Then, the governing equations at the micro-and the macro-scales are presented. The boundary conditions on the microscopic domain are defined through the constraints imposed by the downscaling relations. Finally, the upscaling is performed via equivalence of the virtual powers of the macro-and microscales providing the constitutive forms for the macroscopic quantities. The solution of the coupled diffusion mechanics problem on the fully resolved heterogeneous domain, as shown in Fig. 1a, is referred to as direct numerical simulation (DNS). Due to the computational expense of the DNS problem it is preferred, when possible, to divide the problem into micro and macro scales and solve a homogenized problem in a two-scale manner, as shown in Fig. 1b. The homogenizability of the DNS problem depends on the separation of scales, which is discussed next.

Separation of scales
The separation of scales can be defined through the material properties of the constituents, their characteristic length and time scales, and the characteristic scales of the physical phenomena under consideration [35]. For the coupled problem studied here, coupled scales for the mass diffusion, have to be considered, while for the mechanical phenomena a full separation of scales can be assumed. For more details on separation of scales see for example [35].
Mechanics: For the mechanical problem, a full separation of scales is adopted since the microscopic characteristic length scales ( i < m ) are much smaller than the macroscopic characteristic length scale L, which is typically the length over which the macroscopic fields vary over time i.e.
where m and i are the characteristic lengths of the microstructural components (matrix and inclusions, respectively) and L.
Mass diffusion In the mass diffusion problem, the separation of scales can be quantified based on the characteristic times associated with each material constituent. The characteristic times for the matrix t m and the inclusion t i can be written as where D m and D i are the mass diffusivity coefficients of the matrix and inclusions, respectively. In the present work, a relaxed separation of scales is considered for the diffusion problem, which is a special case of coupled scales. In the regime of relaxed separation of scales, the characteristic diffusion time of the matrix t m is very small compared to the one of the inclusion t i , and the macroscopic loading time T : A relaxed separation of scales is applicable to the homogenization of mass diffusion problems in lithium-ion batteries, where the lithium ions diffuse essentially instantaneously through the electrolyte material (matrix) in contrast to the very slow diffusion in the active particles (inclusions). The relaxed separation of scales has a direct implication for the model reduction presented in "Model reduction leading to an enriched continuum" section, since it allows the decomposition of the microscopic solution fields into the steady-state and transient parts. The separation of scales also indicates whether the transient terms in the conservation laws at the micro-and macro-scales should be included or not. These conservation laws are stated next.

Conservation laws at micro and macroscales
Mass conservation the mass conservation at the macroscale reads: wherej andċ are the macroscopic mass flux and the macroscopic rate of change of the concentration field, respectively. To capture the time dependent mass diffusion behavior inside transient inclusions, the mass conservation is considered at the microscale: where j andċ are the mass flux and the rate of change of concentration at the microscale.

Conservation of linear momentum
For the considered problem, mechanical inertia can be neglected, for which the macroscopic linear momentum balance equation reads: where σ is the macroscopic stress tensor. Given the full separation of scales for the mechanical problem, the conservation of linear momentum at the microscale also does not include transient terms neither and reads where σ is the microscopic stress tensor. The constitutive equations for the macroscopic quantities σ,j andċ are yet unknown; in the computational homogenization, these are obtained through an upscaling procedure. The boundary and initial conditions at the macroscale are given by the particular problem at hand. At the microscale, the constitutive equation for σ, j and c are assumed to be known, as presented in "Coupled diffusion-mechanics formulation" section. The boundary conditions at the microscale are obtained by downscaling relations, which will be presented next.

Downscaling
In first-order computational homogenization, the microscopic fields are approximated as the first-order Taylor's series expansion around a macroscopic pointx. The chemical potential μ in a microscopic domain can then be written as whereμ and ∇μ are the macroscopic chemical potential and its gradient, respectively, andμ is the fluctuation field of the chemical potential at the microscale. The latter is due to the difference in material properties of the constituents, and the transient loading conditions at the macroscale. Similarly, the microscopic displacement field u can also be expressed as the first order Taylor's series expansion around a macroscopic pointx whereū and ∇ū are the macroscopic displacement field and its gradient, respectively, and u is the microfluctuation of the displacement field.
In computational homogenization, downscaling is referred to as the transfer of macroscopic quantities to the microscale, as shown in Fig. 1b. Macroscopic quantities which are to be transferred to the microscale depend on the physical phenomena under consideration. For instance, in first-order transient computational homogenization, both for diffusion processes [36] and dynamics [37,38], both the primary macroscopic field and its gradient are transferred to the microscale. In a steady-state/static computational homogenization scheme, only the gradient information needs to be transferred to the microscale.
In transient computational homogenization, the first constraint on the microscale solution is that the volume average of the microscopic primary field is enforced to be equal to the corresponding macroscopic field which, by using the definitions (37) and (38) a chosen positioning of the microscopic domain such that x −x = 0, requires that the average of the microfluctuations over the microscopic domain vanishes The second constraint on the microscopic solution fields is that the average of the microscopic gradient fields should be equal to the corresponding macroscopic gradients which by using Eqs. (37) and (38) and the identity ∇(x −x) = I can be written as The last terms in the above equations, i.e. the average of the gradient of the microfluctuation fields ∇μ and ∇ũ should vanish to satisfy the requirements (41). After applying Gauss's theorem, these can be written as where n is the outward unit-normal vector to the microscopic boundary ∂ with an infinitesimal surface area d .
Constraints (40) can be applied by prescribing the respective fields at one point in the microscopic domain, along with the elimination of rigid body motion, to the corresponding macroscopic field values. To apply constraints (43), specific types of boundary conditions are used at the microscale. Typical choices for these boundary conditions are (i) zero fluctuation boundary conditions or (ii) periodic fluctuation boundary conditions as used later in this work.

Upscaling
Next, we discuss the upscaling relations which provide the constitutive equations for the macroscopic quantities. In computational homogenization, upscaling refers to the transfer of information from the microscale to the macroscale by requiring equality of the macroscopic and volume averaged microscopic (virtual) powers, known as the (extended) Hill-Mandel conditions in the literature [31,37,39]. The microscopic primary field ansatz e.g. (37) and (38), is then injected in the expression of the virtual power average and the macroscopic quantities are obtained by applying proper boundary conditions.

Mass diffusion
The micro-macro scale equivalence of the virtual power due to mass diffusion Substituting the variation of the microscopic chemical potential δμ using (37) in the right hand side of (44) yields Rearranging the above expression for δμ and δμ yields The last term in the above expression, after applying the chain rule and the divergence theorem, reflects the weak form of the microfluctuation mass conservation The first term on the right hand side of the above expression is the weighted residual of the microscopic conservation of mass (34), whose solution at the microscale should vanish. For the prescribed zero microfluctuation boundary condition or the periodic boundary conditions, the second term also vanishes and Eq. (46) reduces to from where the macroscopic mass flux can be recognized as and the rate of change of the macroscopic concentration aṡ The volume averages in Eqs. (49) and (50) can also be converted to boundary integrals using the divergence theorem and the microscopic mass conservation (34) with j n = j · n the normal outward mass flux.
Mechanics In the absence of inertia effects, the standard Hill-Mandel condition applies for the homogenization of the mechanical problem. Following similar steps as described above, allows identification of the (standard) macroscopic stress which can be converted to an expression in terms of tractions at the microscopic boundary Once the solution to the microscopic problem (34) and (36) is known, the reaction mass fluxes j n and the reaction forces t n can be computed and post-processed to obtain the macroscopic quantitiesj,ċ and σ. Next, we discuss the solution procedure to obtain the reaction fluxes j n and forces t n through a reduced order model, rather than the fully resolved model of the microscopic domain. An alternative homogenization route is to average the dissipation, given in Eq. (7), at the microscale and equating it to an assumed macroscopic dissipation expression. For a firstorder computational homogenization approach, the ansatz in Eqs. (40) and (43) can be inserted into the microscopic dissipation. Expanding and applying the required boundary conditions to eliminate the fluctuation fields, the corresponding macroscopic quantities can be obtained along with the weak forms of the balance laws at the microscale.

Model reduction leading to an enriched continuum
In this section, a model reduction of the microscopic coupled diffusion-mechanics problem is presented. The microscopic chemical potential and displacement fields are first decomposed into their steady-state and transient parts and reduced bases are identified.
The reaction fluxes and tractions, which are required to compute the macroscopic quantities, are written in terms of the coefficients of these reduced bases. Next, the expressions for the macroscopic quantities are derived explicitly. Finally, an emergent macroscopic enriched-continuum formulation, which arises as a consequence of model reduction at the microscale, is presented.

Finite element discretization
Using the finite element discretization, the linear momentum balance (36), the mass conservation (34) and the constitutive models (27)-(29), the discretized coupled diffusionmechanics problem in terms of the unknown nodal values of the chemical potential μ and displacements u can be written as where K μμ , M μμ , K uu and K uμ are the mobility, mass, stiffness and coupling matrices, respectively, and K μu = − K uμ T . The right hand sides j n and t n are the vector of reaction fluxes and reaction forces.
In the computational homogenization framework, once the solution for the microscopic primary fields μ and u is known, the reaction fluxes j n and reaction forces t n can be computed. In a two-scale setting, this is an expensive task, especially in the transient regime, since it requires the solution of a coupled problem at each macroscopic material point at each time step. Hence, an approximate solution based on a model reduction technique is called for.
To apply the model reduction, instead of solving a coupled system of Eqs. (56)-(57), we first analyze each equation separately and then the coupling effect is taken into account when the reduced bases are constructed. The homogenization conditions in Eq. (40) amounts to kinematically constraint the microscale to the macroscopic pointx and requires the macroscopic chemical potentialμ to be the average value of the microscopic chemical potential field μ. In a discrete setting, it can be achieved by prescribing the microscopic fields μ and u degrees of freedom (DOF), at a point x in the microscopic domain, equal to the corresponding reference values of macroscopic fieldsμ andū. It is allowed to fix the displacement field and the chemical potential at a point in the microscale because all the material properties i.e. C, S, M and are independent of these solution fields. Also, the displacement field u at the microscale are defined up to the rigid body motion and the chemical potential μ can also be defined up to a constant since the microscopic flux j, given in Eq. (15), depends on the gradient of chemical potential ∇μ and linear momentum balance, given in Eq. (26), is not affected by adding a constant term to the chemical potential. In this study this point is chosen to be x 1 which is at the lower left corner of the rectangular microscopic domain. The constraint (43) is satisfied by applying the periodic boundary conditions on the fluctuation fieldsμ andũ. Due to the applied periodicity, the DOFs at the other three corner nodes, denoted as points x 2 , x 3 and x 4 , are also fully prescribed, while the rest of the DOFs in the microscopic domain are considered free. More details on applying the boundary conditions in a discrete setting for a scalar field like μ can be found in [31] and for a vector field like u in [32]. The discretized mass diffusion Eq. (56) partitioned into prescribed 'p' and free 'f ' degree of freedoms takes the form Similarly, the mechanical equation after partitioning into its prescribed and free DOF can be written as For the microscopic response, both the chemical potential μ and displacement u are next split into their steady-state and transient parts.

Microscopic fields decomposition
According to the relaxed separation of scales, the transient response of the system evolves independently from the steady-state one. The steady-state response depends on the macroscopic input parameters (μ,ū, ∇μ, and ∇ū) through the prescribed DOFs μ p and u p , whereas the transient response only affects the inclusions that are part of the free DOFs. (In a discrete setting this requires that the prescribed DOFs always reside in the matrix material so that the transient response can evolve independently.) Consequently, the free parts of the microscopic solution fields are decomposed into a steady-state and a transient part. The free part of the chemical potential field can be written as where μ f ss is the steady-state and μ f tr is the transient part. Since the mechanical response is coupled to that of the mass diffusion, the displacement field will also evolve in time due to the change of the chemical potential. The free part of the microscopic displacement field u f is also decomposed into its steady-state u Next, the steady-state and transient reduced bases have to be determined for both the chemical and mechanical fields.

Steady-state response
The steady-state part of the microscale solution follows the macroscale solution instantaneously. To obtain the steady-state response, the discrete systems of Eqs. (58) and (59) are written considering the steady-state contributions μ ss and u ss only.

Mass diffusion
Substituting the steady-state chemical potential field μ ss in the second line of Eq. (58) yields The steady-state part of the chemical potential μ ss can then be expressed in terms of the prescribed DOF μ p as where μμ is the Schur-complement. When multiplied with the macroscopic quantities S fp μμ provides the steady-state homogenized response for the linear material model and thus can be considered as the steady-state reduced basis for the chemical potential field.

Mechanics
Similarly, to obtain the steady-state displacement field u f ss the second line of Eq. (59) is considered Substituting expression (64) for μ f ss in Eq. (65) yields

Transient response
As stated in "Microscopic fields decomposition" section , due to the relaxed separation of scales, to identify the transient reduced basis it is justified to use the free DOFs only. From Eq. (58) with account for (64), the free part of the discrete mass conservation equation can be written as and from Eq. (59) with account for (65), the free part of the discrete conservation of linear momentum can be written as from where where * μ is the matrix containing the columns of the reduced transient functions k μ and N μ q is the number of reduced basis functions for the chemical potential, which is much smaller than the total number of the free DOF N f , i.e. N μ q N f . Similarly, the transient displacement field u tr can also be written in terms of the reduced basis functions k u and their corresponding coefficients η k u as where N u q is the number of reduced basis for the displacement field. We will show later that N μ q and N u q and η u and η μ are the same. The selection criteria for the set of N q basis functions will also be presented later. Next, the reduced basis functions in Eqs. (71) and (72) are identified using a spectral decomposition scheme.

Mass diffusion
where M μμ ff uμ is the coupled mass matrix. The mass conservation (34) is a parabolic partial differential equation which has a natural solution that decays exponentially in time, i.e. μ = k exp (−α k t), substituting it in Eq. (73) yields the eigenvalue problems where k is the k-th eigenvector and α k the associated k-th eigenvalue.
where η k can be interpreted as the modal amplitude, and η is a column of size N q .

Mechanics
Substituting the expression of μ with * u = S ff uμ * . Next, we reconstruct the total solution for the chemical potential and displacement fields from their respective transient and steady-state responses.

Linear superposition
Substituting the expressions for μ ss from (64) and μ tr from (77) into (60), the total chemical potential field at the microscale can be written as where only the reduced basis * is coupled to the microscopic mechanical problem via the coupled mass matrix M μμ ff appearing in the eigenvalue problem (74). Similarly, the total microscopic displacement field u can be reconstructed by substituting the expression for u ss from (67) and u tr from (78) into Eq. (61) i.e.
Both, the steady-state and transient parts of the microscopic displacement field are coupled to the chemical problem as the coupling matrix K uμ appears in the matrices S ff uμ , S ff uμ and * . Equations (79) and (80) shows that the microscopic solution fields, μ and u , are completely given by the chemical potential μ p , the displacement u p at the prescribed DOFs and the coefficients of the transient reduced basis η . Generally, in a two-scale setting, the microscopic fields at the prescribed DOFs, where the microfluctuationsμ andũ are zero, are given by the macroscopic quantities, as can be seen in Eqs. (37) and (38). Therefore, the only remaining unknown fields at the microscale are η which can be obtained from the evolution equation, derived in the next subsection.

Evolution equation
The time evolution of η can be obtained from the free part of Eq. (58) which after using the normalization conditions in (75) and (76) takes the form where Equation (84) is a set of N q decoupled ordinary differential equations (ODEs) which represent the reduced order model for the evolution of diffusion-mechanics behavior at the microscale. The right hand side of (84) acts as the forcing term to the set of ODEs in terms of macroscopic fields present in μ p and u p .

Reaction fluxes and forces
Next, we write the reaction mass fluxes j p n and tractions t p n in terms of the coefficients of the steady-state and transient bases functions.

Reaction fluxes
The reaction mass fluxes j p n can be obtained from the first line of the discrete mass conservation Eq. (58) Substituting the expressions for μ f and u f from Eqs. (79) and (80) Rearranging terms gives the resulting reaction mass flux (91)

Reaction Forces
Similarly, the first part of Eq. (59) provides the expression for the reaction forces t p n at the prescribed DOF Substituting the expressions for μ f and u f from (79) and (80) into (92), gives which after rearranging for terms can be written as where In the expressions of reaction fluxes (90) and reaction forces (94), the only unknown is η which needs to be solved for in combination with the evolution Eq. (84), while μ p and u p are written in terms of the given (prescribed) macroscopic quantities.

Macroscopic quantities
Next, the expressions for the macroscopic quantities σ,j andċ are derived in terms of macroscopic DOF, and the coefficients of the microscopic transient basis η .

Macroscopic flux
In the discretized form, the boundary integral (51) of the macroscopic fluxj can be written asj where x p = (x p −I px ), with I p is the column of ones with dimension (p×1). Substituting the expression for j p n from (90) in (96) gives

Macroscopic stress
Similarly, the expression (55) for the macroscopic stress σ in its discrete form can be written as Substituting the expression for the reaction forces t p n from Eq. (94) provides which after using the discretized μ p and u p from Eqs. (98) and (99) take the following form The coefficients in Eq. (107) are given by where x p T ⊗ K uu pp I p = 0 has been accounted for.

Mode selection criteria
The microscopic fields μ and u , given by Eqs. (79) and (80), can be fully described by the macroscopic fields (μ,ū), their gradients (∇μ, ∇ū) and the coefficients of the reduced bases η . The size of the original eigenvalue problem is equal to the number of free DOF N f present in the system, which provides the complete set of eigenvectors . Owing to the fact that in diffusion problems the lowest eigenvalues α k are the most important ones, the eigenvectors corresponding to the first (several hundreds) lowest eigenvalues could be taken as the reduced basis. However, this would still entails a computationally inefficient scheme, since in a two-scale setting, where η is solved at the macroscopic quadrature points as internal variables, solving hundreds of ordinary differential equations for the internal variables would still require noticeable computational efforts. Therefore, the reduced set of eigenvectors * can be extracted from by taking into account that the right hand side of (84) acts as the forcing term and the modal coordinate η k corresponding to the forcing terms with a higher magnitude will have a higher amplitude and therefore contribute more to the homogenized behavior at the macroscale. Substituting the expressions for the prescribed chemical potential μ p and displacement u p fields, Eqs. (98) and (99) , in the evolution Eq. (84) provides which after using the definition of the coupling terms in (101), (104) and (108) takes the following form where 0 C η = V ( 0 C η ), 1 M η = V ( 1 M η ) and 2 C η = V ( 2 C η ), and takes into account that K μu qp I p = 0 q . The coefficients 0 C η , 1 M η and 2 C η couple the microscopic transient behavior, in terms of η k , to the macroscale fields. The higher the value of a coefficient, the higher the contribution of the respective η k to the macroscale behavior. This information can be exploited to identify a reduced set of eigen vectors * . The eigenvectors associated to 0 Cη k with a relatively high contribution are identified using where | • | is the absolute value of •. Similarly, for each component of 1 Mη k it can be stated and for each component of 2 C η k using Then, a reduced eigenbasis * can be obtained by requiring a minimum threshold e • for a coefficient •, such that * For a macroscopic simulation during an offline stage, individual threshold value signifies the corresponding macroscopic quantity and it should be selected accordingly.

Macroscale enriched continuum
The model reduction at the microscale leads to an enriched continuum formulation, as shown in Fig. 1c Macroscopic concentration rate:ċ = 0 C η Tη + 0 C˙μμ + 1 C ∇μ · ∇μ + 2 C ∇u : ∇u T , Macroscopic momentum conservation: ∇ · σ = 0, Macroscopic stress: Internal variable evolution: Different solution methods can be adopted to solve the coupled diffusion-mechanics enriched continuum problem, depending on whether η is evaluated at the macroscopic quadrature points, leading to an internal variable solution scheme, or at the nodes along withμ andū, which leads to a multi-field solution scheme. Numerical analysis for the solution of the enriched continuum formulation will be discussed in a future contribution.

Numerical examples
In this section, the proposed reduced order homogenization for coupled transient diffusion-mechanics is analyzed at the microscale. First, the problem setting is presented. The coupled transient bases are identified. Then, the microscopic fields and macroscopic quantities computed with the reduced order homogenization are compared with those obtained through the expensive, fully resolved, conventional computational homogenization scheme. Finally, the computational efficiency of the proposed reduced order homogenization is assessed.

Problem setting
Lithium ion battery electrodes are majorly composed of two components: the electrolyte (matrix) and the active particles (inclusions). As an example, in this study a cathodeelectrolyte system is considered, in which the electrolyte is lithium hexa-fluoro phosphate (LiPF 6 ) and the embedded active particles are made of lithium cobalt oxide (LiCO 2 ). For simplicity, it is assumed that the active particles are surrounded by the electrolyte only. All the other materials, e.g. the polymer binders, conductive particles etc., are disregarded following similar simplifications made in [30,40]. The material and geometric parameters are listed in Table 1. All material properties are assumed to be constant and do not change with the chemical potential or stresses in the material. For a given material, the chemical modulus and the mobility coefficient M combine to form the diffusivity coefficient D = M of the material. The diffusivity D m of the lithium ions is much larger in the electrolyte as compared to the diffusivity D i in the active particles, indicating that the relaxed separation of scales (32) holds for the considered problem. It is assumed that the electrolyte material does not swell with the introduction of lithium ions. The active particles are spherical in shape, vary in size and are placed randomly in the electrolyte which creates a poly-disperse heterogeneous medium [23]. In this example (for simplicity reasons), we consider a two dimensional mono-dispersed heterogeneous medium, as shown in Fig. 2, which is generated by a level set based random sequential adsorption method [41].  [28] i = k b T 0 /c 0 10,202 Matrix diffusivity [28] D m 6 × 10 −11 (m 2 s −1 ) Inclusion diffusivity [28] D  In the simulations, all the parameters were non-dimensionalized, the time was normalized with respect to the characteristic diffusion time of the inclusion, i.e.t = t t i , the lengths are normalized with respect to the characteristic length of the microscopic domain, i.e.x = x , the chemical potential is normalized with respect to the maximum attainable chemical potential in the inclusionμ = μ μ max , where μ max = i (c max − c 0 ), the displacement field is normalized with respect to the characteristic length of the microscopic domain, i.e.û = u and the stresses are normalized with respected to the Young's modulus of the inclusion E i .
The microscopic domain is excited by the chemical loading given in terms of the macroscopic chemical potentialμ and the gradient of the macroscopic chemical potential ∇μ as a function of time i.e.
where ω = 2π T is the angular loading frequency, T is the time of one period andμ = μ max and ∇μ = 0.1μ max are assumed. Externally applied mechanical loads to the microscale are neglected here. At the microscale, periodic boundary conditions are used to satisfy the Hill-Mandel conditions for both the mass diffusion and mechanical problems. The microscopic domain, shown in Fig. 2b, is discretized with linear triangular finite elements. For time integration, the backward-Euler method was used with a time step t = 1 × 10 −3 T [s].

Reduced basis identification
After assembling the finite element matrices and applying the boundary conditions at the microscale, the first step of the reduced order homogenization is the solution of the coupled eigenvalue problem (74). This eigenvalue problem is solved for the first two hundred smallest eigenvalues and the corresponding eigenvectors . Then, using the mode selection criteria given by Eq. (114), the reduced basis * ∈ is based on the coupling terms 0 Cη k , 1 Mη k and 2 C η k with the threshold value e C = e M = e C = 0.1. The number of eigenvectors selected in the eigenbasis * depends on the topology of the micro-structure, the strength of the coupling in diffusion-mechanics and the material contrast between the matrix and the inclusions. For each selected eigenvector k μ , there is a corresponding coupled mechanical eigenvector k u = S ff uμ k μ , both * μ and * u are shown for the considered micro-structure in a coupled manner in Fig. 3. The ten modes selected are not the modes corresponding to the 10 consecutive smallest eigenvalues.
The inclusions swell where the chemical potential is high, indicated by the red regions inside the domain and the inclusions shrink where the chemical potential is low, indicated by the blue regions. The modes have contributions in the inclusions only, in accordance with the relaxed separation of scales. If the material properties do not fulfill the requirements of relaxed separation of scales (32)  in the matrix and consequently the proposed reduced homogenization method will not capture the phenomena adequately. Among the selected eigenvectors, shown in Fig. 3, there is one eigenvector with the highest relative importance of the eigenmodes in terms of their contribution to the macroscale, which can be quantified by a measure ξ k where | • | and || • || are the absolute value and Frobenius norm of a quantity •. In this example, eigenvector numbered 7 has the highest contribution to the macroscale, while other eigenvectors in * make only small improvements in capturing the phenomena. For a more detailed mode selection analysis, in the case of heat diffusion, the reader is referred to [31].

Microscale simulations
Next, we compare the microscopic fields computed by the model reduction method and the (expensive) fully resolved finite element calculations. For the fully resolved finite element analysis, the coupled system of Eqs. (56)-(57) is solved for μ and u , directly on the finite element mesh of the considered RVE 2b. For the reduced model, the coefficients η are solved by using Eq. (110); subsequently, the microscopic fields μ and u , are reconstructed (localization operation) by post-processing through Eqs. (79) and (80), respectively. Note that, in a two-scale simulation the post-processing of microscopic fields μ and u , is generally not done, unless the microscopic fields are also the quantities of interest in addition to the macroscopic field. Figure 4 shows the contour plots of the normalized chemical potentialμ and Fig. 5 shows the contour plots of the normalized hydrostatic stresŝ σ hyd =σ 11 +σ 22 3 at timet = 250 t . The minor differences between the fully resolved solution and the reduced order model are due to the approximate nature of the model reduction.

Effective macroscopic quantities
Next, we compare the macroscopic quantitiesj,ċ, σ computed with the conventional transient homogenization and the developed reduced order homogenization. For conventional computational homogenization, the fully resolved finite element analysis of the coupled system of Eqs. (110) is solved for η . Once η is known, the macroscopic quantities are calculated directly from the expressions (100), (103) and (107) for the macroscopic mass fluxj, the macroscopic concentration rateċ and the macroscopic stress σ, respectively. Figure 6 shows the time evolution of the macroscopic quantities computed with the (expensive) conventional transient computational homogenization (CTH) method (shown in red) and the proposed inexpensive reduced computational homogenization (RTH) method (shown in blue). The reduced order homogenization method shows an excellent approximation without any noticeable discrepancies.
The computational gains achieved with the reduced model are substantial. Neglecting the off-line stage, motivated by the fact that for a specific microstructure and set of material parameters the off-line stage only needs to be performed once. Using Matlab 2018b on a computer with Core-i7 4.4 GHz processor and 16Gb memory, for the conventional computational homogenization the coupled problem (56)-(57) takes approximately 5000 times more computational time than the solution of the uncoupled ordinary differential Eq. (110). Next, we asses the proposed reduced model with different value of coupling coefficient γ and the microscopic domain size .

Diffusion-mechanics coupling effect
In diffusion-mechanics, the coupling is governed by the partial molar volume parameter γ in Eq. (28). The higher the value of γ , the higher the coupling will be. The constitutive Eq.  Fig. 7 The effect of the coupling term γ on a the macroscopic rate termċ and b the macroscopic equivalent stress σ eq . The equivalent stress is calculated as σ eq = σ 2 11 + σ 2 22 + 2σ 2 12 can become non-positive definite, which will make the eigenvalues α (76) equal to zero or even negative. However, for the realistic material properties of the cathode in lithium ion batteries, this is not a problem since the upper limit for γ is 21 × 10 −6 (mol m −3 ), which is much greater than the physical value of γ = 3.497 × 10 −6 (mol m −3 ). Figure 7 shows the effect of increasing the γ value. The macroscopic rate termċ and the macroscopic stress σ increase as γ increases in accordance with their microscopic counterparts (23) and (24), respectively. The proposed model reduction scheme captures the full finite element solution very well, and hence for clarity the finite element solution is not shown anymore. Next, we analyze the effect of the microscopic domain size on the macroscopic quantities.

Size effect
To measure microscopic size effect on the macroscopic quantities, the material parameters are kept the same and the characteristic size of the microscopic domain is changed while keeping the inclusions size i the same or scaling it along with the microscopic size. In the first case, the macroscopic quantities do not vary with the changing RVE size. However, if the inclusions size i is scaled along with the characteristic microscopic size the microscopic and the macroscopic quantities change. For the later case, the normalized characteristic length of the microscale is changed fromˆ = 1 × 10 −2 toˆ = 1 × 10 1 and For the larger microstructural sizes the response is clearly size dependent. In particular, the macroscopic stresses σ are much higher for small microstructures compared to the large ones. This is due to the coupling effect and the diffusion rate. For smaller sizes, the chemical potential (and the concentration of a species) increases in the inclusion domain, which causes the inclusions to swell and produce higher stresses on average. Conversely, when the microscale size is increased, within the same time period, the mass diffusion happens to the outer layer of the inclusions, swelling only that part of the inclusions, which creates higher local stresses, as can be seen in Fig. 9. However, due to overall increase in the volume of the microscopic domain and the small relative volume of high stress regions, the macroscopic volume average stresses decrease, as can be observed in Fig. 8.

Conclusions
In this work, a model reduction based homogenization technique for coupled diffusionmechanics problems has been presented. A formulation based on the chemical potential and linear strain field is derived, which eases the implementation since it only requires a C 0 -continuous discretization. This is in contrast with the conventional formulation used in diffusion-mechanics based on the concentration and strain fields, which requires a