Two-scale analysis of transient diffusion problems through a homogenized enriched continuum

Abstract This paper addresses the two-scale problem underlying the enriched continuum for transient diffusion problems, which was previously developed and tested at the single scale level only (Waseem et al., Comp.Mech, 65, 2020). For a linear material model exhibiting a relaxed separation of scales, a model reduction was proposed at the micro-scale that replaces the micro-scale problem with a set of uncoupled ordinary differential equations (ODEs). At the macro-scale, the balance law, the ODEs and the macroscopic constitutive equations collectively represent an enriched continuum description. Examining different discretization techniques, distinct solution methods are presented for the macro-scale enriched continuum. Proof-of-principle examples are solved for a mass diffusion system in which species diffuse slower in the inclusion than in the matrix. The results from the enriched continuum formulation are compared with the computational transient homogenization (CTH) and direct numerical simulations (DNS). Without compromising the solution accuracy, significant computational gains are obtained through the enriched continuum approach.


Introduction
Composite materials are used in a wide range of applications from large structures to micrometer size components [1]. The contrast between the properties of the constituents may differ by several orders of magnitude [2]. In general, it is computationally very expensive to simulate such heterogeneous materials at the micro-scale since the analysis involves a large number of degrees of freedom to capture the physical phenomena. Homogenization can be used instead, which provides the effective properties through an averaging procedure. Various homogenization methods can be found in the literature, i.e. analytical bounds [3,4], self-consistent methods [5] and asymptotic homogenization [6,7,8]. Most of these homogenization methods assume a steady-state condition at the micro-scale to calculate the effective properties. However, in reality, due to the high contrast in material properties and time varying loading conditions, this assumption is often not valid. For example, the steady-state assumption is not appropriate in the case of diffusing species in polycrystaline materials, where the diffusion process in the grain boundaries reaches steady-state while the grain interior still remains transient [9]. Similar phenomena are observed for fluid and solute transport in geo-materials [10], diffusion in porous gels [11] and diffusion of Lithium ions in electrode-electrolyte systems of Li-ion batteries [12].
Computational homogenization (CH) [13,14,15,16] is a more recent and robust homogenization procedure, which has been extended to transient problems [17], including diffusion [18,19]. For a detailed review and perspectives of CH see [20,21,22]. Computational homogenization for transient problems, despite its robustness, suffers from a high computational cost. It requires the solution of a micro-scale problem at each macroscopic material point, which in transient problems must be solved at each time increment. To circumvent this problem, reduced models for the effective diffusion response in transient regimes have recently been proposed [23,24], in which a chemical creep function is obtained. This chemical creep function behaves similar to phenomenological models of viscoelasticity and consists of an infinite number of Kelvin-Voigt units. A reduced model is then obtained by selecting a finite number of Kelvin-Voigt units. In [24], an analytical expression is obtained for a transient inclusion, embedded in an infinitely fast matrix, which is only possible for certain inclusion morphologies.
In a previous paper, an alternative reduced order model was proposed for the micro-scale problem, first for dynamics [25] and then for diffusion problems in the context of heat conduction [26]. It is based on computational homogenization and applies a reduction technique to the whole unit-cell (as opposed to the inclusion only in [24]). This numerical approach accounts for different complex RVE morphologies. It is applicable in the relaxed separation of scales regime in which the characteristic diffusion time τ m of species in the matrix is much smaller than the characteristic diffusion time τ i in the inclusions. Moreover, the macroscopic loading time T is such that the matrix can be assumed in a steady-state regime while the inclusion exhibits transient inertia effects. The relaxed separation of scales can be expressed in terms of the characteristic diffusion times of the constituents and the characteristic loading time as follows For more discussion on the relaxed separation of scales in diffusion problems, see [8,24]. Assuming a linear material model for both the matrix and inclusions, allows to perform an additive decomposition of the microscopic solution field into a steady-state and a transient part. The reduced model is defined by means of a static condensation of the steady-state part of the response, and by projecting the transient part onto a reduced (eigen)basis. The size of the microscopic system of equations is thereby reduced from N coupled finite element degrees of freedom (d.o.f) to just a few N q decoupled d.o.f with N q N . An averaging is performed to obtain the macroscopic constitutive equations in terms of the coefficients of the reduced basis and the emerging macroscopic (internal) variables. The resulting set of equations to be solved at the macro-scale describes what is called an enriched-continuum, consisting of the macroscopic mass balance, constitutive equations obtained by homogenization and an evolution equation for the enrichment variables. The coefficients of the micro-scale reduced basis are called enrichment-variables and the vectors, constituting the micro-scale reduced basis are called enrichment-functions. For a certain microscopic domain and material properties, the enrichment-functions and the coupling terms in the macroscopic constitutive equations are computed once and for all, and the subsequent on-line computation only consists in solving N q decoupled ordinary differential equations along with the macroscopic mass balance equation. Such an enriched-continuum formulation can be applied to a heterogeneous medium which exhibits transient diffusion phenomena within the assumption of linear material properties and the relaxed separation of scales, for example mass diffusion problems in batteries where these assumptions are valid in application some ranges.
The enriched continuum formulation for (heat) diffusion problems was developed in [26] and validated at the micro-scale unit-cell level against full transient computational homogenization results. No macro-scale enriched continuum simulations were used in [26]. The aim of this work is to demonstrate the applicability of the reduced order framework to the solution of transient problems with an underlying microstructure obeying the relaxed separation of scales. Since the resulting macroscopic problem consists of an enriched continuum with enriched variables and corresponding evolution equations, dedicated solution techniques have to be adopted. To this end, two solution methods based on different spatial discretization schemes are analyzed. First, a multi-field solution method is presented, in which both the enrichment-variables and the primary macroscopic field are computed on the nodes and interpolated using finite element shape functions. Next, an internal-variable solution method is presented, in which the enrichment-variables are evaluated at the macroscopic quadrature points and considered as internal-variables, while the macroscopic primary field variable is interpolated in a classical manner using the finite element shape functions, as shown in Figure 1. The advantage of the multi-field solution method is that the macroscopic primary variable and the enrichment-variables are solved in a coupled manner at each time step; the disadvantage is the significantly increased size of the finite element system of equations. In the internalvariable method, the enrichment-variables are eliminated, locally at the Gauss quadrature point, by condensation. Once the solution of the macroscopic variable field is available, the enrichmentvariables are evaluated and stored for the next time increment. This allows to use different time integration schemes for the enrichment variables as opposed to the primary field, e.g. to better capture micro-scale transient effects. The local condensation of the enrichment-variables reduces the size of finite element system of equations to be solved globally. However, on the downside the internal flux vector needs to be evaluated at each time step, in the internalvariable method since it depends on the values of the enrichment-variables from the previous time step. This again makes it less efficient. These two solution methods will be assessed in this paper. The novel contributions of the work presented here are: • the two-scale enriched macroscopic continuum implementation emerging from the model order reduction approach, applied to the solution of diffusion boundary value problems with incorporation of micro-structural transient effects; • the proposed different solution methods for the enriched-continuum for transient diffusion problems; • a comparison between the two-scale diffusion enriched-continuum macroscopic solution, the conventional transient homogenization and direct numerical simulations, allowing to evaluate the associated computational cost reduction.

