Advection improves homogenized models of continuum diffusion in one-dimensional heterogeneous media

TL;DR

This paper introduces an advection term into homogenized diffusion models for one-dimensional heterogeneous media, resulting in more accurate approximations of the original complex system.

Contribution

It presents a novel method to compute effective advection and diffusion coefficients by solving boundary value problems, enhancing homogenized models of heterogeneous diffusion.

Findings

Including advection improves model accuracy

Effective coefficients depend on boundary conditions

Method outperforms standard homogenization approaches

Abstract

We propose an alternative method for one-dimensional continuum diffusion models with spatially variable (heterogeneous) diffusivity. Our method, which extends recent work on stochastic diffusion, assumes the constant-coefficient homogenized equation takes the form of an advection-diffusion equation with effective (diffusivity and velocity) coefficients. To calculate the effective coefficients, our approach involves solving two uncoupled boundary value problems over the heterogeneous medium and leads to coefficients depending on the spatially-varying diffusivity (as usual) as well as the boundary conditions imposed on the heterogeneous model. Computational experiments comparing our advection-diffusion homogenized model to the standard homogenized model demonstrate that including an advection term in the homogenized equation leads to improved approximations of the solution of the original heterogeneous model.

Tables

Table 1. Table 1: Problem descriptions for the nine test cases used in the computational experiments of Section IV . For each test case, this table identifies the spatially-variable diffusivity D ​ ( x ) 𝐷 𝑥 D(x) , initial solution f ​ ( x ) 𝑓 𝑥 f(x) , and boundary values g 0 ​ ( t ) subscript 𝑔 0 𝑡 g_{0}(t) and g L ​ ( t ) subscript 𝑔 𝐿 𝑡 g_{L}(t) appearing in the heterogeneous model ( 6 )–( 8 ).

Case	Description
A	$D\left(x\right)=0.5+0.2\sin \left(20x\right)$, $f\left(x\right)=0$, $g_0\left(t\right)=1$ and $g_L\left(t\right)=0$.
B	$D\left(x\right)=0.8-0.6x+0.2\sin \left(20x\right)$, $f\left(x\right)=0$, $g_0\left(t\right)=1$ and $g_L\left(t\right)=0$.
C	$D\left(x\right)=0.5+0.2\sin \left(x/\epsilon \right)$ with $\epsilon =0.005$, $f\left(x\right)=0$, $g_0\left(t\right)=1$ and $g_L\left(t\right)=0$.
D	$D\left(x\right)=0.1\exp \left(2.2x^2\right)+0.05\sin \left(80x\right)$, $f\left(x\right)=0$, $g_0\left(t\right)=1$ and $g_L\left(t\right)=0$.
E	$D\left(x\right)=0.5+0.5\exp \left(x-1\right)\cos \left(80\left(1-x\right)\right)$, $f\left(x\right)=0$, $g_0\left(t\right)=1$ and $g_L\left(t\right)=0$.
F	$D\left(x\right)$ is piecewise linear on the sub-intervals $(\left(i-1\right)H,iH)$ where $i=1,\dots ,50$ and $H=1/50$. At locations $x_i=iH$ for $i=0,\dots ,50$, $D\left(x_i\right)=0.9-0.8/\left(1+\exp \left(20\left(x-0.5\right)\right)\right)+\epsilon$, where $\epsilon$ is a random number generated from a normal distribution with support $[-0.05,0.05]$. For all $i=1,\dots ,50$, within sub-interval $(\left(i-1\right)H,iH)$ the linear form of $D\left(x\right)$ is constructed to interpolate the values previously assigned at $x_{i-1}=\left(i-1\right)H$ and $x_i=iH$. $f\left(x\right)=0$, $g_0\left(t\right)=1$ and $g_L\left(t\right)=0$.
G	$D\left(x\right)=0.5+0.24\left(\sin \left(20x\right)+\sin \left(80x\right)\right)$, $f\left(x\right)=0$, $g_0\left(t\right)=1$ and $g_L\left(t\right)=0.75\left(1-\exp \left(-25t\right)\right)$.
H	$D\left(x\right)$ is piecewise constant on the sub-intervals $(\left(i-1\right)H,iH)$ where $i=1,\dots ,16$ and $H=1/16$. Within each sub-interval the constant value of $D\left(x\right)$ is assigned randomly from a uniform distribution with support $[0.01,0.99]$. $f\left(x\right)=\exp \left(-30{\left(x-0.5\right)}^2\right)$, $g_0\left(t\right)=0$ and $g_L\left(t\right)={10}^{-6}$.
I	$D\left(x\right)$ is piecewise linear on the sub-intervals $(\left(i-1\right)H,iH)$ where $i=1,\dots ,24$ and $H=1/24$. At locations $x_i=iH$ for $i=0,\dots ,24$, we assign a random value of the diffusivity generated from a uniform distribution with support $[0.01,0.99]$. For all $i=1,\dots ,24$, within sub-interval $(\left(i-1\right)H,iH)$ the linear form of $D\left(x\right)$ is constructed to interpolate the random values previously assigned at $x_{i-1}=\left(i-1\right)H$ and $x_i=iH$. $f\left(x\right)=2x$ if $x\in [0,0.5]$ otherwise $f\left(x\right)=2-2x$ if $x\in [0.5,1]$, $g_0\left(t\right)=0$ and $g_L\left(t\right)=0.5$.

