Macroscale modelling of multilayer diffusion: Using volume averaging to correct the boundary conditions

Applied Mathematical Modelling 47 (2017) 600–618

Contents lists available at ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

Macroscale modelling of multilayer diffusion: Using volume averaging to correct the boundary conditions E.J. Carr a,∗, I.W. Turner a,b, P. Perré c a

School of Mathematical Sciences, Queensland University of Technology (QUT), Brisbane, Australia Australian Research Council Centre of Excellence for Mathematical and Statistical Frontiers (ACEMS), Queensland University of Technology (QUT), Brisbane, Australia c LGPM, CentraleSupélec, Université Paris-Saclay, Châtenay-Malabry, France b

a r t i c l e

i n f o

Article history: Received 14 July 2016 Revised 13 January 2017 Accepted 21 March 2017 Available online 31 March 2017 Keywords: Multilayer diffusion Composite medium Macroscale Volume averaging Homogenization Boundary conditions

a b s t r a c t This paper investigates the form of the boundary conditions (BCs) used in macroscale models of PDEs with coefficients that vary over a small length-scale (microscale). Specifically, we focus on the one-dimensional multilayer diffusion problem, a simple prototype problem where an analytical solution is available. For a given microscale BC (e.g., Dirichlet, Neumann, Robin, etc.) we derive a corrected macroscale BC using the method of volume averaging. For example, our analysis confirms that a Robin BC should be applied on the macroscale if a Dirichlet BC is specified on the microscale. The macroscale field computed using the corrected BCs more accurately captures the averaged microscale field and leads to a reconstructed microscale field that is in excellent agreement with the true microscale field. While the analysis and results are presented for one-dimensional multilayer diffusion only, the methodology can be extended to and has implications on a broader class of problems. © 2017 Elsevier Inc. All rights reserved.

1. Introduction Mathematical models taking the form of partial differential equations (PDEs) involving small-scale variation in coefficients are ubiquitous in engineering and scientific disciplines. Popular applications include heat conduction in composite materials [1], groundwater flow and contaminant transport in heterogeneous aquifers [2,3] and oil recovery in subsurface reservoirs [4], where the spatial heterogeneity of the conductivity/permeability dictates the flow behaviour. For such problems, direct numerical solution of the governing equations via a discretisation of the domain yields a prohibitively large number of mesh elements and hence a discrete system which is exorbitantly expensive to solve [5]. Such numerical issues motivate the need for a macroscale approach. In this work, the microscale refers to the usual scale at which transport phenomena is described, namely the classical macroscopic or continuum scale, where for example the laws of Fourier and Darcy hold. On the other hand, the macroscale refers to a larger scale where the heterogeneous medium can be replaced by an equivalent or effective homogeneous medium [6]. We assume that the boundary conditions (BCs) at the microscale are known or can be controlled (i.e., the parameters appearing in the BCs are measurable) with our goal being to develop a macroscale model, whose solution U(x, t) provides a good approximation to the smoothed/averaged behaviour of the microscale field u(x, t) (in some sense). ∗

Corresponding author. E-mail addresses: [email protected], [email protected] (E.J. Carr), [email protected] (I.W. Turner), [email protected] (P. Perré).

http://dx.doi.org/10.1016/j.apm.2017.03.044 0307-904X/© 2017 Elsevier Inc. All rights reserved.

E.J. Carr et al. / Applied Mathematical Modelling 47 (2017) 600–618

601

Consider, for example, the diffusion equation in a heterogeneous medium:

  ∂u ∂ ∂u κ (x ) = , l0 < x < lm , ∂t ∂ x ∂x

(1)

where the diffusivity κ (x) varies over a microscale which is small relative to the length of the domain lm − l0 . Formulating a macroscale model of (1) can be achieved in several ways. The simplest approach is to assume that the diffusivity κ (x) can be replaced by an averaged or effective value κ eff , which is constant (or varies slowly with the spatial coordinate x):

  ∂U ∂ ∂U = , l0 < x < lm . κ ∂t ∂ x eff ∂ x

(2)

Note that solving macroscale equations such as (2) require significantly less computational effort compared with solving microscale equations such as (1) over the entire domain as the grid resolution is not constrained by the small-scale heterogeneous geometry. A more rigorous approach is to use upscaling methods such as homogenization [7] and volume averaging [8] to derive the macroscale equation. A recent and comprehensive comparison between both approaches is given by Davit et al. [9]. For volume averaging in one spatial dimension, we define the averaging operator:

u(x, t ) =

1 2δ



x+δ

x −δ

u (ζ , t ) d ζ ,

(3)

which replaces u(x, t) by its local average in a small interval of length 2δ around each point x. The central idea of volume averaging is to obtain an approximation to the averaged field u(x, t) (without computing the microscale field u(x, t)) by deriving a macroscale equation through direct spatial averaging of the microscale equation via application of the operator (3). Various theorems for interchanging averaging and differential operators then permit a macroscale equation to be determined, the solution of which approximates u(x, t). The idea behind homogenization is less intuitive. Given that there are two inherent length scales in the problem, the idea is to introduce a second spatial variable y = x/ε (scale of the heterogeneities) and consider the family of problems parameterised by ε . Assuming the solution uε (x, y, t) has an expansion of the form:

uε (x, y, t ) = u0 (x, y, t ) + ε u1 (x, y, t ) + ε 2 u2 (x, y, t ) + · · · ,

(4)

homogenization seeks the asymptotic solution as ε → 0 (i.e., letting the heterogeneities tend to zero). Ultimately, for the diffusion equation (1), both approaches lead to the same macroscale Eq. (2), where the the effective diffusion coefficient κ eff is given by the harmonic average of the microscale diffusivity [10]. In the method of volume averaging, the macroscale field U(x, t) is an approximation to the averaged field u(x, t) while in the homogenization method U(x, t) is equal to the leading order approximation u0 (x, t) in the expansion (4) (which can be shown to be independent of the small length scale y). In this paper, we aim to answer the following important question: given a set of BCs on the microscale, what form should the BCs on the macroscale take? In other words, what choices for the macroscale BCs give the best match between the macroscale and microscale fields? Traditionally, the microscale BCs do not play any role in the derivation of the macroscale equation in neither volume averaging nor homogenization. As pointed out by Pavliotis and Stuart [10], the microscale BCs are “irrevelant” in the homogenization procedure as they do not enter into the analysis: exactly the same macroscale equation is generated regardless of the form of the microscale BCs. For both approaches (volume averaging and homogenization), usually the macroscale BCs imposed are of the same type as the microscale BC. For example, if a Dirichlet BC is applied on the microscale field then the same Dirichlet condition is imposed on the macroscale field [5,9–14]. If the microscale flux at the boundary is specified and assigned a constant value then the macroscale flux at the same boundary is assumed equal to that constant [5,11,12,15]1 . In the context of homogenization by asymptotic expansion, for certain differential operators and sets of BCs it is possible to rigorously prove convergence of the sequence of microscale solutions uε (x, y, t) towards u0 (x, t) in the limit ε → 0 [7,10]. This, of course, assumes that the ratio of the microscale to macroscale length scales is infinitesimally small, which is often a flawed assumption in many real-world problems. As a means of validating the macroscale model and corresponding BCs, in this work we use the solution of the microscale model as the reference solution (referred to by Davit et al. [9] as validation by direct numerical simulation)2 . Two examples that raise concerns about the common choices for the macroscale BCs are given in Fig. 1. In Example 1, three snapshots of the microscale field u(x, t) are given with a Dirichlet BC imposed at the left boundary. Overlayed is an averaged (macroscale) field computed by averaging u(x, t) via the operator (3). Observe that u(x, t) satisfies a similar evolution law to the microscale field, however, it clearly does not satisfy the Dirichlet BC specified on the microscale at the left boundary3 . In Example 2, the Dirichlet condition at the right boundary is modified to a Robin condition. Applying the same Robin 1

For example, for the diffusion equations (1) and (2) the microscale BCs u(x, t ) = a and κ1 ∂∂ ux = q at x = 0 traditionally lead to macroscale BCs of