Outline
Section 2 summarizes the enriched-continuum equations with expressions for the effective mass flux and the effective concentration rate. Section 3 presents the analyzed solution methods for the macro-scale enriched-continuum. First, the time integration schemes are presented and, subsequently, the multi-field and the internal-variable solution methods are derived. In section 4, proof-of-concept numerical examples are solved for species diffusion in a material with high contrast material properties. This will demonstrate the ability of the enriched-continuum formulation to reliably reproduce the results obtained by conventional transient homogenization (CTH) or direct numerical simulations (DNS), which are both will be shown to be outperformed in terms of computational efficiency.

Symbols and Notations
The (homogenized) macroscopic quantities are shown with a bar on top, e.g. a scalar, a vector and a second-order tensor are represented with a, a, A. There is no bar on top of a quantity belonging to the heterogeneous micro-scale problem, i.e., a scalar, a vector and a second-order tensor are written as a, a, A. The same Cartesian coordinate system is used at the micro-and macro-scale. Standard calculus operators are used in this work. For example, the dot product between two vectors is a·b = a i b i , between a second-order tensor and a vector is The tensorial dyadic product between two vectors is a ⊗ b = a i b j e i ⊗ e j and between a secondorder tensor and a vector is where e i represent the vectors of the Cartesian basis. The gradient of a scalar, a vector and a second-order tensor is defined as ∇a := ∂a /∂x i e i , ∇a := ∂a i /∂x j e i ⊗ e j and ∇A := ∂A ij /∂x k e i ⊗ e j ⊗ e k , respectively. Similarly, the divergence operates as ∇ · a := ∂a i /∂x i and ∇ · A := ∂A ij /∂x i e j . For linear algebra notation, columns are represented with a tilde underneath a lowercase letter, e.g. a , and matrices with a bar underneath an uppercase letter, e.g. 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 Ω .