Table 2. Table 2: Mean absolute errors ( 33 )–( 34 ) and effective coefficients for the new homogenized model ( 9 )–( 11 ) and standard homogenized model ( 29 )–( 32 ).

	Homogenized advection-diffusion model			Homogenized diffusion-only model
Case	$D_{\mathrm{eff}}$	$v_{\mathrm{eff}}$	Error	${\tilde{D}}_{\mathrm{eff}}$	Error
A	$0.497$	$0.079$	$1.48\text{×}{10}^{-02}$	$0.464$	$1.81\text{×}{10}^{-02}$
B	$0.409$	$0.805$	$2.36\text{×}{10}^{-02}$	$0.345$	$1.38\text{×}{10}^{-01}$
C	$0.462$	$0.009$	$1.51\text{×}{10}^{-03}$	$0.459$	$1.91\text{×}{10}^{-03}$
D	$0.160$	$-0.319$	$1.36\text{×}{10}^{-02}$	$0.162$	$7.22\text{×}{10}^{-02}$
E	$0.318$	$0.374$	$1.91\text{×}{10}^{-02}$	$0.324$	$6.57\text{×}{10}^{-02}$
F	$0.186$	$-0.622$	$2.06\text{×}{10}^{-02}$	$0.193$	$1.20\text{×}{10}^{-01}$
G	$0.342$	$0.097$	$2.93\text{×}{10}^{-02}$	$0.306$	$3.18\text{×}{10}^{-02}$
H	$0.312$	$-0.317$	$1.69\text{×}{10}^{-02}$	$0.215$	$4.32\text{×}{10}^{-02}$
I	$0.427$	$0.062$	$1.51\text{×}{10}^{-02}$	$0.381$	$1.87\text{×}{10}^{-02}$

Equations

$$\frac{\partial u}{\partial t}=\frac{\partial }{\partial x}\left(D(x)\frac{\partial u}{\partial x}\right),$$

$$\frac{\partial U}{\partial t}=\frac{\partial }{\partial x}\left(D_{\mathrm{eff}}\frac{\partial U}{\partial x}\right),$$

$$D_{\mathrm{eff}}=\frac{L}{{\int }_0^{L}D(x{)}^{-1}\mathrm{d}x}.$$

$$\frac{\partial u}{\partial t}=D(x)\frac{{\partial }^{2}u}{\partial x^2}+D^{\prime }(x)\frac{\partial u}{\partial x}.$$

$$\frac{\partial U}{\partial t}=D_{\mathrm{eff}}\frac{{\partial }^{2}U}{\partial x^2}-v_{\mathrm{eff}}\frac{\partial U}{\partial x},$$