U (x, t ) = a and κeff ∂∂Ux = q being proposed at x = 0, respectively. 2 For example, for the diffusion problem, to validate the macroscale model, the macroscale field U(x, t) is compared to the microscale field u(x, t). 3 Unless δ = 0, the average of the microscale field over the interval [0, δ ] is less than 1. Also, for the macroscale field to pass through the centre of the oscillations, u(0, t) cannot be equal to 1 for all t > 0.

602

E.J. Carr et al. / Applied Mathematical Modelling 47 (2017) 600–618

Fig. 1. Solution u(x, t) of the multilayer diffusion problem (5) with l0 = 0, lm = 1, m = 20 layers, n = 10 periods (periodicity p = 2), period diffusivities [κ1 , κ2 ] = [0.1, 1], initial condition ui (x, 0 ) = 0 (i = 1, . . . , m), left BC u(0, t ) = 1 and constant layer width d = 0.05 such that the interfaces are located at x = li = id (i = 1, . . . , m − 1). The right BC is ∂∂ ux (1, t ) = 0 in Example 1 and u(1, t ) + ∂∂ ux (1, t ) = 1 in Example 2. For Example 1, the averaged microscale field

u(x, t) is computed at discrete points: x = 0, d2 , 32d , . . . , (m − 12 )d, 1 using the operator (3) with δ = d/2. For Example 2, the macroscale solution U(x, t) is shown, computed using the macroscale Eq. (12a) subject to the same initial and boundary conditions imposed on the microscale: U (x, 0 ) = 0, U (0, t ) = 1 and U (1, t ) + ∂∂Ux (1, t ) = 1.

condition at the macroscopic level leads to a poor match between the macroscale and microscale fields (U(x, t) and u(x, t)). What should the corresponding macroscale BCs be in this case? Recently, Chen [16] employed centre-manifold theory and modern dynamical systems theory to derive a macroscale model of a one-dimensional lattice diffusion system, where the diffusivities are constant between lattice points and repeat periodically throughout the domain. A key result of the work is the derivation of non-traditional macroscale BCs for different choices of microscale BCs. For example, a macroscale BC of Robin type is derived for a microscale BC of Dirichlet type and is shown to improve the accuracy of the macroscale model. The methodology is also applicable to other types of BCs and other two-scale problems, such as a wave equation [16, Chapter 5], a nonlinear heat exchanger [17] and a ‘two-strand’ lattice diffusion system [18]. The task of determining macroscale BCs is closely related to but distinct from the well-established concepts of boundary layers and edge effects in the homogenization literature (see, e.g., [19–22]). In such approaches, boundary layers are introduced into the asymptotic expansion (4), which satisfy elliptic problems defined on a semi-infinite strip and exhibit exponential decay in the interior of the domain to a non-zero limit known as the boundary layer tail [19]. Primarily, the idea is to improve the accuracy of the asymptotic expansion (4) taken beyond leading order [20–22]. An exception is the work by Allaire and Amar [19] on the Dirichlet problem, who demonstrate how these boundary layer tails can be incorporated into the homogenized (macroscale) equation by replacing the Dirichlet BC with a suitably-defined Robin condition. In the current paper, we corroborate and extend the work of Chen [16] and Allaire and Amar [19] using a volume averaging approach via the usual ‘average plus perturbation’ decomposition of the microscale field [9]. By accounting for the perturbation term, we derive a set of corrected macroscale BCs for, not only Dirichlet microscale BCs, but for several choices of microscale BCs. Apart from some special cases, the corrected BCs differ from both the original microscale BCs imposed and classical choices for the macroscale BCs. Specifically, our analysis confirms that a macroscale BC of Robin type emerges from a microscale BC of Dirichlet type. Additionally, we generalise the corrected BCs to microscale BCs of Neumann, Flux and Robin type obtaining explicit expressions for the (corrected) macroscale BCs in each case. While the analysis presented applies to the specific problem of one-dimensional multilayer diffusion we believe that it is best to consider such a simple prototype problem to convey the ideas. Nonetheless, multilayer diffusion is an important research problem in its own right as it finds application to a wide array of industrial and biological applications (see, e.g., Hickson et al. [23]). This research presented in this paper follows the recent work of two of the current authors [24], where a semi-analytical solution for one-dimensional multilayer diffusion was derived and applied to an application involving macroscale modelling with standard macroscale BCs. The remaining sections of this paper are outlined as follows. In the next section, we present the microscale model and then in Section 3 describe the macroscale model including the forms of the macroscale equation, effective diffusivity, corrected macroscale BCs and the reconstructed microscale field. The majority of the derivation of the results of Section 3 is included as an appendix (Appendix A). In Section 5, several test cases are presented that demonstrate that the corrected BCs improve the accuracy of both the macroscale field and reconstructed microscale field. 2. Microscale model Consider diffusion in a one-dimensional composite medium on the interval l0 < x < lm consisting of m layers. Let κ i denote the constant diffusivity in the ith layer and denote by x = li the location of the interface between layers i and i + 1.

E.J. Carr et al. / Applied Mathematical Modelling 47 (2017) 600–618

603

Fig. 2. Example periodic composite medium with m = 10 layers, n = 5 periods (periodicity p = 2) and period diffusivities [κ 1 , κ 2 ]. The macroscale approach replaces the heterogeneous medium with a fictitious homogeneous medium.

We assume the medium has a periodic structure with n periods/cells, each consisting of p layers so that the total number of layers m is equal to np. The period has length L = (lm − l0 )/n = l p − l0 and the diffusivities in the period, denoted by the array [κ1 , . . . , κ p ], have periodicity p, so that: κk+ j p = κk for k = 1, . . . , p and j = 1, . . . , n − 1. The configuration considered is given in Fig. 2 for p = 2 and m = 10. Note that the non-periodic problem is recovered if p = m. The microscale model consists of the diffusion equation:

  ∂ ui ∂ ∂u = κi i , ∂t ∂x ∂x

(5a)

in each layer (i = 1, . . . , m), subject to the initial and external boundary conditions:

ui (x, 0 ) = uint ,

i = 1, . . . , m,

∂ u1 ( l , t ) = cL , ∂x 0 ∂ um aR um ( lm , t ) + bR ( l , t ) = cR , ∂x m aL u1 ( l0 , t ) − bL

(5b) (5c) (5d)

and internal BCs at the interfaces:

ui+1 (li , t ) = ui (li , t ),

κi+1

∂ ui+1 ∂u (l , t ) = κi i (li , t ), ∂x i ∂x

(5e) (5f)

for i = 1, . . . , m − 1, where li−1 < x < li , t > 0 and uint is a constant. The scalars aL , bL , aR and bR are all assumed to be non-negative constants satisfying aL + bL = 0 and aR + bR = 0. Henceforth, we will refer to (5) as the microscale model. 3. Macroscale model 3.1. Formulation A key technique employed in the method of volume averaging (see, e.g., Davit et al. [9]) is to decompose the microscale solution ui (x, t) into two components: an average component u(x, t) and a perturbation component  ui (x, t ), such that

ui (x, t ) = u(x, t ) +  ui (x, t ),

(6)

for i = 1, . . . , m. Averaging the microscale Eq. (5a) over the interval [x − δ, x + δ ] using the operator (3) and introducing the above decomposition, one can show that under various assumptions (see Appendix A) the average component approximately satisfies the macroscale equation:

  ∂ u ∂ ∂ u  , l0 < x < lm , κ ∂t ∂ x eff ∂ x

(7)

and that the perturbation term is approximately given by:

 ui (x, t )  ψi (x )

∂ u (x, t ), li−1  x  li , ∂x

(8)

for i = 1, . . . , m. For k = 1, . . . , p, the functions ψ k (x) satisfy the periodic cell problem:

 d    κ ψ ( x ) + 1 = 0, dx k k

(9a)

604

E.J. Carr et al. / Applied Mathematical Modelling 47 (2017) 600–618

ψ1 ( l0 ) = ψ p ( l p ),     κ1 ψ1 (l0 ) + 1 = κ p ψ p (l p ) + 1 ,

(9b)

ψk (lk ) = ψk+1 (lk ), k = 1, . . . , p − 1,     κk ψk (lk ) + 1 = κk+1 ψk+1 (lk ) + 1 , k = 1, . . . , p − 1,