Enriched Continuum Formulation
In this section, the transient mass diffusion problem is presented briefly, first for a heterogeneous domain and then for a homogenized domain, for which the two-scale computational homogenization and the resulting enriched continuum will be considered.

Direct Numerical Simulation
The diffusion of species in solids is governed by the mass balance equation, where,ċ is the rate of change of the concentration field c, j = −M · ∇µ is the mass flux, µ = Λ(c − c 0 ) is the chemical potential, Λ is the chemical modulus, c 0 is the reference concentration of a specie, M is the mobility tensor; Ω is the heterogeneous macro-scale domain, as shown in Figure 2(a). The numerical procedure to solve the diffusion equation (4) in Ω, with corresponding initial and boundary conditions is referred to as the direct numerical simulation (DNS). In here, ∂Ω µ and ∂Ω j are the Dirichlet and Neumann sub parts of the boundary ∂Ω, respectively, and its n outward unit-normal vector.

Computational Homogenization
Given its exorbitant computational cost, the DNS problem is commonly replaced by a computationally homogenized problem which represents an equivalent homogeneous problem by solving a coupled two-scale (macro-scale and a micro-scale) problem, as shown in Figure 2(b). The heterogeneous domain, in Figure 2(a), is replaced by a homogeneous one Ω, to which at each point, a microscopic domain Ω is attached. This microscopic domain is a representative volume element (RVE) or a unit-cell in the case of a periodic medium. In transient computational homogenization, at the macro-scale a transient mass balance equation is solved which is complemented by the initial and boundary conditions i.e.
Where, ∂Ωμ and ∂Ω j are the Dirichlet and Neumann sub parts of macroscopic boundary ∂Ω, respectively, n its outward unit-normal vector,ĵ is the prescribed mass influx andμ is the prescribed chemical potential. The constitutive responses for the macroscopic mass flux j and the macroscopic rate of change of concentration fieldċ are obtained by solving a micro-scale problem for the microscopic chemical potential field µ, which in the first-order computational homogenization is based on the following ansatz where,μ is the fluctuation in the chemical potential at the micro-scale due to the difference in material properties of the microconstituents properties and transient loading conditions, imposed through the macroscopic chemical potentialμ and its gradient ∇μ, see Figure 2(b). The transient mass diffusion problem which is solved at the micro-scale is accompanied by specific boundary conditions that satisfy the extended Hill-Mandel condition. In computational homogenization it is assumed that the average of the microscopic chemical potential equals to the macroscopic chemical potential and the gradient of the microscopic chemical potential to 6 be equal to the gradient of the macroscopic chemical potential.
The constitutive response of the microscopic mass flux j = −M · ∇µ and µ = Λ(c − c 0 ) are the same as used in the DNS. By using the equivalence of virtual power between the micro-scale and the macroscopic material point x (the extended Hill-Mandel conditions) the macroscopic flux and concentration can be obtained from the micro-scale fields as Despite of its robustness, computational homogenization is still very expensive in solving a transient mass diffusion problem, since it requires the solution of a microscopic problem (8) at each macroscopic point x at each time instance (increment). For more details on the transient computational homogenization and its solution schemes, the reader is referred to [18,19].