$$\frac{\partial u}{\partial t}=\frac{\partial }{\partial x}\left(D(x)\frac{\partial u}{\partial x}\right),$$

$$u(x,0)=f(x),$$

$$u(0,t)=g_0(t),u(L,t)=g_L(t),$$

$$\frac{\partial U}{\partial t}=D_{\mathrm{eff}}\frac{{\partial }^{2}U}{\partial x^2}-v_{\mathrm{eff}}\frac{\partial U}{\partial x},$$

$$U(x,0)=f(x),$$

$$U(0,t)=g_0(t),U(L,t)=g_L(t).$$

$${\int }_0^{L}s(x)\text{d}x={\int }_0^{L}S(x)\text{d}x,$$

$$\frac{\mathrm{d}}{\mathrm{d}x}\left(D(x)\frac{\mathrm{d}s}{\mathrm{d}x}\right)=0,$$

$$s(0)=\underset{t\to \infty }{\lim }{g}_0(t),s(L)=\underset{t\to \infty }{\lim }{g}_L(t),$$

$$D_{\mathrm{eff}}\frac{{\mathrm{d}}^{2}S}{\mathrm{d}{x}^2}-v_{\mathrm{eff}}\frac{\mathrm{d}S}{\mathrm{d}x}=0,$$

$$S(0)=\underset{t\to \infty }{\lim }{g}_0(t),S(L)=\underset{t\to \infty }{\lim }{g}_L(t).$$

$${\int }_0^{L}w(x)\text{d}x={\int }_0^{L}W(x)\text{d}x,$$

$$w(x)$$

$$W(x)$$

$$\frac{\mathrm{d}}{\mathrm{d}x}\left(D(x)\frac{\mathrm{d}w}{\mathrm{d}x}\right)=r(x),$$

$$w(0)={\int }_0^{\infty }\left[s(0)-g_0(t)\right]\mathrm{d}t,$$

$$w(L)={\int }_0^{\infty }\left[s(L)-g_L(t)\right]\mathrm{d}t,$$

$$D_{\mathrm{eff}}\frac{{\mathrm{d}}^{2}W}{\mathrm{d}{x}^2}-v_{\mathrm{eff}}\frac{\mathrm{d}W}{\mathrm{d}x}=R(x),$$

$$W(0)={\int }_0^{\infty }\left[S(0)-g_0(t)\right]\mathrm{d}t,$$

$$W(L)={\int }_0^{\infty }\left[S(L)-g_L(t)\right]\mathrm{d}t,$$

$$\mathbf{F}(\mathbf{c})={\left(F_1(\mathbf{c}),F_2(\mathbf{c})\right)}^T=0,$$

$$F_1(\mathbf{c})$$

$$F_2(\mathbf{c})$$

$$\frac{\partial \tilde{U}}{\partial t}={\tilde{D}}_{\mathrm{eff}}\frac{{\partial }^2\tilde{U}}{\partial x^2},$$

$$\tilde{U}(x,0)=f(x),$$

$$\tilde{U}(0,t)=g_0(t),\tilde{U}(L,t)=g_L(t),$$

$$\text{with Deff=∫0LD(x)−1dxL}.$$

$$\displaystyle\text{Error}=\frac{1}{N_{t}N_{x}}\sum_{j=1}^{N_{t}}\sum_{k=1}^{N_{x}}\bigl{|}U_{k}^{j}-u_{k}^{j}\bigr{|},$$

$$\displaystyle\text{Error}=\frac{1}{N_{t}N_{x}}\sum_{j=1}^{N_{t}}\sum_{k=1}^{N_{x}}\bigl{|}\widetilde{U}_{k}^{j}-u_{k}^{j}\bigr{|},$$

Full text

I Introduction