(9d)

(9c)

p  1  lk ψk (x ) dx = 0, l p − l0 lk−1

(9e) (9f)

k=1

where p is the periodicity defined in the microscale model (5). The effective diffusivity appearing in the macroscale (averaged) Eq. (7) is found to be (see Appendix A):

κeff =

p    1  lk κk ψk (x ) + 1 dx. l p − l0 lk−1

(10)

k=1

Solving (9) and substituting the solution into the integral expression for the effective diffusivity (10) yields the classical result for one-dimensional layered media, namely, the harmonic average of the individual diffusivities:



κeff =

p 1  l p − l0



lk − lk−1

k=1

 −1

.

κk

We remark that the solutions for the periodic cell problems are given by ψk (x ) = αk + βk x, where βk = (κeff /κk ) − 1. Note further that these functions are periodic (Assumption 5 of Appendix A) so ψk+ j p (x + jL ) = ψk (x ) for k = 1, . . . , p and j = 1, . . . , n − 1, where lk−1  x  lk and L = l p − l0 is the length of the period. Note that similar results, to those presented in this section, apply in the case of a more general diffusivity function κ (x) [9,14]. 3.2. Corrected macroscale boundary conditions To obtain the correct form of the BCs to impose on the macroscale Eq. (7), the “average plus perturbation” decomposition (6) is introduced into the microscale BCs. For example, introducing (6) and (8) with i = 1 into the left microscale BC (5c) gives:



aL u(l0 , t ) + aL u˜1 (l0 , t ) − bL

 ∂ u˜1 ∂u (l , t ) + ( l , t ) = cL . ∂x 0 ∂x 0

(11)

Using the form of the perturbation component (8) and Assumption 3 of Appendix A, which is used to derive the macroscale Eq. (7), we have:

∂ u1 ∂ u  ψ1 (x ) . ∂x ∂x Using this result in Eq. (11), yields:

aL

    ∂ u ∂ u ∂ u ( l0 , t ) − bL (l0 , t ) + ψ1 (l0 ) ( l0 , t )  cL , u(l0 , t ) + ψ1 (l0 ) ∂x ∂x ∂x  ∂ u aL u(l0 , t ) − bL (1 + ψ1 (l0 )) − aL ψ1 (l0 ) ( l , t )  cL . ∂x 0

The same analysis applied to the right microscale BC (5d) produces:



aR u(l0 , t ) + bR (1 + ψ p (lm )) − aR ψ p (lm )

∂ u ( l , t )  cR , ∂x m

where we note that ψm (lm ) = ψ p+(n−1 ) p (l p + (n − 1 )L ) = ψ p (l p ) due to periodicity. With U ࣃ u, we obtain the macroscale model, which is summarised as follows:

  ∂U ∂ ∂U = , κ ∂t ∂ x eff ∂ x

U (x, 0 ) = uint , aL U ( l0 , t ) −  bL

(12a) (12b)

∂U ( l , t ) = cL , ∂x 0

(12c)

E.J. Carr et al. / Applied Mathematical Modelling 47 (2017) 600–618

aR U ( lm , t ) +  bR

∂U ( l , t ) = cR , ∂x m

605

(12d)

where l0 < x < lm , t > 0, and

 bL = bL (1 + ψ1 (l0 )) − aL ψ1 (l0 ),

(12e)

 bR = bR (1 + ψ p (l p )) + aR ψ p (l p ).

(12f)

Except for some special cases, the corrected macroscale BCs (Eqs. 12c and 12d) are not the same as the microscale BCs (Eqs. 5c and 5d): (i) Dirichlet condition For a Dirichlet microscale BC (aL = 1 and bL = 0):

u 1 ( l0 , t ) = cL ,

(13)

the corrected macroscale BC is of Robin type, and takes the form:

U ( l0 , t ) + ψ1 ( l0 )

∂U ( l , t ) = cL . ∂x 0

(14)

This BC is the same as the microscale BC for the special case ψ1 (l0 ) = 0. (ii) Neumann condition For a Neumann microscale BC (aL = 0 and bL = −1):

∂ u1 ( l , t ) = cL , ∂x 0 the corrected macroscale BC is also a Neumann condition, and takes the form:

(1 + ψ1 (l0 ))

∂U ( l , t ) = cL , ∂x 0

or equivalently,

κeff ∂ U ( l , t ) = cL , κ1 ∂ x 0

(15)

since ψ1 (x ) = (κeff /κ1 ) − 1. This BC is the same as the microscale BC if either κeff = κ1 or cL = 0 (zero Neumann condition). (iii) Flux condition Consider a microscale BC that prescribes a microscale flux at x = l0 equal to cL in the negative x-direction (aL = 0, bL = −κ1 and outward unit normal n = −1):

−κ1

∂ u1 ∂ u1 · n := κ1 ( l , t ) = cL . ∂x ∂x 0

The corrected macroscale BC takes the form:

κ1 (1 + ψ1 (l0 ))

∂U ( l , t ) = cL , ∂x 0

or, equivalently:

κeff

∂U ( l , t ) = cL , ∂x 0

which states that the macroscale flux at x = l0 in the negative x-direction must also be equal to cL . (iv) Newton condition Consider a Newton microscale BC, where the microscale flux is proportional to the difference between the boundary value u(l0 , t) and a surrounding value u∞ (aL = σ , bL = κ1 and cL = σ u∞ ):

κ1

∂ u1 ( l , t ) = σ ( u1 ( l0 , t ) − u∞ ). ∂x 0

(16)

The corrected macroscale BC takes the form:



κ1 (1 + ψ1 (l0 )) − σ ψ1 (l0 )

∂U (l , t ) = σ (U (l0 , t ) − u∞ ), ∂x 0

or, equivalently:

[κeff − σ ψ1 (l0 )]

∂U (l , t ) = σ (U (l0 , t ) − u∞ ). ∂x 0

Note that the macroscale flux is not proportional to the difference between the boundary value U(l0 , t) and the surrounding value u∞ , as is the case for the microscale BC (16).

606

E.J. Carr et al. / Applied Mathematical Modelling 47 (2017) 600–618

A natural question to ask is what happens for a homogeneous medium? If κk = κ (k = 1, . . . , p) then κeff = κ and the solution of the cell problem (9) is ψk (x ) = 0 (k = 1, . . . , p). Therefore, for a homogeneous medium, the corrected macroscale BCs (Eqs. 12c and 12d) reduce to the same BCs imposed on the microscale (Eqs. 5c and 5d), namely:

aL U ( l0 , t ) − bL

∂U ∂U ( l , t ) = cL , aR U ( lm , t ) + bR ( l , t ) = cR . ∂x 0 ∂x m

3.3. Reconstructed microscale field The “average plus perturbation” decomposition (6) can be used to compute a reconstructed microscale field that provides a good approximation to the true microscale field u(x, t) [9]. Substituting U ࣃ u and ui  ui into (6), and using (8) yields:

ui (x, t ) = U (x, t ) + ψi (x )

∂U (x, t ), ∂x

(17)

where li−1 < x < li for i = 1, . . . , m. In words, the reconstructed microscale field is computed by adding to the macroscale field a microscale perturbation, which is found by replicating the solution of the cell problem (9) periodically throughout the domain and multiplying by the derivative of the macroscale field U(x, t) (12) at each point x. 3.4. Layers of equal width Consider the special case of each layer having width d := (lm − l0 )/m so that the interfaces li = l0 + id for i = 0, . . . , m. For periodicity p = 2, the corrected macroscale BCs (12c) and (12d) simplify to:



 κeff d κ1 − κ2 ∂ U aL U ( l0 , t ) − bL − aL ( l0 , t ) = cL , κ1 2 κ1 + κ2 ∂ x   κ d κ1 − κ2 ∂ U aRU (lm , t ) + bR eff + aR ( lm , t ) = cR . κ2 2 κ1 + κ2 ∂ x

(18a) (18b)

For Dirichlet microscale BCs (aL = aR = 1 and bL = bR = 0), these BCs reduce further to:

U ( l0 , t ) −