Enriched Continuum
To alleviate the computational cost of the transient computational homogenization scheme, for linear problems, a model reduction technique was presented for heat conduction in [26], which provides an enriched continuum at the macro-scale. It replaces the microscopic problem (8) with a set of uncoupled ordinary differential equations by projecting the weak Galerkin form of (8) onto an orthogonal set of reduced basis functions Φ , here called enrichment functions. The microscopic chemical potential field is written in terms of the macroscopic quantities and the enrichment functions as where η (k) , called enrichment variables, are the coefficients of the corresponding enrichment function Φ (k) . N q is the number of enrichment functions used, which is much smaller than the number of degrees of freedom N used to solve the microscopic problem (8) in computational homogenization, i.e. N q N . S (p) represents the steady-state counter part of the Φ (k) and it is obtained through static-condensation. N p N are the prescribed d.o.f. in the discretized domain Ω . This reduction in the number of the degrees of freedom at each macroscopic point x provides a significant reduction in computational time. The effective macroscopic transient mass balance equation remains consistent with the one used in transient computational homogenization scheme, i.e.equation (6). However, the expression for the macroscopic mass flux j and the macroscopic rate of change of concentration fieldċ are now written in terms of the (reduced) enrichment variables η i.e.
The coefficients in (11) and (12) are computed only once for a given microstructure in an off-line computation stage and used at the macroscale during the on-line stage where µ and η are solved for. A concise derivation of the expressions for the coefficients in equation (11) and (12) are provided in Appendix A. For more details on the enriched continuum formulation and its derivation the readers are directed to our prior work on heat conduction [26]. Through model reduction of the micro-scale problem, evolution equations for η emerge, which are a set of (uncoupled) ordinary differential equations, given bẏ where α is the diagonal matrix of size (N q ×N q ) containing the eigenvalues which are obtained by solving a generalized eigenvalue problem at the micro-scale during the off-line stage, expressions for * a and * d are given in Appendix A. Equations (6), (11), (12) and (13) represent the enrichedcontinuum model capturing the diffusion problem in the relaxed separation of scales regime (1). Next, we discuss the solution methods for the set of equations for the enriched-continuum.

Solution Methods
In this section, first the time integration schemes are presented for the macroscopic variableμ and the enrichment variables η . Then, the multi-field and the internal-variable solution methods are derived.

Time Integration Schemes
The rate of change of the macroscopic chemical potentialμ is integrated in time using a backward-Euler time integration scheme. For example, the general equatioṅ can be discretized asμ In the enriched-continuum formulation, the microscopic transient effects are captured at the macro-scale through the enrichment variables η . This allows to compute the transient effects, present at the micro-scale, on a coarse macroscopic mesh. In the internal-variable solution method, the internal-variables η are condensed out of the final system of equations, therefore different time integration schemes can be used for equations (6) and (13). Following the literature on visco-elasticity [27], as recalled in [28], a one-step second-order accurate time integrator can be obtained by writing equation (13) in a convolution form Using the semi-group property of an exponential and the additive property of an integral, the approximation in (16) can be written as A midpoint approximation of the integral in equation (17) provides a one-step second-order accurate time integration scheme to evaluate η n+1 , which, as argued in [28], is unconditionally stable. This second-order accurate scheme will be compared to the first-order accurate backward-Euler scheme for which the approximation for η n+1 is given by