Many heat and mass transport modelling problems involve heterogeneous media exhibiting spatial variability in material properties. Some examples include water and pollutant transport in groundwater aquifers composed of soils and rocks Chen and Ren (2008), heat and moisture transport within wood during drying Perré and Turner (2002) and oxygen transport in biological tissues Matzavinos and Ptashnyk (2016). When such problems involve material properties that vary rapidly relative to the size of the problem domain, direct computation of the solution is expensive since one has to discretise the domain using a very fine mesh in order to capture the heterogeneity. A popular strategy for alleviating such computational issues is to formulate a simpler homogenized model with constant coefficients that smoothes out the heterogeneity while remaining accurate Chen and Ren (2008); Huysmans and Dassargues (2007); Auriault (1991); Roberts (2010); Davit et al. (2013); Carr et al. (2017); Abdulle and E (2003); Samaey et al. (2005).

In this paper, we consider the one-dimensional diffusion equation in a heterogeneous medium $(0,L)$:

[TABLE]

where $D(x)>0$ is the spatially varying diffusivity. Our goal is to approximate the smooth or average behaviour of $u(x,t)$ by the solution of a simpler equation with spatially constant coefficients. The natural approach is to use a simpler equation of the form

[TABLE]

where $D_{\mathrm{eff}}>0$ is a constant effective, equivalent or homogenized diffusivity chosen so that $U(x,t)$ provides an accurate approximation to $u(x,t)$. The standard choice for $D_{\mathrm{eff}}$ is the harmonic average of $D(x)$ (see, e.g., Carr et al. (2017); Crank (1975); Hornung (1997); Holmes (2013); Pavliotis and Stuart (2008); Roberts (2010); Ray et al. (2018)):

[TABLE]

This definition possesses a strong theoretical foundation being the result of applying the methods of homogenization by asymptotic expansion and volume averaging Carr et al. (2017); Davit et al. (2013); Holmes (2013); Pavliotis and Stuart (2008) to the heterogeneous equation (1). The definition (3) is also physically intuitive, as it can be derived by considering the boundary value problem consisting of the steady state analogue of the diffusion equation (1) paired with the boundary conditions $u(0)=0$ and $u(L)=L\frac{\partial U}{\partial x}$ Hornung (1997), which impose a (constant) macroscopic gradient of $\frac{\partial U}{\partial x}$ over the medium. Solving this boundary value problem yields a solution exhibiting a constant (homogenized) flux over the medium taking the form of $q=-D_{\mathrm{eff}}\frac{\partial U}{\partial x}$, where $D_{\mathrm{eff}}$ is defined as in Eq (3) Hornung (1997).

In Figure 1, we plot the solution of the homogenized equation (2)–(3) and the heterogeneous equation (1) for two choices of $D(x)$ and a particular set of initial and boundary conditions. For the first problem, the homogenized solution $U(x,t)$ provides a good approximation to the heterogeneous solution $u(x,t)$ (Figure 1(a)). However, for the second problem, which exhibits advective behaviour in the positive $x$ direction (Figure 1(b)), $U(x,t)$ poorly approximates $u(x,t)$ since the homogenized equation (2) cannot capture such behaviour. Note this behaviour becomes apparent when differentiating the diffusive flux in Eq (1):

[TABLE]

To address this issue, in this paper, we present a new homogenization approach for Eq (1). Motivated by the results in Figure 1, the observation (4) and recent work on homogenization of random-walks through heterogeneous media Carr and Simpson (2019), our approach assumes the homogenized equation takes the form of an advection-diffusion equation:

[TABLE]

where $v_{\mathrm{eff}}$ is an effective velocity. To calculate the effective coefficients ($D_{\mathrm{eff}}$ and $v_{\mathrm{eff}}$), our proposed method involves applying two constraints enforcing equality of appropriate measures of the spatial and temporal behaviour of the respective homogenized and heterogeneous models. This procedure requires the solution of two uncoupled boundary value problems over the heterogeneous medium and in contrast to the harmonic average definition (3) leads to effective coefficients depending on the boundary conditions imposed within the heterogeneous model.

The rest of this paper is structured as follows. In section II, we outline the specific heterogeneous model considered in this work including initial and boundary conditions. We then present the advection-diffusion homogenized model and describe how the effective coefficients are calculated (section III). In section IV, computational experiments assessing the accuracy of the homogenized model are presented. Conclusions and a summary of the work are given in section V along with possible avenues for future research.