d κ2 − κ1 ∂ U ( l0 , t ) = cL , 2 κ1 + κ2 ∂ x

(19a)

U ( lm , t ) +

d κ1 − κ2 ∂ U ( lm , t ) = cR , 2 κ1 + κ2 ∂ x

(19b)

which are the same as the microscale BCs (u(l0 , t ) = cL and u(lm , t ) = cR ) if κ1 = κ2 (homogeneous medium) or in the limit as d → 0 (m → ∞). In general, for any periodicity p, the corrected macroscale BCs (12c) and (12d) can be expressed as:



 (2k − p − 1 )κk−1 ∂ U ( l , t ) = cL , p ∂x 0 κ −1 k=1 k   p −1 κ d k=1 (2k − p − 1 )κk ∂U aRU (lm , t ) + bR eff + aR ( l , t ) = cR . p −1 κp 2 ∂x m κ k=1 k aL U ( l0 , t ) − bL

κeff d − aL κ1 2

p

k=1

(20a) (20b)

For Dirichlet BCs (aL = aR = 1 and bL = bR = 0) the above BCs reduce to those derived by Chen [16] for the onedimensional lattice diffusion system (in their approach the lattice/grid points are located at the interfaces between layers so d also denotes the constant spacing between adjacent lattice points). Recall that the main idea of the homogenization approach is to let the heterogeneities tend to zero. In this sense, with d → 0 (or m → ∞), the corrected BCs (20) simplify to:

κeff ∂ U ( l , t ) = cL , κ1 ∂ x 0 κ ∂U aRU (lm , t ) + bR eff ( l , t ) = cR . κp ∂ x m aL U ( l0 , t ) − bL

(21a) (21b)

We will refer to the above BCs, which recover the usual macroscale BCs imposed when the microscale BC is of Dirichlet or Flux type, as the homogenization limit BCs due to their method of construction. For the simplest case of periodicity p = 2, Table 1 lists various types of microscale BCs at the left-boundary (x = l0 ) together with the corresponding homogenization limit and corrected macroscale BCs. The two choices of macroscale BCs differ when a Dirichlet or Newton microscale BC is imposed. The homogenization limit BCs (21) are also recovered in the following two special cases. For p even, if κ j = κ p− j+1 for p all j = 1, . . . , p/2, then (2k − p − 1 )κk−1 = 0 (in Eq. (20)) due to cancellation of the k = j and k = p − j + 1 terms in k=1 p the summation for all j = 1, . . . , p/2. Similarly, for p odd, if κ j = κ p− j+1 for all j = 1, . . . , ( p − 1 )/2, then ( 2k − p − k=1 1 )κk−1 = 0 (in Eq. (20)) due to the zero k = ( p + 1 )/2 term as well as cancellation of the k = j and k = p − j + 1 terms in the summation for all j = 1, . . . , ( p − 1 )/2.

E.J. Carr et al. / Applied Mathematical Modelling 47 (2017) 600–618

607

Table 1 Form of the homogenization limit and corrected macroscale BCs corresponding to various choices of the left microscale BC (5c) for the simple case of periodicity p = 2 and constant layer width d.

Dirichlet

Microscale BC

u1 ( l0 , t ) = cL

Macroscale BC (homogenization limit)

U ( l0 , t ) = cL

Macroscale BC (corrected)

U ( l0 , t ) −

Microscale BC Neumann

Macroscale BC (homogenization limit) Macroscale BC (corrected) Microscale BC

Flux

Macroscale BC (homogenization limit) Macroscale BC (corrected) Microscale BC

Newton

Macroscale BC (homogenization limit) Macroscale BC (corrected)

d κ2 − κ1 ∂ U ( l0 , t ) = cL 2 κ1 + κ2 ∂ x

∂ u1 ( l , t ) = cL ∂x 0 κeff ∂ U ( l , t ) = cL κ1 ∂ x 0 κeff ∂ U ( l , t ) = cL κ1 ∂ x 0 ∂u κ1 1 ( l 0 , t ) = c L ∂x ∂U κeff ( l , t ) = cL ∂x 0 ∂U κeff ( l , t ) = cL ∂x 0 ∂ u1 κ1 ( l , t ) = σ ( u1 ( l0 , t ) − u∞ ) ∂x 0 ∂U κeff (l0 , t ) = σ (U (l0 , t ) − u∞ )  ∂x  d κ1 − κ2 ∂ U κeff − σ (l0 , t ) = σ (U (l0 , t ) − u∞ ) 2 κ1 + κ2 ∂ x

4. Computation of macroscale field 4.1. Analytical solution The classical analytical solution of the macroscale model (12), is well-known [25,26], and given by:

U (x, t ) = w(x ) +

∞ 

cn e−t λn κeff φn (x ),

(22)

n=1

where w(x) is the linear steady state solution satisfying aL w(l0 ) −  bL w (l0 ) = cL and aR w(lm ) +  bR w (lm ) = cR and the coefficients in the expansion are defined as

cn =

uint − w(x ), φn (x ) ,  φn ( x ) , φn ( x ) 

where  f (x ), g(x ) =

(23)

 lm

f (x )g(x ) dx and φm (x ), φn (x ) = 0 if m = n. The eigenvalues λn (numbered and ordered such that l0 λ1 < λ2 < · · · ) and corresponding eigenfunctions φ n (x) satisfy the Sturm–Liouville problem

Lφn = λn φn ,

L=−

∂2 , ∂ x2

a L φn ( l 0 ) −  bL φn (l0 ) = 0,

(24a) a R φn ( l m ) +  bR φn (lm ) = 0.

(24b)

4.2. Negative eigenvalues The constants  bR and  bL (defined in Eqs. (12e) and (12f)) appearing in the corrected BCs (Eqs. (12c) and (12d)) are not necessarily non-negative, which means that the eigenvalues of the differential operator L (24) are not guaranteed to be positive [25]. For example, if we have a Robin BC at x = l0 and a Neumann BC at x = lm (aL > 0,  bL = 0, aR = 0 and  bR = 0) there is one negative eigenvalue if  bL < 0 [25]. Note that negative eigenvalues are particularly worrisome since the factor e−t λκeff in the solution expansion (22) grows exponentially over time. To explore further, consider the special case of the corrected BCs (18) with aL = 1, bL = 0, aR = 0, cL = 1, bR = 1, cR = 0:

∂U U ( l0 , t ) −  bL ( l , t ) = 1, ∂x 0

∂U ( l , t ) = 0, ∂x m

(25)

where  bL = d (κ2 − κ1 )/(2(κ1 + κ2 )). For these BCs, if κ 2 > κ 1 (referred to as Case 1) all the eigenvalues are positive while if κ 2 < κ 1 (referred to as Case 2) there is one negative eigenvalue and all the rest are positive. These two cases are demonstrated in Fig. 3 for an example involving m = 20 layers. Observing the evolution of the averaged field u(x, t) for Case 2 (Fig. 3b), it is clear that it is not possible to capture this behaviour using the macroscale Eq. (12a). Since ∂∂Ut > 0 and the macroscale Eq. (12a) holds, it must be true that ∂∂ xU2 > 0 for 0 < x < 1. However, the macroscale (averaged) field in Case 2 2

608

E.J. Carr et al. / Applied Mathematical Modelling 47 (2017) 600–618

Fig. 3. Solution u(x, t) of the multilayer diffusion problem (5) with the same conditions as Example 1 of Fig. 1. The perturbation field is defined as  u(x, t ) = u(x, t ) − u(x, t ). A vertical line is drawn at x = l1 to depict the first layer.

(Fig. 3b) possesses an inflection point near the left boundary thus contradicting the requirement of ∂∂ xU2 > 0 for 0 < x < 1. Hence, negative eigenvalues seem to indicate that the macroscale Eq. (12a) is no longer suitable. Indeed, as can be seen in Fig. 3b in the first layer, Case 2 violates Assumption 3 ( u1 exhibits nonlinear behaviour in the first layer, which is impossible ∂ u ∂ u if  u1 = ψ1 ∂x  with ψ 1 linear and ∂x  constant) and Assumption 6 ( u  0 clearly doesn’t hold), both of which are used to derive the macroscale Eq. (12a) in Appendix A. We will not pursue modification of the macroscale Eq. (12a) in the present work. Instead we study the structure of the macroscale field (22) with the negative eigenvalues removed. Suppose the eigenvalue problem (24) possesses r negative eigenvalues (r = 0, 1 or 2), that is, λn < 0 for n ≤ r and λn > 0 for n > r. In this case to compute the macroscale field, we exclude the negative eigenvalues in the solution expansion (22) by defining the sum over the non-negative spectrum of L only: 2