Multi-Field Method
In a multi-field approach, the primary fieldμ as well as the enriched field η , are discretized on a finite element mesh. A combined weak residual Q(μ, δμ, η , δη ) is formulated by multiplying the macroscopic mass balance equation (6) and the evolution equation (13) with the appropriate admissible test functions δμ and δη , respectively where, R(μ, δμ) and E(η , δη ) are the individual residuals for equations (6) and (13). Integrating by parts and applying the divergence theorem on R(μ, δμ), provides the weak form as follows whereĵ = −j · n is the mass flux through the macroscopic boundary ∂Ω j , n is the outward unit normal, dv and da is the small differential volume and area elements associated to the domain Ω. Substituting the constitutive expressions (11) and (12) for the macroscopic flux j and the macroscopic chemical concentration rateċ yields a coupled system in terms ofμ and η , Q(μ, δμ; η , δη ) = Ω ∇δμ · a Tη + B · ∇μ + cμ + C · ∇μ dv

Finite element implementation
In the multi-field solution method, the macroscopic variableμ and the enrichment-variables η are interpolated using the finite element nodal shape functions where, N I is the shape function value associated to the I-th node of an element and M is the number of nodes in an element. N μ and N η are the column and matrix, respectively, containing the respective shape functions andμ is the column of the degrees of freedom ofμ, Nq } is the column vector of degrees of freedom of the vector fields η containing the nodal enriched degree of freedoms. Equation (22), can then be written in discrete form as The matrices in the above equation are given by Taking into account that equation (24) should hold for all admissible δμ and δ{η} provides In the multi-field method,{η} appears in the final equation together withμ, which suggests the same time integration scheme for both variables. Selecting the backward-Euler time integration scheme (15) for bothμ and{η}, we obtain Rearranging terms to gather all the unknowns at time t n+1 on the left hand side, the coupled system of equations to be solved at each time step can be written as

Internal-Variable method
In the internal-variable method, the macroscopic chemical potentialμ is interpolated using finite element nodal shape functions, while the evolution equation (13) forη is integrated at the macroscopic Gauss quadrature points. The residual R(μ, δμ) is built by multiplying equation (6) with an admissible test function δμ which, after applying integration by parts and the divergence theorem, takes the following form Substituting the expressions for the macroscopic flux from equation (11) and the concentration rate from equation (12), the weak form of equation (30) can then be written as

Finite element implementation
Here, as an example, the derivation is performed using the backward-Euler time integration scheme for bothμ and η ; a similar derivation applies to the integral form given by equation (18). Using the finite element discretization (23) only forμ, equation (31) takes the following form where the matrices are calculated using Next, substitution of the expression for η n+1 from (19) into the above equation leads to where the finite element matrices are recognized as here A = (I + ∆tα ) −1 . Finally, taking into account that equation (33) should hold for all δµ , the system of linear equations to be solved at each time step is given by Once the solution for the chemical potential fieldμ n+1 is known, η n+1 is calculated using equation (19) (or (18) if the second-order accurate scheme is used forη ). By substituting the expression of the time discretized enrichment-variables η n+1 in equation (31), these are condensed out and the resulting finite element system of equations (34) does not contain η as unknowns. This is the advantage of the internal-variable method, since the number of unknowns is significantly smaller as compared to the multi-field method. However, the downside of the internal-variable method is that the internal flux vector F int n , which contains the enrichment-variables η n defined at previous time step in the matrices A n , D n , A 3 n and A 6 n , has to be evaluated at each time increment, which as compared to the multi-field method, is not very efficient.
Accordingly, the internal-variable method and the multi-field method are more or less equivalent for moderate problem sizes. The internal-variable method satisfies the evolution equation for η (19) in an exact manner locally at the Gauss quadrature points, whereas the multi-field method satisfies equation (19) in an average, global, sense over the whole macroscale domain Ω. To demonstrate the equivalence, piece-wise constant finite element shape functions N η can be used for the enrichment-variables and the condensation can be performed to eliminate {η} n+1 from equation (28). Thereafter, the obtained system of equations consists of the same terms as in (34), with η weakly approximated at the nodes.