II Heterogeneous model

We consider the heterogeneous diffusion model:

[TABLE]

where $u(x,t)$ is the unknown scalar field (heterogeneous or benchmark solution), $D(x)$ is the specified spatially-varying diffusivity, $L$ is the length of the medium, $f(x)$ provides the initial solution and $g_{0}(t)$ and $g_{L}(t)$ are known functions specifying the boundary values of the solution. We make the assumption that the limits $\lim_{t\rightarrow\infty}g_{0}(t)$ and $\lim_{t\rightarrow\infty}g_{L}(t)$ exist.

III Homogenized model

Our homogenization approach seeks to approximate the solution $u(x,t)$ of the heterogeneous model (6)–(8) by the solution $U(x,t)$ of the following advection-diffusion model with homogenized (constant) coefficients $D_{\mathrm{eff}}$ and $v_{\mathrm{eff}}$:

[TABLE]

To determine the two unknown effective coefficients, $D_{\mathrm{eff}}$ and $v_{\mathrm{eff}}$, we specify two constraints. First, we enforce equality of the spatial-average of the steady-state solutions of the heterogeneous (6)–(8) and homogenized (9)–(11) models. We express this constraint as follows:

[TABLE]

where $s(x)$ is the steady state solution of the heterogeneous model (6)–(8) satisfying the boundary value problem:

[TABLE]

and $S(x)$ is the steady state solution of the homogenized model (9)–(11) satisfying the boundary value problem:

[TABLE]

The constraint (12) ensures $U(x,t)$ accurately exhibits the averaged spatial behaviour of $u(x,t)$ at steady-state. To match the temporal behaviour of $U(x,t)$ and $u(x,t)$ we enforce:

[TABLE]

where $w(x)$ and $W(x)$ are defined as:

[TABLE]

These quantities measure the (signed) distance between the steady-state and transient solutions of the heterogeneous and homogenized models over time, respectively, with the averaged values of $w(x)$ and $W(x)$ providing a simple way to quantify the temporal behaviour of each model between initial and steady state. Attractively, both $w(x)$ and $W(x)$ can be calculated without explicit calculation of $u(x,t)$ and $U(x,t)$ appearing in the definitions Carr (2019). Following previous work Carr (2019); Carr and Simpson (2018); Ellery et al. (2012), applying the linear operator $\mathcal{L}\varphi:=\frac{\partial}{\partial x}\left(D(x)\frac{\partial\varphi}{\partial x}\right)$ to both sides of the definition (18) and making use of the differential equations (6) and (13) yields the following boundary value problem satisfied by $w(x)$:

[TABLE]

where $r(x)=f(x)-s(x)$. We remark that the boundary conditions (21)–(22) follow directly from the definition (18) and the heterogeneous model boundary conditions (8) with the imposed values of $w(0)$ and $w(L)$ calculated by directly evaluating the integrals since the values of $s(0)$ and $s(L)$ (14) and boundary functions $g_{0}(t)$ and $g_{L}(t)$ are known. In a similar manner to that carried out for $w(x)$, except now with the linear operator defined by $\mathcal{L}\varphi:=D_{\mathrm{eff}}\frac{\partial^{2}\varphi}{\partial x}-v_{\mathrm{eff}}\frac{\partial\varphi}{\partial x}$, the following boundary value problem is derived for $W(x)$:

[TABLE]

where $R(x)=f(x)-S(x)$.

Since the solutions $S(x)$ and $W(x)$ will depend nonlinearly on $D_{\mathrm{eff}}$ and $v_{\mathrm{eff}}$, the two constraints (12) and (17) together define a pair of coupled nonlinear equations:

[TABLE]