U (x, t ) = w(x ) +

∞ 

cn e−t λn κeff φn (x ).

(26)

n=r+1

Although the above strategy is rather crude, it is necessary to produce physically reasonable solutions. We note that this strategy is similar to one employed by Chen et al. [18], where exponentially growing modes were excluded in the left boundary layer when deriving macroscale BCs. The macroscale field (26) satisfies the corrected BCs (Eqs. (12c) and (12d) but not the initial condition (12b). In fact, one can show that (26) is nothing more than the solution of the macroscale model (12) with the following modified initial condition:

U (x, 0 ) = uint −

r  n=1

c n φn ( x ) = w ( x ) +

∞ 

c n φn ( x ) ,

n=r+1

which is the orthogonal projection of the macroscale initial condition (12b) onto the affine space

A = w(x ) + span{φr+1 (x ), φr+2 (x ), . . .}.

(27)

E.J. Carr et al. / Applied Mathematical Modelling 47 (2017) 600–618

609

The solution of the macroscale Eq. (12a) subject to (27), (12c) and (12d) is

U (x, t ) = w(x ) +

∞ 

 cn e−t λn κeff φn (x ),

(28)

n=1

where the coefficients in the expansion are given by

   uint − w(x ) − rk=1 ck φk (x ), φn (x ) 0 n≤r  cn = = cn n > r.  φn ( x ) , φn ( x ) 

Hence, it is clear that (28) and (26) are equal. Example. Consider the problem depicted in Fig. 3b with an arbitrary number of layers m. For this problem, there is one negative eigenvalue (r = 1), which can be identified as λ1 = −μ2 , where μ is the unique positive root of the equation 9 tanh(μ ) = −1/(μ bL ), where  bL = d (κ2 − κ1 )/(2(κ1 + κ2 )) = − 22 d (κ1 = 1 and κ2 = 0.1) and d = 1/m is the constant layer width. The corresponding eigenfunction is φ1 (x ) = cosh(μ(1 − x )). With these results, we have

c1 =

uint − w(x ), φ1 (x ) 4(uint − cL )(e−3μ − e−μ ) = ,  φ1 ( x ) , φ1 ( x )  1 − e−4μ + 4μe−2μ

and hence the modified initial condition (27) can be found to be

U (x, 0 ) = uint − c1 φ1 (x ) = uint +

2(cL − uint )(1 − e−2μ )(eμ(x−2) + e−μx ) . 1 − e−4μ + 4μe−2μ

The value of μ increases with increasing m with the asymptotic relation μ ∼ −1/ bL = m. This leads to the following asymptotic relation:

22 9d

accurate even for small values of

U (x, 0 ) ∼ uint + 2(cL − uint )e−μx = uint + 2(cL − uint )e− 9d x , 22

for the modified initial condition (27) as the layer width d tends to zero. This analysis indicates that there is a boundary layer of thickness O (d ) at the left boundary: U(x, 0) is essentially indistinguishable from uint everywhere except near x = 0. Therefore, we expect that at short times the accuracy of the macroscale field (26) (and hence reconstructed microscale field (17)) will be poor. 5. Results In this section, we assess the accuracy of the macroscale model (12) using the microscale model (5) as the benchmark. Four test cases are considered, each involving a composite medium on the interval 0  x  1 (l0 = 0 and lm = 1) with an even number of layers m, constant layer width di = 1/m (i = 1, . . . , m) and periodicity p = 2. In addition to visually comparing the solution fields, we also record the following root mean square errors:



er rha :=

 er rca

:=

m 2 1  u(xi , t ) − U (xi , t ) m

(homogenization limit BCs),

(29a)

m 2 1  u(xi , t ) − U (xi , t ) m

(corrected BCs),

(29b)

k=1

k=1

  Nx m   1   2 m er rh :=  u(xi, j , t ) − U (xi, j , t ) (homogenization limit BCs), mNx

  Nx m   1   2 m er rc :=  u(xi, j , t ) − U (xi, j , t ) (corrected BCs), mNx

(29d)

i=1 j=1

  Nx m   1   2 m er rr :=  u(xi, j , t ) − u(xi, j , t ) (corrected BCs). mNx

(29c)

i=1 j=1

(29e)

i=1 j=1

The first and second error metrics measure the accuracy of the macroscale field U(x, t) using the averaged microscale field

u(x, t) as the reference solution, while the third and fourth error metrics instead use the true microscale field u(x, t) as the reference solution. The fifth error metric measures the accuracy of the reconstructed microscale field u(x, t ) using the true microscale field u(x, t) as the reference solution. In the first two error metrics, the sum is taken over the number of layers

610

E.J. Carr et al. / Applied Mathematical Modelling 47 (2017) 600–618

m, where xi = (li−1 + li )/2 is the midpoint of the ith layer, and u(xi , t) is the average value of the true microscale solution u(x, t) over the ith layer. In the third, fourth and fifth error metrics, Nx is the number of uniformly-spaced grid points in each layer, at which the fields are evaluated, such that xi, j = li−1 + j (li+1 − li )/(Nx − 1 ) (i = 1, . . . , m and j = 1, . . . , Nx ). Note that for the macroscale field U(x, t) (26) the averaged microscale field u(x, t) is used as the reference solution while for the reconstructed microscale field u(x, t ) (17) the true microscale field u(x, t) (5) is used. The four individual test cases are summarised below: (a) Case A imposes Dirichlet and Neumann BCs at the left and right boundaries, respectively:

∂ um ( 1, t ) = 0, ∂x

u1 ( 0, t ) = 1,

(30)

and an initial condition of ui (x, 0 ) = 0 (i = 1, . . . , m). The period diffusivities are taken to be [κ1 , κ2 ] = [0.1, 1.0] giving an an effective diffusivity of κeff = 2/11 ≈ 0.181818. The corresponding corrected macroscale BCs (18) take the form:

∂U U ( 0, t ) −  bL ( 0, t ) = 1, ∂x

∂U ( 1, t ) = 0, ∂x

(31)

where  bL = −ψ1 (0 ) = d (κ2 − κ1 )/(2(κ1 + κ2 )) = 9/440 ≈ 0.0204545. The corresponding homogenization limit BCs (21) are the same as those applied on the microscale (30):

∂U ( 1, t ) = 0. ∂x

U ( 0, t ) = 1,

(32)

Note that this is the same problem to the one shown in Fig. 3a. (b) Case B is the same as Case A except the order of the diffusivities are reversed: [κ1 , κ2 ] = [1.0, 0.1]. This produces a single change from Case A: the corrected BCs now take the form of (31) with  bL = −ψ1 (0 ) = d (κ2 − κ1 )/(2(κ1 + κ2 )) = −9/440 ≈ −0.0204545. Note that this is the same problem to the one shown in Fig. 3b. (c) Case C is the same as Case A with the exception that a Newton BC is applied at x = 1:

−κm

u1 ( 0, t ) = 1,

∂ um ( 1, t ) = σ ( um ( 1, t ) − u∞ ), ∂x

(33)

where σ = 2.0 and u∞ = 1.0. The corrected macroscale BCs (18) are given by:

∂U U ( 0, t ) −  bL ( 0, t ) = 1, ∂x

− κeff

∂U (1, t ) = σ (U (1, t ) − u∞ ), ∂x

(34)

where  bL is as defined in Case A and  κeff = κeff + σ ψ p (0.1 ) = 31/220 ≈ 0.140909. The corresponding homogenization limit BCs (21) take the form:

U ( 0, t ) = 1,

−κeff

∂U (1, t ) = σ (U (1, t ) − u∞ ). ∂x

(35)

(d) Case D imposes Dirichlet BCs at both boundaries:

u1 ( 0, t ) = 0.5,

−κm

∂ um (1, t ) = −0.4, ∂x

(36)

and an initial condition of ui (x, 0 ) = 1 (i = 1, . . . , m). The corrected macroscale BCs (18) are given by:

∂U U ( 0, t ) −  bL ( 0, t ) = 0.5, ∂x

−κeff

∂U (1, t ) = −0.4, ∂x

(37)

where  bL is as defined in Case A. The corresponding homogenization limit BCs (21) take the form:

U ( 0, t ) = 0.5,

−κeff

∂U (1, t ) = −0.4. ∂x

(38)

Computed macroscale and reconstructed microscale fields for Cases A–D with m = 20 layers are shown in Figs. 4–8 at selected points in time (the first 200 terms are taken in the solution expansion (26)). The microscale field (solution of the microscale model (5)) is computed using the semi-analytical method proposed by Carr and co-authors [24,27]. Two observations are immediately obvious from Figs. 4–8: (i) the macroscale field with corrected BCs passes through the centre of the microscale field oscillations, accurately interpolating the averaged microscale field u(x, t) (Figs. 4–7). (ii) the reconstructed microscale field computed using the corrected BCs is in excellent agreement with the true microscale field (Fig. 8). Interestingly, the macroscale field with corrected BCs tends to interpolate the microscale field at the centre of each layer (x = 0.025, 0.075, . . . , 0.975) while the homogenization limit BCs produce a macroscale field that favours interpolation at the boundaries and at every second interface (x = 0, 0.1, . . . , 1.0). As observed for Case A in Fig. 4, the corrected macroscale BC of Robin type at x = 0 effectively softens the Dirichlet microscale BC at x = 0, providing a smooth transition from the initial

E.J. Carr et al. / Applied Mathematical Modelling 47 (2017) 600–618

611

Fig. 4. Macroscale fields U(x, t) at t = 0.1tc , tc , 0.2, 1.0 for Case A with either homogenization limit or corrected BCs, where tc = d12 /κ1 is the characteristic time of diffusion in the first layer.

Fig. 5. Macroscale fields U(x, t) at t = 0.1tc , tc , 0.2, 1.0 for Case B with either homogenization limit or corrected BCs, where tc = d12 /κ1 is the characteristic time of diffusion in the first layer.

Fig. 6. Macroscale fields U(x, t) at t = 0.1tc , tc , 0.2, 1.0 for Case C with either homogenization limit or corrected BCs, where tc = d12 /κ1 is the characteristic time of diffusion in the first layer.

612

E.J. Carr et al. / Applied Mathematical Modelling 47 (2017) 600–618

Fig. 7. Macroscale fields U(x, t) at t = tc , 0.2, 1.0 for Case D with either homogenization limit or corrected BCs, where tc = d12 /κ1 is the characteristic time of diffusion in the first layer.

Table 2 Relative errors (29) for Cases A–D , where tc = d12 /κ1 is the characteristic time of diffusion in the first layer. The error metrics are defined in Eq. (29) with Nx = 5 used to produce the results shown. er rha

Case

t

er rha

er rca

er rrm

A

0.1tc tc 0.2 1.0

1.24e−02 4.00e−02 2.80e−02 2.03e−02

3.56e−02 7.41e−03 1.49e−03 5.56e−04

1.21e−02 4.06e−03 7.82e−04 2.77e−04

0.3 5.4 18.8 36.4

B

0.1tc tc 0.2 1.0

7.79e−02 1.05e−01 3.12e−02 2.16e−02

6.78e−02 1.24e−02 1.47e−03 5.68e−04

7.93e−02 7.89e−03 7.26e−04 2.67e−04

1.1 8.4 21.3 38.0

C

0.1tc tc 0.2 1.0

1.31e−02 4.22e−02 3.35e−02 1.10e−02

3.57e−02 9.41e−03 2.33e−03 6.04e−04

1.37e−02 5.59e−03 1.18e−03 3.16e−04

0.4 4.5 14.4 18.2

D

tc 0.2 1.0

2.01e−02 1.40e−02 7.83e−03

4.19e−03 1.49e−03 9.52e−04

2.33e−03 7.86e−04 4.66e−04

4.8 9.4 8.2

er rca

value of u = 0 to u = 1 instead of the instantaneous change that occurs with the Dirichlet condition. Without this softening it is impossible for the macroscale field to pass through the centre of the microscale field oscillations. Imposing the corrected BCs (18) on the macroscale field is more accurate than imposing the homogenization limit BCs (21). This can be seen in the figures but is also reflected in the tabulated relative errors (Table 2). In the table, error ratios are also given, which provide a measure of the reduction in error. For example, for Case A at t = 0.2 imposing the corrected BCs instead of the homogenization limit BCs reduces the error of the macroscale field by a factor of 18.8. Overall, the corrected BCs offer a reduction in the error by a factor of 8.1, on average across all test cases, for the macroscale field, relative to the homogenization limit BCs. For each test case it is interesting to study the eigenvalue problem (24) corresponding to the corrected BCs. For Cases A, C and D, all the eigenvalues are positive while for Case B there is one negative eigenvalue. By excluding this negative eigenvalue in the solution expansion (as outlined in Section 4.2), the resulting macroscale field (26) performs remarkably well, as observed in Figs. 5 and 8 b. However, as predicted at the end of Section 4.2, the macroscale field performs poorly at short times, which is particularly obvious at the shortest time (t = 0.1tc ) (Figs. 5 and 8 b). Indeed the homogenization limit BCs are more accurate at short times for Case B as evident in Table 2. This inaccuracy at short times is even more prominent for a small number of layers as shown in Fig. 9. Clearly, this poor approximation is due to the large critical time [28] (time required for the microscale field in each layer to be approximately linear) after which Assumption 2 in Appendix A holds, a condition that must be true for the microscale perturbation to be proportional to the solution of the cell problem (8). The macroscale field in Case B (Fig. 5) performs well in the interior of the domain but not so well near the boundaries. In order to pass through the center of the microscale field oscillations, the macroscale field must be greater than one near

E.J. Carr et al. / Applied Mathematical Modelling 47 (2017) 600–618

613

Fig. 8. Reconstructed microscale fields u(x, t ) at t = 0.1tc , tc , 0.2, 1.0 for Cases A–D with corrected BCs, where tc = d12 /κ1 is the characteristic time of diffusion in the first layer.

Fig. 9. Reconstructed microscale fields u(x, t ) at t = 0.1tc , tc , 0.2, 1.0 for Cases A and B with corrected BCs and m = 4 layers, where tc = d12 /κ1 is the characteristic time of diffusion in the first layer.

the left boundary. The reason for this occurring is clear from the form of the corrected macroscale BCs (31) and (37). For example for Case B, with  bL < 0 and ∂∂Ux (0, t ) < 0 we have:

∂U U ( 0, t ) = 1 +  bL ( 0, t ) > 1. ∂x

614

E.J. Carr et al. / Applied Mathematical Modelling 47 (2017) 600–618

Fig. 10. Plot of the error for Case A versus the scale parameter ε = (lm − l0 )/m with m = 22 , 23 , . . . , 210 layers at t = 0.2 using for the benchmark solution: (a) the true microscale solution u(x, t), and (b) the averaged microscale solution u(x, t). The error metrics are defined in Eq. (29) with Nx = 5 used to produce the results shown.

Although it is not evident in the reconstructed microscale fields (Fig. 8b), this behaviour of the macroscale field is clearly non-physical for many real-world problems4 . The reason for this inconsistency comes back to the averaged microscale field depicted in Fig. 3b (this is the same problem as Case B), which does not obey the evolution law (12a) as argued in Section 4.2. It is for this reason that negative eigenvalues seem to infer that the macroscale Eq. (12a) should be modified. Fig. 10 plots the error for Case A versus the scale parameter ε = (lm − l0 )/m (equivalent to the layer width d) for a specified time t, showing that the error decreases for decreasing ε as expected. Using the slopes of the error curves, we can deduce convergence rates for the difference approaches. In particular, using the microscale field as the reference solution: • • •

the error of the reconstructed microscale field is O(ε 2 ), the error of the macroscale field with homogenization limit BCs is O(ε ), and the error of the macroscale field with corrected BCs is O(ε ),

while using the averaged microscale field as the reference solution: • •

the error of the macroscale field with homogenization limit BCs is O(ε ), and the error of the macroscale field with corrected BCs is O(ε 2 ).

The benefit of the corrected BCs is now obvious, providing superior convergence to the averaged microscale field as the scale parameter ε tends to zero. Specifically, the improvement in accuracy grows from an order of magnitude for ε = 2−4 (16 layers) to three orders of magnitude for ε = 2−10 (1024 layers). Fig. 11 depicts slight increasing errors for Case A as the problem becomes more heterogeneous, that is, as the heterogeneity ratio κ 2 /κ 1 increases above one. As is visible in Fig. 11a the reconstructed microscale field attains an error that is an order of magnitude less than both macroscale fields, across each value of the heterogeneity ratio tested. Moreover, using the microscale field as the reference solution, the macroscale field offers only a small improvement in accuracy under corrected BCs compared with homogenization limit BCs. Recall, however, that the aim of the macroscale field is to approximate the averaged microscale field. In this case, as can be seen in Fig. 11b, the corrected BCs provide an order of magnitude improvement in accuracy. 6. Summary and conclusions This paper has studied the form of the boundary conditions (BCs) utilized in the classical macroscale model of onedimensional multilayer diffusion. Using the method of volume averaging and an “average plus perturbation” decomposition of the microscale field, we derived a set of corrected BCs for the macroscale field. For microscale BCs of Dirichlet and Robin type and a finite number of layers, the corrected BCs differ from traditionally-imposed macroscale BCs, taking the form of appropriately-defined Robin conditions; however, both BCs are equivalent in the limit as the number of layers tends to infinity. Numerical results demonstrated that the corrected BCs produce a macroscale field that substantially improves approximation of the averaged microscale field, and leads to a reconstructed microscale field that is in excellent agreement with the true microscale field. Finally, the analysis and simulations presented in this paper suggest that for some problems 4 For example, in Case B suppose that u(x, t) denotes the temperature of a rod at position x and time t. With the given initial and boundary conditions, the temperature must be bounded as follows 0 ≤ u(x, t) ≤ 1 for all 0 ≤ x ≤ 1 and t ≥ 0.

E.J. Carr et al. / Applied Mathematical Modelling 47 (2017) 600–618

615

Fig. 11. Plot of the error for Case A versus the heterogeneity ratio κ 2 /κ 1 for m = 20 layers. The error is calculated at t = 5d12 /κ1 using as the benchmark solution: (a) the true microscale solution u(x, t), and (b) the averaged microscale solution u(x, t). The diffusivity κ 2 is fixed while κ 1 varies: κ2 = 1 and κ1 = 2−8 , 2−7 , . . . , 2−2 . The error metrics are defined in Eq. (29) with Nx = 5 used to produce the results shown.

Fig. A.12. Averaging interval [x − δ, x + δ ] consisting of r layers.

the classical macroscale equation may not be suitable and should be modified. This will be interesting to pursue in the future. Acknowledgements This research was funded by Australian Research Council (ARC) grant DE150101137. The third author acknowledges the School of Mathematical Sciences at QUT for travel support during his visit to Brisbane in 2015. Appendix A. Volume averaging In this appendix, we derive the macroscale (averaged) Eq. (7) by applying the method of volume averaging to the microscale model (5). For the most part, we follow the presentation of [9]. Consider an averaging interval (Fig. A.12), of length 2δ , centered at the point x and comprising r layers such that the interfaces between adjacent layers satisfy:

x − δ < lq < lq+1 < · · · < lq+r−2 < x + δ. For a piecewise function:

φ (x ) =

⎧ φ1 ( x ) ⎪ ⎪ ⎨φ2 (x )

l0  x  l1 l1  x  l2

. ⎪ ⎪ ⎩ . φm ( x )

lm−1  x  lm ,

.

616

E.J. Carr et al. / Applied Mathematical Modelling 47 (2017) 600–618

define the following averaging operator over the averaging interval:



1 φ(x, t ) = 2δ



lq

x −δ

φq ( ζ , t ) d ζ +

q+r−2  l  i i=q+1

li−1

φi ( ζ , t ) d ζ +



x+δ

lq+r−2

φq+r−1 (ζ , t ) dζ ,

or, equivalently, with ξ = ζ − x:



1 φ(x, t ) = 2δ



lq −x

−δ

φq ( x + ξ , t ) d ξ +

q+r−2  l −x  i i=q+1

li−1 −x

φi ( x + ξ , t ) d ξ +

 δ lq+r−2 −x

(A.1)

φq+r−1 (x + ξ , t ) dξ .

Averaging the microscale Eq. (5a) over the interval [x − δ, x + δ ] (Fig. A.12), we have:



    ∂u ∂u ∂ = . κ (x ) ∂t ∂x ∂x

where

u(x, t ) =

⎧ u (x, t ) ⎪ ⎪ 1 ⎨ u2 (x, t ) ⎪ ⎪ ⎩

l0  x  l1 l1  x  l2

.. .

um (x, t )

κ (x ) =

lm−1  x  lm ,

⎧ κ ⎪ ⎪ 1 ⎨ κ2

l0 < x < l1 l1 < x < l2

. ⎪ ⎪ ⎩. κm

lm−1 < x < lm .

.

Applying Leibniz’s rule for differentiation under the integral sign yields:

   q+r−2  ∂ u ∂ ∂u 1  ∂u ∂u κ (x ) κi i (li , t ) − κi+1 i+1 (li , t ) . = + ∂t ∂x ∂x 2δ ∂x ∂x i=q

The second term on the right-hand side of the above equation is equal to zero due to the interface conditions (5f), and hence the averaged microscale equation reduces to:

  ∂ u ∂ ∂u = . κ (x ) ∂t ∂x ∂x

(A.2)

We follow Davit et al. [9] and introduce the following decomposition of u(x, t) into an average component u(x, t) and a perturbation component  u(x, t ):

u(x, t ) = u(x, t ) +  u(x, t ),

(A.3)

or equivalently

ui (x, t ) = u(x, t ) +  ui (x, t ),

(A.4)

where li−1 < x < li (i = 1, . . . , m). Introducing the “average plus perturbation” decomposition (A.4) into the averaged Eq. (A.2) and the microscale Eq. (5a), we obtain the following two equations:

    ∂ u ∂ ∂ u ∂ ∂ u = + , l0 < x < lm , κ (x ) κ (x ) ∂t ∂x ∂x ∂x ∂x     ∂ u ∂  ui ∂ ∂ u ∂ ∂ ui + = + , li−1 < x < li , κ κ ∂t ∂t ∂x i ∂x ∂x i ∂x

(A.5)

(A.6)

for i = 1, . . . , m. Assumption 1.

∂ u is constant over the averaging interval: ∂x

∂ u ∂ u  (x + ξ , t )  (x, t ) ∂x ∂x for ξ ∈ (−δ, δ ). Under Assumption 1, the averaged Eq. (A.5) becomes

    ∂ u ∂ ∂ u ∂ ∂ u = + . κ (x ) κ(x ) ∂t ∂x ∂x ∂x ∂x

(A.7)

Subtracting (A.7) from (A.6) produces

      ∂ ui ∂ ∂ u ∂ ∂ u ∂ ∂ u κi i − κ (x ) = + . [κi − κ(x )] ∂t ∂x ∂x ∂x ∂x ∂x ∂x

(A.8)

E.J. Carr et al. / Applied Mathematical Modelling 47 (2017) 600–618

617

Assumption 2. The perturbation  ui (x, t ) is quasi-stationary:

∂ ui  all other terms in equation A.8. ∂t Assumption 3.

∂ u ∂ 2 u is constant in each layer or equivalently  0 in each layer. ∂x ∂ x2

Under Assumptions 2 and 3, Eq. (A.8) becomes:

    ∂ ∂ ui ∂ ∂ u d ∂ u  κ κ (x ) − + = 0. [κi − κ(x )] ∂x i ∂x ∂x ∂x dx ∂x

(A.9)

Let the perturbation component take the form:

 u(x, t ) = ψ (x )

∂ u  (x, t ), ∂x

(A.10)

where ψ (x) is piecewise-continuous:

ψ (x ) =

⎧ ψ1 ( x ) l0  x  l1 ⎪ ⎪ ⎨ψ2 (x ) l1  x  l2 .

. ⎪ ⎪ ⎩ . ψm (x ) lm−1  x  lm ,

so that the perturbation component in the ith layer is defined as

 ui (x, t ) = ψi (x )

∂ u (x, t ). ∂x

Substituting (A.10) into (A.9) and using Assumption 3 yields:

 ∂ u  d  − κi ψi (x ) dx ∂x

  ∂ ∂ u d ∂ u   κ ( x )ψ ( x ) + = 0. [κi − κ(x )] ∂x ∂x dx ∂x

(A.11)

Applying Assumptions 1 produces:

 ∂ u  d   ∂ u d d  ∂ u − + = 0, κi ψi (x ) κ ( x )ψ  ( x ) [κi − κ(x )] dx ∂x dx ∂x dx ∂x

(A.12)

and hence:

 d  d d  κi ψi (x ) − κ (x )ψ  (x ) + [κi − κ(x )] = 0. dx dx dx Assumption 4.

(A.13)

d κ(x )  0 . dx

Under Assumption 4, Eq. (A.13) reduces to:

 d   d    κi ψi (x ) + 1 − κ ( x )ψ  ( x ) = 0. dx dx

(A.14)

Assumption 5. The solutions ψ1 , . . . , ψm of (A.14) can be solved over the unit-cell [l0 , lp ] (see Section 2) with periodic BCs. Assumption 6.  u(x, t )  0. Under Assumption 5, Eq. (A.14) reduces to:

 d    κ ψ ( x ) + 1 = 0, dx k k

(A.15)

for k = 1, . . . , p. We have derived the periodic cell problem (9a)–(9e), the solution of which is unique up to an additive constant. Under Assumption 6 and noting the form of the perturbation component (A.10) we obtain the additional constraint (9f). Consider again the macroscale (averaged) Eq. (A.5). Substituting (A.10) and using Assumption 3 yields:

    ∂ u ∂ ∂ u ∂ ∂ u = + , κ (x ) κ ( x )ψ  ( x ) ∂t ∂x ∂x ∂x ∂x

and hence

    ∂ u ∂ u ∂  κ (x ) ψ  (x ) + 1  , ∂t ∂x ∂x

(A.16)

618

E.J. Carr et al. / Applied Mathematical Modelling 47 (2017) 600–618

under Assumption 1. Since both κ (x) and ψ (x) are periodic, the operator · reduces to the average over the unit-cell:

∂ u ∂  ∂t ∂x



p    1  lk κk ψk (x ) + 1 dx l p − l0 lk−1 k=1



 ∂ u . ∂x

(79)

κeff

Hence, we have derived the macroscale (averaged) Eq. (7) and the form of the effective diffusivity κ eff . References [1] J.-L. Auriault, Macroscopic modelling of heat transfer in composites with interfacial thermal barrier, Int. J. Heat Mass Tranf. 37 (18) (1994) 2885–2892. [2] E.J. Carr, I.W. Turner, Two-scale computational modelling of water flow in unsaturated soils containing irregular-shaped inclusions, Int. J. Numer. Methods Eng. 98 (3) (2014) 157–173. [3] J. Bear, A.-D. Cheng, Modeling Groundwater Flow and Contaminant Transport, Springer, 2010. [4] P. Jenny, S.H. Lee, H.A. Tchelepi, Adaptive multiscale finite-volume method for multiphase flow and transport in porous media, Multiscale Model. Simul. 3 (1) (2004) 50–64. [5] E.J. Carr, P. Perré, I.W. Turner, The extended distributed microstructure model for gradient-driven transport: a two-scale model for bypassing effective parameters, J. Comput. Phys. 327 (2016) 810–829. [6] J.-L. Auriault, Heterogeneous medium: is an equivalent macroscopic description possible? Int. J. Eng. Sci. 29 (7) (1991) 785–795. [7] A. Bensoussan, J.L. Lions, G. Papanicolaou, Asymptotic analysis for periodic structures, North-Holland Publishing Company, 1978. [8] S. Whitaker, Simultaneous heat, mass and momentum transfer in porous media: a theory of drying, in: J.P. Hartnett, T.F. Irvine (Eds.), Advances in Heat Transfer, 13, Elsevier, 1977, pp. 119–203. [9] Y. Davit, C.G. Bell, H.M. Byrne, L.A.C. Chapman, L.S. Kimpton, G.E. Lang, K.H.L. Leonard, J.M. Oliver, N.C. Pearson, R.J. Shipley, S.L. Waters, J.P. Whiteley, B.D. Wood, M. Quintard, Homogenization via formal multiscale asymptotics and volume averaging: how do the two techniques compare? Adv. Water Resour. 62 (2013) 178–206. [10] G.A. Pavliotis, A.M. Stuart, Multiscale Methods: Averaging and Homogenization, Springer, 2008. [11] C.C. Mei, B. Vernescu, Homogenization Methods for Multiscale Mechanics, World Scientific Publishing Co., New Jersey, USA, 2010. [12] M. Prat, On the boundary conditions at the macroscopic level, Transp. Porous Med. 4 (1989) 259–280. [13] U. Hornung, Homogenization and Porous Media, Springer–Verlag, New York, 1997. [14] M.H. Holmes, Introduction to Perturbation Methods, Springer-Verlag, 1995. [15] J. Lewandowska, A. Szymkiewicz, K. Burzynksi, M. Vauclin, Modeling of unsaturated water flow in double-porosity soils by the homogenization approach, Adv. Water Resour. 27 (2004) 283–296. [16] C. Chen, Multiscale Modelling of Continuum and Discrete Dynamics in Materials with Complicated Microstructure, The University of Adelaide, 2015 Ph.D. thesis. [17] C. Chen, A.J. Roberts, J.E. Bunder, Macroscale boundary conditions for a non-linear heat exchanger, ANZIAM J. 56 (CTAC2014) (2015) C16–C31. [18] C. Chen, A.J. Roberts, J.E. Bunder, The macroscale boundary conditions for diffusion in a material with microscale varying diffusivities, ANZIAM J. 55 (EMAC2013) (2014) C218–C234. [19] G. Allaire, M. Amar, Boundary layer tails in periodic homogenization, ESAIM Control Optim. Calc. Var. 4 (1999) 209–243. [20] A. Matine, N. Boyard, P. Cartraud, G. Legrain, Y. Jarny, Modeling of thermophysical properties in heterogeneous periodic media according to a multi-scale approach: effective conductivity tensor and edge effects, Int. J. Heat Mass Transf. 62 (2013) 586–603. [21] D. Gérard-Varet, N. Masmoudi, Homogenization and boundary layers, Acta Mathematica 209 (2012) 133–178. [22] M. Neuss-Radu, The boundary behaviour of a composite material, Math. Model. Numer. Anal. 35 (3) (2001) 407–435. [23] R.I. Hickson, S.I. Barry, G.N. Mercer, Critical times in multilayer diffusion. Part 1: exact solutions, Int. J. Heat Mass Transf. 52 (2009) 5776–5783. [24] E.J. Carr, I.W. Turner, A semi-analytical solution for mulitlayer diffusion in a composite medium consisting of a large number of layers, Appl. Math. Model. 40 (15–16) (2016) 7034–7050. [25] W.A. Strauss, Partial Differential Equations: An Introduction, John Wiley & Sons, Inc., 1992. [26] D.W. Trim, Applied Partial Differential Equations, PWS-Kent Pub. Co, 1990. [27] E.J. Carr, N.G. March, A semi-analytical solution to multilayer diffusion problems with transient boundary conditions and general interface conditions, Submitted for publication (2016). [28] R.I. Hickson, S.I. Barry, H.S. Sidhu, G.N. Mercer, A comparison of critical time definitions in multilayer diffusion, ANZIAM J. 52 (2011) 333–358.