Numerical Examples
Numerical examples are conducted for species diffusion in a heterogeneous material with high contrast in properties, with slower diffusion in the inclusion than in the matrix. The microscopic domain is chosen to be a square unit-cell with a single centered circular inclusion, as shown in  (11) and (12) in the off-line stage. The material properties and geometrical parameters are provided in Table 1.
µ max are the maximum achievable chemical potential in the matrix, given byμ max = Λ(c max −c 0 ), where Λ = k b T /c 0 is the chemical modulus, c max and c 0 = 0.19c max is the maximum and reference concentration values in the matrix, respectively. Using a non-zero parameter 'a' in equation (35) implies that the chemical potential field can be linearly varied along ∂Ωμ; in (35) ω = 2π /T is the angular loading frequency and T is the total loading time. The material properties, length scales and the loading conditions are such that the relaxed separation of scales (1) is satisfied. The default parameters for the geometry, material and mesh are given in Table 1. At the micro-scale, a converged unit-cell mesh consisting of ∼ 4400 linear triangular elements (∼ 2200 nodes) was used. The DNS domain was discretized accordingly, with approximately the same number of linear triangular elements in each unit-cell. The macroscopic homogenized domain Ω, replacing the DNS domain, is discretized with rectangular four node bi-linear elements for bothμ and η in the multi-field method, or forμ only in the internal-variable method. The effect of the macroscopic mesh size on the homogenized results will be discussed in section 4.4. The smallest and the largest element size in the macroscopic domain are and 10 , respectively. A default time step size of ∆t = 3.6[s] is used in all the simulations unless otherwise stated.

Enrichment-Functions Selection and Solution Accuracy
To obtain the enrichment-functions and construct the coefficient terms for equations (11), (12) and (13), according to the expressions given in Appendix A, an eigenvalue problem with periodic boundary conditions is solved on the unit-cell. The eigenvectors corresponding to the relatively high values of the coupling terms d and a constitute the reduced basis set. The selected enrichment-functions Φ (k) for the considered unit-cell are shown in Figure 4. These six enrichment-functions will be used in the simulations, unless otherwise stated.   Figure 4: The selected enrichment-functions Φ (k) obtained by solving the generalized eigenvalue problem (44). An associated decay/rise time τ (k) = 1 /α (k) , with α (k) the k-th eigenvalue, is also given for each enrichment function Φ (k) .
A detailed discussion on the selection of the reduced basis set and the accuracy of the micro-scale solution with respect to the increasing size of the reduced basis set can be found in previous work [26]. Here, the focus will be on the accuracy of the macro-scale solution obtained with the two methods. The accuracy of the macro-scale chemical potential field with respect to the number of enrichment-functions is analyzed here. The macroscopic chemical potential fieldμ (Nq) obtained using the complete reduced basis set Φ with N q = 6, as shown in Figure  4, is considered as the reference solution. Calculations were performed, for both the internalvariable and multi-field methods, on a mesh of nel x 1 × nel x 2 = 30 × 3 elements, where nel x 1 and nel x 2 are the number of elements in x 1 and x 2 directions, respectively. As shown in Figure  5, both the internal-variable and multi-field method perform equally well. With the increase of the number of enrichment-functions the accuracy increases as well. Figure 5(a) shows the time evolution of the relative L 2 -error, while Figure 5(b) shows the time averaged relative L 2 -error.