whose solution $\mathbf{c}=\left[D_{\mathrm{eff}},v_{\mathrm{eff}}\right]^{T}$ provides the effective coefficients supplied to the homogenized model (9)–(11). The form of the component functions, $F_{1}$ and $F_{2}$, are formulated by solving the various boundary value problems numerically. In this work, we employ a vertex-centered finite volume method on a uniform grid consisting of $N_{x}$ nodes with node spacing $h=L/(N_{x}-1)$. Let $s_{k}$, $S_{k}$, $w_{k}$ and $W_{k}$ denote the numerical approximation to the solutions of the boundary value problems (13)–(14), (15)–(16), (20)–(22) and (23)–(25) at $x=x_{k}:=(k-1)h$ for $k=1,\ldots,N_{x}$. Using these solutions, a simple trapezoidal rule is applied to evaluate the integrals in the constraint equations (12) and (17) yielding the component functions:

[TABLE]

Note the values of $S_{k}$ and $W_{k}$ depend nonlinearly on $D_{\mathrm{eff}}$ and $v_{\mathrm{eff}}$ so $F_{1}$ and $F_{2}$ are nonlinear functions of $\mathbf{c}$.

IV Computational experiments

We now compare the solution of the advection-diffusion homogenized model (9)–(11), $U(x,t)$, to the benchmark solution, $u(x,t)$, of the heterogeneous model (6)–(8). We also compare these solutions to the solution of the standard diffusion-only homogenized model, $\widetilde{U}(x,t)$, where the effective equation takes the form of the diffusion equation with constant harmonic-averaged effective diffusivity:

[TABLE]

To obtain the the effective coefficients $D_{\mathrm{eff}}$ and $v_{\mathrm{eff}}$ for the advection-diffusion homogenized model (9)–(11), we solve the nonlinear system (26)–(28) using MATLAB’s in-built fsolve function. To solve the heterogeneous model (6)–(8), the advection-diffusion homogenized model (9)–(11) and the diffusion-only homogenized model (29)–(32), we use a vertex-centered finite volume method on a uniform grid consisting of $N_{x}$ nodes with node spacing $h=L/(N_{x}-1)$ (as in Section III). The resulting system of differential equations is then solved numerically using MATLAB’s in-built ode15s function. Further implementation details can be found in our code available on GitHub: https://github.com/elliotcarr/Carr2019b. Let $u_{k}^{j}$, $U_{k}^{j}$ and $\widetilde{U}_{k}^{j}$ denote the resulting numerical approximations to $u(x,t)$, $U(x,t)$ and $\widetilde{U}(x,t)$ at $x=x_{k}:=(k-1)h$ for $k=1,\ldots,N_{x}$ and $t=t_{j}:=j\tau$ for $j=1,\ldots,N_{t}$, where $\tau>0$ and $N_{t}\in\mathbb{N}^{+}$. Using these solutions we define the following mean absolute errors:

[TABLE]

for the the advection-diffusion homogenized model (9)–(11) and

[TABLE]

for the diffusion-only homogenized model (29)–(32).

In our computational experiments, we consider the nine test cases described in Table 1. In Figure 2 and Table 2, we report results for $N_{x}=1001$, $N_{t}=100$ and $\tau=0.01$. An immediate observation from these results is that the advection-diffusion homogenized model (9)–(11) is superior at capturing the smooth behaviour of the heterogeneous model (6)–(8) compared to the diffusion-only homogenized model (29)–(32). This is demonstrated in Figures 2(j)–(r) by the observation that $U(x,t)$ agrees with $u(x,t)$ better than $\widetilde{U}(x,t)$ does, and in Table 2 by the smaller values of the mean absolute error for the advection-diffusion homogenized model across all nine test cases.

Several other interesting observations are evident from the individual test cases. Cases A and B involve a constant and linearly decreasing diffusivity perturbed by a sinusoidal function, respectively, and are identical to the two problems shown in Figure 1. For case A, both homogenized models provide a good approximation to the solution behaviour of the heterogeneous model (see Figure 2(j)). However, for case B, the advection-diffusion homogenized model significantly outperforms the standard diffusion-only homogenized model. In this test case, the benchmark heterogeneous model produces a solution exhibiting advective behaviour in the positive $x$ direction, which is especially evident from the steady-state solution (see Figure 2(k)). This behaviour is accurately captured by the advection-diffusion equation (9) (note the large value of $v_{\mathrm{eff}}$ for case B compared to case A) but cannot be captured by the standard homogenized equation (29), where the steady-state solution is linear regardless of the value of $\widetilde{D}_{\mathrm{eff}}$ (Figure 2(k)). For case B, the new homogenized model also provides a superior match at earlier times (Figure 2(k)). Case C is a classical homogenization problem with $D(x)$ expressible as a function of $x/\varepsilon$, where $\varepsilon$ is a small parameter ($\varepsilon=0.005$ in this case). For this test case, the standard homogenized model (29)–(32) is obtained in the homogenization limit $\varepsilon\rightarrow 0$ with $U(x,t)$ and $\widetilde{U}(x,t)$ matching almost precisely with one another (Figure 2(l)). Case C demonstrates the self-averaging property of the diffusion equation where the solutions appear homogeneous when the characteristic length of the heterogeneities is very small relative to the length of the heterogeneous medium Puzko and Merzlikin (2019). The results for cases A, B and C demonstrate that differences between the advection-diffusion homogenized model and the standard diffusion-only homogenized model are accentuated when $D^{\prime}(x)$ in Eq (4) is large in magnitude relative to $D(x)$. Cases D, E and F further illustrate this observation with large values of $v_{\mathrm{eff}}$ in Table 2 and pronounced differences between the two homogenized models in Figures 2(m)–(o). Finally, cases G, H and I demonstrate that our homogenization approach performs well for time-dependent boundary conditions, non-uniform initial conditions, and non-monotone transitions from initial to steady-state, respectively (see Figure 2(p)–(r) and Table 2).

V Conclusions

Determining an effective homogeneous medium providing the smooth/averaged behaviour of a continuum transport process across a heterogeneous medium is a classical problem in many scientific and engineering disciplines. In this paper, we have presented a homogenization approach for the one-dimensional diffusion equation with spatially variable (heterogeneous) diffusivity. Our approach is noteworthy for the inclusion of an effective advection term in the homogenized equation in addition to the standard effective diffusion term. To identify the effective diffusivity and effective velocity we enforce equality of the spatial average of the steady-state solution of the homogenized and heterogeneous models and equality of the spatial average of a quantity characterising the temporal behaviour of both models. Our proposed homgenization approach requires the solution of two uncoupled boundary value problems Eqs (13)–(14) and (20)–(22), over the heterogeneous medium and the solution of a small system of nonlinear equations for the effective parameters (26)–(28). Due to the appearance of $g_{0}(t)$ and $g_{L}(t)$ in the boundary value problems, e.g. Eqs (23)–(25), the computed effective coefficients depend not only on the spatially-varying diffusivity but also on the boundary conditions imposed in the heterogeneous model.

Computational experiments carried out in Section IV demonstrate, for nine test cases, that the inclusion of an effective advection term leads to an improved approximation of the smooth behaviour of continuum diffusion in a one-dimensional heterogeneous medium. Possible avenues for future work include extension to two or three-dimensional problems or non-Dirichlet boundary conditions. For higher-dimensional problems, additional constraints are required to identify the additional diffusivities and velocities present in two and three dimensions. One way forward could be to note that the quantities $w(x)$ (18) and $W(x)$ (19) can be interpreted as the zeroth temporal moment of $s(x)-u(x,t)$ and $S(x)-U(x,t)$, respectively, and introduce constraints on the corresponding higher-order moments Carr and Simpson (2018). For non-Dirichlet boundary conditions, additional thought is required as the constraint enforcing equality of the spatial average of the steady-state solution (12) provides no information for certain choices of boundary conditions (8), e.g., for $u(0,t)=1$ and $\frac{\partial u}{\partial x}(L,t)=0$, the steady-state solutions, $s(x)$ and $S(x)$, are uniformly equal to one for all $D(x)$, $D_{\mathrm{eff}}$ and $v_{\mathrm{eff}}$.