Accuracy of the Time Integration Schemes for the Internal-Variable Method
Different time integration schemes and associated convergence analyses have been studied in the literature for parabolic PDEs [29] and viscoelastic materials [28]. However, the present emergent enriched-continuum is an unusual description of a macroscopic continuum, since the enrichment-variables η appear both in the diffusion and (more importantly) in the capacitance terms as well. Therefore, it is justified to investigate the convergence behavior of different time integration schemes for this enriched continuum. In the internal-variable method, due to local condensation of η at the Gauss quadrature points, the time integration forμ andη can be chosen independently. In section 3, two different time integration schemes for the approximation ofη were presented, i.e. a second-order accurate approximation in equation (18) and a first-order accurate approximation in equation (19). In this regard, simulations were performed using the internal-variable method on a mesh of nel x 1 × nel x 2 = 50 × 5 elements. For each time integration scheme, a reference solution η   , where D ef f 11 is the first component of the effective macroscopic diffusivity tensor D ef f = ΛB, the effective mobility B is given in (46). The results are shown in Figure 6. As expected, the backward-Euler first-order time integration scheme converges with a rate of convergence equal to one, while the time integration scheme given in equation (18) converges with a rate of convergence of almost two, providing a more accurate approximation ofη . However, the second-order method changes towards to first-order accuracy, which due to the backward-Euler time integration scheme used at the macroscale which limits the overall convergence rate. For example, for the default time step size of ∆t = 3.6[s], the second-order scheme provides three orders more accurate result than the backward-Euler time integration method. Since the second-order scheme given by equation (18) is a one-step time integration scheme, no extra memory or computational costs are associated with it. Therefore, it is suggested to use the second-order accurate time integration scheme whenever possible to capture the transient effects more accurately.

Homogenized Solution
Next, the results computed using the homogenized enriched continuum formulation are compared to the direct numerical simulations (DNS). To do so, the DNS solution µ(x, t) is averaged over the heterogeneous unit-cell sized domain λ around the point x = (2.5, 2.5) × 10 −2 [m], as shown in Figure 3 with a small square. In this example, the homogenized domain was discretized with nel x 1 × nel x 2 = 10 × 3 elements, with four Gauss quadrature points (ngp) per element. Figure 7 shows the time evolution of the macroscopic concentration rateċ(x, t), given by equation (12), and the averaged DNS solution μ Λ λ . These results are in perfect agreement with each other, which demonstrates the ability of the enriched-continuum formulation to capture the transient DNS solution.

Macroscopic Discretization Effect
In this section, the effect of the macroscopic mesh size (in both x 1 -and x 2 -directions), on the accuracy of the macroscopic solution,μ is examined. The DNS solution is taken as the reference. The homogenization level 'h' is defined as the ratio between the number of unitcells 'nuc' composing the DNS domain to the number of integration points to be solved in the homogenized domain. This provides an indicator for the trade-off between the accuracy and the computational cost. A higher h value indicates more homogenization, and thus lower computational costs. In a two-dimensional setting, with a quadrilateral discretization, the homogenization level can be written as In the subsequent simulations, the parameters, nuc = 1000 and ngp = 4 per element are fixed. The number of elements in x 1 -and x 2 -directions, nel x 1 and nel x 2 , are varied, resulting in a different value of h. In figure 8, the heterogeneous chemical potential field µ, calculated with the direct numerical simulation (DNS) is compared with the homogenized chemical potential field µ, obtained by the enriched-continuum formulation solved with the multi-field method. The homogenization level h given in equation (37) is increased by decreasing the number of elements in the macroscopic domain, while keeping the ratio nel x 2 = At h = 1, where the number of unit-cells is the same in both the heterogeneous and homogenized domains, the enriched-continuum response captures the DNS almost perfectly. In this example, due to a non-homogeneous boundary condition on ∂Ω µ , as given in equation (35), the homogenization problem requires high finite element density in x 2 -direction in the vicinity of ∂Ω µ , lower element density is needed in x 1 -direction.

Computational Efficiency
Next, the efficiency of the proposed approach compared to the direct numerical simulations (DNS) and conventional computational transient homogenization (CTH) scheme is investigated. All simulations were performed on a computer with a core-i5@3.2GHz processor and 8Gb RAM using Matlab 2018b. It can be seen in Figure 9, that the CPU time increases for all the homogenization methods as h → 1. The CTH is the most expensive among the considered methods, since (when implemented in its standard form) for each time increment at each macroscopic Gauss quadrature point it requires the solution of a micro-scale finite element problem (one computation of the effective flux and two sensitivity problems) followed by the assembly of the internal flux vector. For the complete solution procedure of the CTH see [18]  and [19]. For low levels of homogenization, it becomes even more expensive than the DNS problem.
Since the enriched-continuum formulation replaces the microscopic finite element problem with a set of ordinary differential equations, the computational gains are remarkably high. The internal-variable solution method is a little more expensive than the multi-field method because it requires the assembly of the internal flux vector at each time increment using the data from the previous time step. In the internal-variable method, the enrichment-variables are computed at the macroscopic Gauss quadrature points entailing higher memory requirements than the multi-field method, since the total number of integration points in the model is usually (much) larger than the total number of finite element nodes. The multi-field method is the least expensive because the finite element assembly is performed only once and it does not require the assembly of the internal flux vector at each time step. Moreover, if a reordering algorithm is used, such as Cuthill-McKee [30], the direct linear solvers can be optimized resulting in even lower computation times. The CPU time comparison has been performed for the same time integration scheme, i.e. backward-Euler, for all considered methods. The computational gains, shown in Figure 9, are expected to be even more prominent in a three-dimensional setting, consistent with the homogenization literature for micro-scale simulations [31,32] and coupled two-scale simulations [33].
To further elaborate on the relative computational costs between the multi-field and the internal-variable method, simulations were performed with an increasing number of enrichmentvariables at a homogenization level h = 2.78 provided by nel x 1 × nel x 2 = 30 × 3 elements. With the increase in the number of enrichment-variables, solved at the macro-scale, the computational cost of the multi-field method increases relative to internal-variable method, see Figure  10. It was also observed that the computational cost of the internal-variable method remained nearly constant for a changing number of enrichment-variables. This suggests that, in the internal-variable method, the major cost incurred is due to the element level calculations and the assembly of the internal flux vector F int n . Notice, however, that as discussed before, the internal-variable methods allows for more flexibility, e.g. in the choice of different time integration schemes. In this study, for N q = 6 the multi-field method still remains less expensive than the internal-variable method, however, it can be inferred from Figure 10 that it can become more expensive if N q increases due to a change in the microstructural morphology or material properties.

Conclusions
In this work, the enriched continuum formulation, which emerges from model reduction at the micro-scale, is applied to the solution of a heterogeneous diffusion boundary value problem with pronounced micro structural transient effects. Different solution methods for the enrichedcontinuum transient diffusion problem formulation are discussed. These solution methods make use of two discrete spatial discretization schemes for the enrichment-variables. The primary macroscopic field variable is always interpolated using the finite element shape functions, while the field of enrichment-variables can either be discretized using finite elements, leading to a multi-field solution method, or evaluated at the Gauss quadrature points, leading to an internal-variable solution method. To capture the transient effects more accurately, a one-step second-order accurate time integration scheme is proposed for the internal-variable method. Results for the enriched-continuum formulation with the proposed solution methods are compared with the conventional computational transient homogenization (CTH) scheme and direct numerical simulations (DNS). The proposed solution methods provide the same accuracy with respect to DNS as the expensive CTH with high gains in computational times. The CPU time and the memory requirements for the multi-field method was the lowest among the proposed solution methods. eigenvectors Φ . The expressions for the macroscopic flux and rate term in equation (9) can be converted to boundary integrals by using the divergence theorem, which in their discrete form can be written as where ∆x p = (x p − x) and I p is the column of ones of length (p × 1). By substituting the expression for j p n from (39) into (45) and rearranging the terms, the macroscopic constitutive form in equations (11) and (12) where V is the volume of the microscopic domain Ω as shown in Figure 2(b).