A fast algorithm for semi-analytically solving the homogenization boundary value problem for block locally-isotropic heterogeneous media

TL;DR

This paper introduces a semi-analytical algorithm that efficiently computes effective diffusivity in heterogeneous media, reducing computational time while maintaining accuracy, thus enabling faster diffusion simulations in complex geometries.

Contribution

The paper presents a novel semi-analytical solution method for the boundary value problem in homogenization, outperforming standard finite volume methods in speed and accuracy.

Findings

Achieves equivalent accuracy to finite volume methods

Reduces computational time for homogenization calculations

Effectively handles complex heterogeneous geometries

Abstract

Direct numerical simulation of diffusion through heterogeneous media can be difficult due to the computational cost of resolving fine-scale heterogeneities. One method to overcome this difficulty is to homogenize the model by replacing the spatially-varying fine-scale diffusivity with an effective diffusivity calculated from the solution of an appropriate boundary value problem. In this paper, we present a new semi-analytical method for solving this boundary value problem and computing the effective diffusivity for pixellated, locally-isotropic, heterogeneous media. We compare our new solution method to a standard finite volume method and show that equivalent accuracy can be achieved in less computational time for several standard test cases. We also demonstrate how the new solution method can be applied to complex heterogeneous geometries represented by a two-dimensional grid of rectangular blocks. These results indicate that our new semi-analytical method has the potential to significantly speed up simulations of diffusion in heterogeneous media.

Figures

Tables

Table 1. Table 1: Effective diffusivity for cases 1–4 as calculated using the semi-analytical and finite volume methods and compared to the results presented in [ 7 ] . The principal directions of diffusion are overlaid upon each geometry. All results are reported to three significant figures, consistent with [ 7 ] .

	Case 1	Case 2
Geometry
Szymkiewicz [7]
Finite Volume
Semi-Analytical
	Case 3	Case 4
Geometry
Szymkiewicz [7]
Finite Volume
Semi-Analytical

Table 2. Table 2: Error and runtimes of the effective diffusivities calculated using the semi-analytical and finite volume methods.

Case A	Semi-Analytical		Finite Volume
$N_x$	Relative Error ($𝐄$)	Time (s)	Relative Error ($𝐄$)	Time (s)
4		0.00805		0.00938
8		0.0112		0.0280
16		0.0338		0.113
32		0.0630		0.531
64		0.269		2.93
Case B	Semi-Analytical		Finite Volume
$N_x$	Relative Error ($𝐄$)	Time (s)	Relative Error ($𝐄$)	Time (s)
4		0.00747		0.00923
8		0.0109		0.0277
16		0.0331		0.115
32		0.0629		0.530
64		0.270		2.92
Case C	Semi-Analytical		Finite Volume
$N_x$	Relative Error ($𝐄$)	Time (s)	Relative Error ($𝐄$)	Time (s)
4		0.0702		0.124
8		0.151		0.549
16		0.619		2.98
Case D	Semi-Analytical		Finite Volume
$N_x$	Relative Error ($𝐄$)	Time (s)	Relative Error ($𝐄$)	Time (s)
4		0.0696		0.0650
8		0.151		0.388
16		0.615		2.43

Table 3. Table 3: Effective diffusivity for four randomly generated aggregated media consisting of a 10 10 10 by 10 10 10 , 20 20 20 by 20 20 20 , 50 50 50 by 50 50 50 and 100 100 100 by 100 100 100 array of blocks, respectively, as calculated using the semi-analytical method. In all four cases N x = N y = 8 subscript 𝑁 𝑥 subscript 𝑁 𝑦 8 N_{x}=N_{y}=8 and the number of blocks with diffusivity 0.1 (dark grey blocks) is equal to the number of blocks with diffusivity 1 (light grey blocks). The aggregated geometries were inspired by those presented in [ 23 ] .

	$10\times 10$	$20\times 20$
Geometry
Effective diffusivity
	$50\times 50$	$100\times 100$
Geometry
Effective diffusivity

Table 4. Table 4: Pixelated geometries corresponding to Figure 1 (a) with corresponding effective diffusivities ( 𝐃 eff subscript 𝐃 eff \mathbf{D}_{\text{eff}} ), relative errors ( 𝐄 𝐄 \mathbf{E} ), number of rows of blocks in the pixelated geometry( r 𝑟 r ), number of rows in the corresponding linear system ( 53 ) ( N 𝑁 N ) and time taken to solve the linear system (Time(s)).

Benchmark
Geometry	${𝐃}_{\text{eff}}$	$𝐄$	$N$	Time(s)
			$147456$	328
$r=128$
			$36864$	12.5
$r=64$
			$9216$	0.912
$r=32$
			$2304$	0.0852
$r=16$

Equations

$$\frac{\partial }{\partial t}u(\mathbf{x},t)+\nabla \cdot \left(-D(\mathbf{x})\nabla u(\mathbf{x},t)\right)=0,\mathbf{x}\in \Omega ,$$

$$\frac{\partial }{\partial t}U(\mathbf{x},t)+\nabla \cdot \left(-{\mathbf{D}}_{\mathrm{eff}}\nabla U(\mathbf{x},t)\right)=0,\mathbf{x}\in \Omega ,$$

$$[{\mathbf{D}}_{\mathrm{eff}}{]}_{(:,1)}=\frac{1}{A}{\int }_{y_0}^{y_m}{\int }_{x_0}^{x_n}D(x,y)\nabla ({\psi }^{(x)}(x,y)+x)dxdy,$$

$$[{\mathbf{D}}_{\mathrm{eff}}{]}_{(:,2)}=\frac{1}{A}{\int }_{y_0}^{y_m}{\int }_{x_0}^{x_n}D(x,y)\nabla ({\psi }^{(y)}(x,y)+y)dxdy,$$

$$0=\nabla \cdot (D(x,y)\nabla ({\psi }^{(\xi )}(x,y)+\xi )),x_0<x<x_n,y_0<y<y_m,$$

$${\psi }^{(\xi )}(x,y)$$

$${\int }_{y_0}^{y_m}{\int }_{x_0}^{x_n}{\psi }^{(\xi )}(x,y)dxdy=0.$$

$$[{\mathbf{D}}_{\mathrm{eff}}{]}_{(:,1)}=\frac{1}{A}\sum _{i=1}^m\sum _{j=1}^n{\int }_{y_{i-1}}^{y_i}{\int }_{x_{j-1}}^{x_j}{D}_{i,j}\nabla ({\psi }_{i,j}^{(x)}(x,y)+x)dxdy,$$

$$[{\mathbf{D}}_{\mathrm{eff}}{]}_{(:,2)}=\frac{1}{A}\sum _{i=1}^m\sum _{j=1}^n{\int }_{y_{i-1}}^{y_i}{\int }_{x_{j-1}}^{x_j}{D}_{i,j}\nabla ({\psi }_{i,j}^{(y)}(x,y)+y)dxdy,$$

$$0=\nabla \cdot (D_{i,j}\nabla ({\psi }_{i,j}^{(\xi )}(x,y)+\xi )),x_{j-1}<x<x_j,y_{i-1}<y<y_i.$$

$${\psi }_{i,j}^{(\xi )}(x,y_i)={\psi }_{i+1,j}^{(\xi )}(x,y_i),$$

$$D_{i,j}\frac{\partial }{\partial y}\left({\psi }_{i,j}^{(\xi )}+\xi \right){\text{vline}}_{y=y_i}=D_{i+1,j}\frac{\partial }{\partial y}\left({\psi }_{i+1,j}^{(\xi )}+\xi \right){\text{vline}}_{y=y_i},$$

$${\psi }_{i,j}^{(\xi )}(x_j,y)={\psi }_{i,j+1}^{(\xi )}(x_j,y),$$

$$D_{i,j}\frac{\partial }{\partial x}\left({\psi }_{i,j}^{(\xi )}+\xi \right){\text{vline}}_{x=x_j}=D_{i,j+1}\frac{\partial }{\partial x}\left({\psi }_{i,j+1}^{(\xi )}+\xi \right){\text{vline}}_{x=x_j},$$

$${\psi }_{1,j}^{(\xi )}(x,y_0)={\psi }_{m,j}^{(\xi )}(x,y_m),$$

$$D_{1,j}\frac{\partial }{\partial y}\left({\psi }_{1,j}^{(\xi )}+\xi \right){\text{vline}}_{y=y_0}=D_{m,j}\frac{\partial }{\partial y}\left({\psi }_{m,j}^{(\xi )}+\xi \right){\text{vline}}_{y=y_m},$$

$${\psi }_{i,1}^{(\xi )}(x_0,y)={\psi }_{i,n}^{(\xi )}(x_n,y),$$

$$D_{i,1}\frac{\partial }{\partial x}\left({\psi }_{i,1}^{(\xi )}+\xi \right){\text{vline}}_{x=x_0}=D_{i,n}\frac{\partial }{\partial x}\left({\psi }_{i,n}^{(\xi )}+\xi \right){\text{vline}}_{x=x_n},$$

$$\frac{1}{A}\sum _{i=1}^m\sum _{j=1}^n{\int }_{y_{i-1}}^{y_i}{\int }_{x_{j-1}}^{x_j}{\psi }_{i,j}^{(\xi )}(x,y)dxdy=0.$$

$$v_{i,j}^{(\xi )}(x,y)={\psi }_{i,j}^{(\xi )}(x,y)+\xi ,$$

$$D_{i,j}\left(\frac{{\partial }^2{v}_{i,j}^{(\xi )}}{\partial x^2}+\frac{{\partial }^2{v}_{i,j}^{(\xi )}}{\partial y^2}\right)=0,$$

$$v_{i,j}^{(\xi )}(x,y_i)$$

$$v_{i,j}^{(\xi )}(x_j,y)$$

$$v_{i,1}^{(\xi )}(x_0,y)$$

$$v_{1,j}^{(\xi )}(x,y_0)$$

$$D_{i,1}\frac{\partial v_{i,1}^{(\xi )}}{\partial x}(x_0,y)=D_{i,n}\frac{\partial v_{i,n}^{(\xi )}}{\partial y}(x_n,y),i=1,\dots ,m,$$

$$D_{1,j}\frac{\partial v_{1,j}^{(\xi )}}{\partial y}(x,y_0)=D_{m,j}\frac{\partial v_{m,j}^{(\xi )}}{\partial y}(x,y_m),j=1,\dots ,n,$$

$$D_{i,j}\frac{\partial v_{i,j}^{(\xi )}}{\partial y}(x,y_i)=D_{i+1,j}\frac{\partial v_{i+1,j}^{(\xi )}}{\partial y}(x,y_i),i=1,\dots ,m-1,j=1,\dots ,n,$$

$$D_{i,j}\frac{\partial v_{i,j}^{(\xi )}}{\partial x}(x_j,y)=D_{i,j+1}\frac{\partial v_{i,j+1}^{(\xi )}}{\partial x}(x_j,y),i=1,\dots ,m,j=1,\dots ,n-1,$$

$$\sum _{i=1}^m\sum _{j=1}^n{\int }_{y_{i-1}}^{y_i}{\int }_{x_{j-1}}^{x_j}{v}_{i,j}^{(\xi )}(x,y)dxdy$$

$$D_{i,1}\frac{\partial v_{i,1}^{(\xi )}}{\partial x}(x_0,y)=D_{i,n}\frac{\partial v_{i,n}^{(\xi )}}{\partial y}(x_n,y):=g_{(i-1)n+1}(y),i=1,\dots ,m,$$

$$D_{1,j}\frac{\partial v_{1,j}^{(\xi )}}{\partial y}(x,y_0)=D_{m,j}\frac{\partial v_{m,j}^{(\xi )}}{\partial y}(x,y_m):=q_{(j-1)m+1}(x),j=1,\dots ,n,$$

$$D_{i,j}\frac{\partial v_{i,j}^{(\xi )}}{\partial x}(x_j,y)=D_{i,j+1}\frac{\partial v_{i,j+1}^{(\xi )}}{\partial x}(x_j,y):=g_{(i-1)n+j+1}(y),i=1,\dots ,m,j=1,\dots ,n-1,$$

$$D_{i,j}\frac{\partial v_{i,j}^{(\xi )}}{\partial y}(x,y_i)=D_{i+1,j}\frac{\partial v_{i+1,j}^{(\xi )}}{\partial y}(x,y_i):=q_{(j-1)m+i+1}(x),i=1,\dots ,m-1,j=1,\dots ,n,$$

$$D_{i,j}\left(\frac{{\partial }^2{v}_{i,j}^{(\xi )}}{\partial x^2}+\frac{{\partial }^2{v}_{i,j}^{(\xi )}}{\partial y^2}\right)=0,x_{j-1}<x<x_j,y_{i-1}<y<y_i,$$

$$D_{i,j}\frac{\partial v_{i,j}^{(\xi )}}{\partial x}(x_{j-1},y)=g_{(i-1)n+j}(y),D_{i,j}\frac{\partial v_{i,j}^{(\xi )}}{\partial x}(x_j,y)=g_{(i-1)n+j+1}(y),y_{i-1}<y<y_i,$$

$$D_{i,j}\frac{\partial v_{i,j}^{(\xi )}}{\partial y}(x,y_{i-1})=q_{(j-1)m+i}(x),D_{i,j}\frac{\partial v_{i,j}^{(\xi )}}{\partial y}(x,y_i)=q_{(j-1)m+i+1}(x),x_{j-1}<x<x_j,$$

$$v_{i,j}^{(\xi )}(x,y)$$

$$-h_i\sum _{k=1}^{\infty }\frac{a_{i,j,k}^{(\xi )}}{{\gamma }_{i,j,k}}\cosh \left[\frac{k\pi (x-x_j)}{h_i}\right]\cos \left(\frac{k\pi (y-y_{i-1})}{h_i}\right)$$

$$+h_i\sum _{k=1}^{\infty }\frac{b_{i,j,k}^{(\xi )}}{{\gamma }_{i,j,k}}\cosh \left[\frac{k\pi (x-x_{j-1})}{h_i}\right]\cos \left(\frac{k\pi (y-y_{i-1})}{h_i}\right)$$

$$-l_j\sum _{k=1}^{\infty }\frac{c_{i,j,k}^{(\xi )}}{{\mu }_{i,j,k}}\cosh \left[\frac{k\pi (y-y_i)}{l_j}\right]\cos \left(\frac{k\pi (x-x_{j-1})}{l_j}\right)$$

Full text

1 Introduction

At the continuum scale, diffusion in heterogeneous media is described by the following equation:

[TABLE]

where $D(\mathbf{x})>0$ is a scalar, spatially-dependent, isotropic diffusivity. When paired with appropriate initial and boundary conditions, equation (1) is infeasible to solve numerically if the scale at which the diffusivity $D(\mathbf{x})$ changes is small compared to the size of the domain $\Omega$, due to the prohibitively fine mesh required to capture the fine-scale heterogeneity [1]. One method of overcoming this problem is to homogenize equation (1) by replacing the spatially-varying diffusivity $D(\mathbf{x})$ by an effective diffusivity tensor $\mathbf{D}_{\mathrm{eff}}$, yielding the homogenized equation:

[TABLE]

where $U(\mathbf{x},t)$ is a smoothed approximation to $u(\mathbf{x},t)$. The initial and boundary conditions for the homogenized equation (2) may take a similar form as the initial and boundary conditions for the fine-scale equation (1), or they may have to be modified, as is the case in layered geometries [2]. Nevertheless, the attraction of equation (2) over equation (1) is that $\mathbf{D}_{\mathrm{eff}}$ is constant across the entire domain $\Omega$ or varies across $\Omega$ at a coarse-scale compared to $D(\mathbf{x})$. This means a coarser mesh can be used to solve the homogenized equation (2) leading to more computationally efficient simulations. The efficiency of solving the homogenized equation, however, is negated to some extent if the overhead of computing the effective diffusivity is high. For some problems, this may not be a significant issue, as the calculation of the effective diffusivities can be performed once before the simulation and stored thereafter, however the computational expense of calculating the effective diffusivities could be significant for other problems, such as those in which the geometry of $\Omega$ changes throughout the simulation. This can occur in areas such as fracturing in shale gas reservoirs [3] and cell wall collapse in wood drying [4], as each change in the geometry requires the effective diffusivity to be recomputed. To address this issue, the aim of this paper is to develop an accurate and efficient method for solving the boundary value problem required to calculate $\mathbf{D}_{\mathrm{eff}}$. While we consider two-dimensional, rectangular domains only in this work, where $\mathbf{x}=(x,y)$, the solution strategy at the core of our method can be extended to three-dimensional and circular domains111As we discuss in later sections of this paper, our semi-analytical method relies on representing the fluxes at the interfaces between adjacent blocks as unknown functions of space and integrating these functions. In three-dimensional, rectangular media double integrals would be required and in circular media these integrals would be computed using polar coordinates [5]..

According to [4], for a two-dimensional periodic domain $\Omega$ with unit cell $\mathcal{C}=[x_{0},x_{n}]\times[y_{0},y_{m}]$, the first and second columns of the corresponding effective diffusivity tensor can be computed using the following formulae:

[TABLE]

where $x_{n}>x_{0}$, $y_{m}>y_{0}$, $A=(x_{n}-x_{0})(y_{m}-y_{0})$, $\xi=x,y$ and $\psi^{(\xi)}(x,y)$ satisfies the boundary value problem:

[TABLE]

We note that the formulae for the effective diffusivity (3)–(4) and the boundary value problem (5)–(6) are equivalent to other formulations such as those presented in [6] and [7]. For a layered heterogeneous medium [8], the boundary value problem (5)–(6) has an exact solution yielding an effective diffusivity consisting of area-weighted arithmetic and harmonic averages of the layer diffusivities in the directions parallel and perpendicular to the layers, respectively [8]. In general, however, the boundary value problem cannot be solved analytically and one must turn to numerical schemes such as the finite volume method [9] or finite element method [10].

Each of these methods require the generation of a mesh over the unit cell $\mathcal{C}$ (see Figure 1(d)) that resolves the fine-scale heterogeneity, which leads to a discretization of equations (5)–(6) taking the form of a large linear system that is computationally expensive to solve. A spectral method for solving equations (5)–(6) was presented in [6], however they noted that these spectral methods are generally not well-suited to composite media with discontinuous diffusivities, such as those considered in this work, as errors caused by Gibbs phenomena can lead to inaccurate effective diffusivities.

In this paper, we develop a new semi-analytical method for solving the boundary value problem (5)–(6). Our solution method is applicable to a heterogeneous medium taking the form of a rectangular array of blocks of any size, where the diffusivity is positive, isotropic and constant within each individual block. We note that these geometries are not smooth as the diffusivity can be discontinuous across the interface between adjacent blocks. While our semi-analytical solution method is restricted to a block geometry, it can be used to homogenize complex heterogeneous geometries by pixellating and representing the geometry as a grid of blocks, as demonstrated in Figures 1(a)–(b). To calculate the effective diffusivity, our semi-analytical solution method involves reformulating the boundary value problem (5)–(6) on the heterogeneous medium into a family of boundary value problems on the homogeneous individual blocks by introducing unknown functions representing the flux at the interfaces between adjacent blocks. This strategy was recently utilised in one spatial dimension [11] and allowed accurate semi-analytical solutions to be derived for time-dependent multilayer diffusion problems involving a large number of heterogeneities (layers). In this paper, we extend these ideas to the significantly more challenging problem of steady-state diffusion (5) in two spatial dimensions where interfaces are no longer points but lines. This is achieved by solving the boundary value problems on each homogeneous block using standard techniques which yield analytical solution expressions that depend on integrals involving the unknown interface functions. By applying an appropriate numerical quadrature rule (see Figure 1(e)) to these integrals and enforcing continuity of the solution across each interface, we show how the solution can be computed in each block. The net result is a semi-analytical method for solving the well-established periodic boundary value problem (5)–(6) and hence calculating the effective diffusivity for a pixelated, locally-isotropic, heterogeneous domain.

Our new semi-analytical solution method has several desirable properties:

it is able to compute, in less computational time, a solution to the boundary value problem (5)–(6) to the same level of accuracy as a standard finite volume method (see Table 2 for details) (the reasons for choosing a finite volume method as our benchmark method are stated in section 4);

• 2.

it does not require representing the diffusivities as a function, and therefore it is not affected by the errors caused by the Gibbs phenomenon arising in standard spectral methods [6];

• 3.

it requires quadrature abscissas only on the interfaces and boundaries of the unit cell whereas standard numerical methods require a mesh across the entire unit cell (see Figure 1(d)–(e)), leading to a smaller discretized linear system (see sections 3 and 4 for details);

• 4.

it can capture complex heterogeneous geometries, as seen in Figure 1 and later in section 5.3; and

• 5.

it is an analytical method, which allows for the integrals defining the effective diffusivity (3)–(4) to be calculated analytically, whereas standard numerical methods require numerical integration to compute the effective diffusivity [9].

We stress here for clarity that our new semi-analytical method is not a new approach to homogenize heterogeneous media but an alternative for solving the boundary value problem (5)–(6) and calculating the effective diffusivity (3)–(4) [12]. Finally, we note that our method may be considered to be a form of Trefftz method according to the definition stated in [13], given that the domain $\Omega$ is partitioned into a number of smaller blocks and the solution across $\Omega$ is taken as a combination of solutions on these individual blocks. However, our method differs from Trefftz methods [13] in two key ways. Firstly, our method is fundamentally an analytical method, as we present a closed form solution to the governing PDE (5) in each block and can thus evaluate the solution $\psi^{(\xi)}$ at any point in the domain $\Omega$, whereas Trefftz methods are typically considered to be numerical methods [13]. Secondly, Trefftz methods are typically classed as a type of finite element method [13], which require the weak form of the governing PDE (5), whereas our semi-analytical method handles the strong form directly.

The remaining sections of this paper are organised as follows. In section 2 we redefine the boundary value problem (5)–(6) and effective diffusivity (3)–(4) for a block heterogeneous medium. In section 3, we develop our semi-analytical method for solving the boundary value problem (5)–(6) and computing the effective diffusivity (3)–(4), describe the implementation of the method and discuss the conditions required to ensure that it provides an accurate solution. In section 4, we discuss a standard finite volume method for solving the boundary value problem and computing the effective diffusivity. In section 5, we compare the semi-analytical and finite volume methods in terms of accuracy and efficiency for a variety of test cases from the literature, apply the semi-analytical method to the homogenization of complex heterogeneous geometries, demonstrate how our semi-analytical method can be applied to complicated curved geometries by first pixelating them and discuss extensions of our semi-analytical method to further increase its efficiency. We conclude the paper in section 6 with a summary of the key findings of the work.

2 Definition of effective diffusivity

Consider a rectangular heterogeneous domain $\mathcal{C}=[x_{0},x_{n}]\times[y_{0},y_{m}]$ consisting of an $m$ by $n$ array of blocks. Let the $(i,j)$th block have domain $[x_{j-1},x_{j}]\times[y_{i-1},y_{i}]$ and scalar diffusivity $D_{i,j}$ with vertical interfaces between blocks at $x=x_{j}$ ($j=1,\ldots,n-1$) and horizontal interfaces between blocks at $y=y_{i}$ ($i=1,\ldots,m-1$) (see Figure 1(b)). For this domain the first and second columns of the effective diffusivity $\mathbf{D}_{\mathrm{eff}}$ are given by:

[TABLE]

where $A=(x_{n}-x_{0})(y_{m}-y_{0})$ and $\psi_{i,j}^{(\xi)}$ satisfies the following PDE on the $(i,j)$th block:

[TABLE]

At the interfaces between adjacent blocks both the solution and flux are continuous [8], leading to the following interface conditions:

[TABLE]

for $x_{j-1}<x<x_{j}$, $i=1,\ldots,m-1$ and $j=1,\ldots,n$ and

[TABLE]

for $y_{i-1}<y<y_{i}$, $i=1,\ldots,m$ and $j=1,\ldots,n-1$, while periodic conditions hold at the external boundaries of the domain:

[TABLE]

for $x_{j-1}<x<x_{j}$ and $j=1,\ldots,n$ and

[TABLE]

for $y_{i-1}<y<y_{i}$ and $i=1,\ldots,m.$ The boundary value problem (10)–(18) has an infinite number of solutions [9] as any arbitrary constant may be added to the solution $\psi_{i,j}^{(\xi)}$ and it will still satisfy the PDE (10) as well as the boundary conditions (11)–(18). The zero mean condition (7) can be rewritten in the following form:

[TABLE]

With $\psi_{i,j}^{(\xi)}(x,y)$ known for $\xi=x$ and $\xi=y$, the effective diffusivity, $\mathbf{D}_{\mathrm{eff}}$, can be computed from equations (8)–(9). Interpretation of the entries of the symmetric positive definite matrix $\mathbf{D}_{\mathrm{eff}}$ comes from its diagonalization, $\mathbf{D}_{\mathrm{eff}}=\mathbf{P}\boldsymbol{\Lambda}\mathbf{P}^{T}$. Here, $\mathbf{P}\in\mathbb{R}^{2\times 2}$ is an orthonormal matrix whose columns are unit vectors pointing in the principal directions of diffusion and $\boldsymbol{\Lambda}\in\mathbb{R}^{2\times 2}$ is a diagonal matrix whose entries represent the apparent diffusivities in these directions [14].

3 Semi-analytical method

3.1 Solution approach

We now describe our new semi-analytical solution approach for solving the boundary value problem (10)–(19) which enables the effective diffusivity (8)-(9) to be determined. We first transform equation (10) into Laplace’s equation via the transformation:

[TABLE]

for $i=1,\ldots,m$ and $j=1,\ldots,n$. Combining equation (20) with equations (10) yields:

[TABLE]

for $x_{j-1}<x<x_{j}$ and $y_{i-1}<y<y_{i}$ and combining equation (20) with equations (11), (13), (15) and (17) yields the following conditions on the solution $v_{i,j}^{(\xi)}$ at the interfaces and boundaries:

[TABLE]

for $x_{j-1}<x<x_{j}$ and $y_{i-1}<y<y_{i}$ and combining equation (20) with equations (12), (14), (16) and (18) yields the following conditions on the flux at the interfaces and boundaries:

[TABLE]

for $x_{j-1}<x<x_{j}$ and $y_{i-1}<y<y_{i}$. Combining the zero mean condition (19) with the transformation (20) yields:

[TABLE]

Our semi-analytical approach for solving this transformed boundary value problem involves first setting the fluxes at the interfaces and boundaries in equations (26)–(29) to be unknown spatial functions [15] as follows:

[TABLE]

for $x_{j-1}<x<x_{j}$ and $y_{i-1}<y<y_{i}$. This allows us to reformulate the heterogeneous boundary value problem (21)–(29) into a family of boundary value problems on each of the $mn$ blocks (see Figure 2):

[TABLE]

where the unknown interface and boundary functions, introduced in equations (31)–(34), are constrained to satisfy the interface and boundary conditions (22)–(25). The boundary conditions (36)–(37) are valid for all blocks except those in either the bottom row ($i=m$) or right column $(j=n)$. The boundary conditions for these blocks are similar to equations (36)–(37) and detailed in A.

The solution of the boundary value problem (35)–(37) is given by [16]222For three-dimensional problems, the equivalent solution would be more complicated, requiring a double summation and the coefficients to be computed in terms of double integrals.:

[TABLE]

where:

[TABLE]

and $K_{i,j}$ is an arbitrary constant. While for the solution of the boundary value problem (35)–(37) the constant $K_{i,j}$ is arbitrary, in the context of the solution across the entire unit cell $\mathcal{C}$ consisting of all $mn$ boundary value problems, the constants are not arbitrary as they are required to ensure continuity of the solution across the interfaces and boundaries. Note that the coefficients $a_{i,j,k}^{(\xi)}$, $b_{i,j,k}^{(\xi)}$, $c_{i,j,k}^{(\xi)}$ and $d_{i,j,k}^{(\xi)}$ appearing in equation (38) depend on the values of $i$ and $j$. For $i=1,\ldots,m-1$ and $j=1,\ldots,n-1$:

[TABLE]

where $k=0,1,\ldots,\infty$, while for the bottom row ($i=m$) and right column ($j=n$) of blocks the coefficients are similar and given in B.

Due to the boundary conditions (36)–(37), the solution (38) requires an additional solvability condition to ensure a net flux of zero is imposed on each block [16]. This condition takes the form:

[TABLE]

for $i=1,\ldots,m-1$ and $j=1,\ldots,n-1$. The solvability conditions for the bottom row ($i=m$) and right column ($j=n$) of blocks have a similar form to equation (45) and are given in C.

Through the coefficients (41)–(44), the solution (38) in each block involves the integral of the unknown interface and boundary functions, introduced in equations (31)–(34), and the unknown constants $K_{i,j}$. By applying the as yet unused interface and boundary conditions (22)–(25), we now derive a linear system of equations whose solution allows the integrals appearing in the coefficients (41)–(44) to be approximated. This is achieved by first applying a quadrature rule to these integrals. For example using the expressions for $a_{i,j,k}^{(\xi)}$ and $c_{i,j,k}^{(\xi)}$ in equations (41)–(44), we have:

[TABLE]

[TABLE]

where $w_{x}^{(j,p)}$ and $w_{y}^{(i,p)}$ are the weights of the $p$th abscissas, $x^{(j,p)}$ and $y^{(i,p)}$, respectively and $N_{x}$ and $N_{y}$ denote the number of abscissas. The weights $w_{x}^{(j,p)}$ and $w_{y}^{(i,p)}$ and abscissas $x^{(j,p)}$ and $y^{(i,p)}$ depend on the type of quadrature rule chosen, as will be discussed in the following section. There are a total of $2mn$ unknown flux functions, $q_{(j-1)m+i}(x)$ and $g_{(i-1)n+j}(y)$ for $i=1,\ldots,m$ and $j=1,\ldots,n$, with the integrals involving these functions approximated using a quadrature rule with $N_{x}$ and $N_{y}$ function evaluations, respectively. In total, this yields $mn(N_{x}+N_{y})$ unknown evaluations at the abscissas: $q_{(j-1)m+i}(x^{(j,p)})$ for $i=1,\ldots,m$, $j=1,\ldots,n$ and $p=1,\ldots,N_{x}$ and $g_{(i-1)n+j}(y^{(i,p)})$ for $i=1,\ldots,m$, $j=1,\ldots,n$ and $p=1,\ldots,N_{y}$. To compute these evaluations, which appear in equations (46)–(47), we make use of the boundary and interface conditions (22)–(25) evaluated at each of the abscissas:

[TABLE]

for $x_{j-1}<x^{(j,p)}<x_{j}$ and $p=1,\ldots,N_{x}$,

[TABLE]

for $y_{i-1}<y^{(i,p)}<y_{i}$ and $p=1,\ldots,N_{y}$. Equations (48)–(51) define a system of $mn(N_{x}+N_{y})$ equations linear in the $mn(N_{x}+N_{y})$ unknown evaluations at the abscissas. In addition to these unknowns, we have an additional $mn$ unknowns corresponding to the constants $K_{i,j}$ appearing in equation (38). These additional unknowns are accounted for by imposing the solvability conditions (45), with the integrals replaced with quadrature approximations similar to those used in equations (46)–(47). The solvability conditions ensure a zero net flux on each block and as periodic conditions are imposed on the external boundaries of the unit cell $\mathcal{C}$, the solvability condition in any one of the $mn$ blocks is redundant. To account for this redundancy, we replace the solvability conditions on block $mn$ with the zero mean condition (19). In the solution expression (38), all terms with the exception of the constants $K_{i,j}$ integrate to zero over their domain of integration $[x_{j-1},x_{j}]\times[y_{i-1},y_{i}]$. This allows for the zero mean condition (30) to be written in the following form:

[TABLE]

Assembling equations (48)–(52) yields a linear system:

[TABLE]

where $\mathbf{A}_{\text{S}}\in\mathbb{R}^{N\times N}$, $\mathbf{x}_{\text{S}}\in\mathbb{R}^{N}$, $\mathbf{b}_{\text{S}}\in\mathbb{R}^{N}$ and $N=mn(N_{x}+N_{y}+1)$. The entries of $\mathbf{A}_{\text{S}}$ and $\mathbf{b}_{\text{S}}$ are determined via the solvability conditions (45), the continuity conditions (48)–(51) and the zero mean condition (52). The vector $\mathbf{x}_{\text{S}}$ contains the function evaluations $q_{(j-1)m+i}(x^{(j,p)})$ for $i=1,\ldots,m$, $j=1,\ldots,n$ and $p=1,\ldots,N_{x}$ and $g_{(i-1)n+j}(y^{(i,p)})$ for $i=1,\ldots,m$, $j=1,\ldots,n$ and $p=1,\ldots,N_{y}$ as well as the constants $K_{i,j}$ for $i=1,\ldots,m$ and $j=1,\ldots,n$. Note that when constructing the linear system, we truncate the summations appearing in the solution expression (38) after a finite number of terms, $k=N_{\text{eig}}$. Solving the linear system (53) then allows the approximations to the integrals (46) and (47) to be computed and thus the solution (38) can be evaluated. With $v_{i,j}^{(\xi)}(x,y)$ computed, the solution of the original boundary value problem (10)–(19) can be calculated as $\psi^{(\xi)}_{i,j}(x,y)=v_{i,j}^{(\xi)}(x,y)-\xi$.

In terms of the transformed solution $v_{i,j}^{(\xi)}(x,y)$, the four entries of the effective diffusivity are given by:

[TABLE]

As we have an analytical expression for $v_{i,j}^{(\xi)}(x,y)$ (38), we can evaluate the derivatives and integrals appearing in equations (54)–(57) directly, yielding:

[TABLE]

where $A_{i,j}=(x_{j}-x_{j-1})(y_{i}-y_{i-1})$. Finally, we remark that the coefficient matrix $\mathbf{A}_{\text{S}}$ in the linear system (53) is identical for both $\xi=x$ and $\xi=y$ , however the vector $\mathbf{b}_{\text{S}}$ differs. This means that the linear system for $\xi=x$ and the linear system for $\xi=y$ can both be solved with only one matrix factorisation.

3.2 Choice of quadrature method

There are several constraints governing the choice of quadrature method utilized in the approximations of the integrals (46)–(47):

*Abscissas cannot include the vertices of blocks.

*If a vertex of a block is used as an abscissa in the quadrature rules (46)–(47), then the set of four continuity conditions (48)–(51) evaluated at the vertex will be linearly dependent. If all $mn$ vertices are used in the quadrature rules (46)–(47), then the linear system (53) is deficient in rank by $mn$. As the vertices of the blocks are the limits of integration for all integrals appearing in the coefficients (41)–(44), to ensure linear independence we implement a quadrature rule that uses at most one of the limits of integration as an abscissa. While we could implement different quadrature rules that use endpoints for some integrals and do not use endpoints for other integrals, we take the consistent option of using no endpoints.

• 2.

*The quadrature rule must use the same abscissas for all frequencies $k=0,\ldots,N_{\text{eig}}$ appearing in the coefficients (41)–(44).

*The linear system (53) is generated by evaluating the boundary and interface conditions (48)–(51) at the abscissas used in the quadrature formulae (46)-(47). In order to minimise the size of the the linear system (53), the same abscissas $x^{(j,p)}$ and $y^{(i,p)}$ should be used for all $k=0,\ldots,N_{\text{eig}}$. Hence, we restrict the choice of quadrature rules to those that use the same abscissas for all frequencies $k=0,\ldots,N_{\text{eig}}$.

• 3.

*The quadrature rule must be able to compute the integrals appearing in equations (41)–(44) for all frequencies $k=0,\ldots,N_{\text{eig}}$ for $N_{\text{eig}}=\text{max}(N_{x},N_{y})-1$ to a sufficient level of accuracy.

*The boundary and interface conditions (48)–(51) result in the coefficient matrix $\mathbf{A}_{S}$ in the linear system (53) being rank deficient if $N_{\text{eig}}=\text{max}(N_{x},N_{y})-1$. As there is a lower limit on $N_{\text{eig}}$, the quadrature rule must be able to accurately compute the integrals appearing in equations (41)–(44) for all frequencies $k=0,\ldots,N_{\text{eig}}$, for $N_{\text{eig}}=\text{max}(N_{x},N_{y})-1$.

• 4.

*The quadrature rule should be more accurate for low frequencies than for high frequencies.

*The entries in the coefficient matrix $\mathbf{A}_{S}$ appearing in the linear system (53) are generated through a summation from $k=0,\ldots,N_{\text{eig}}$ of evaluations of the solution (38) at the abscissas used in the quadrature rule. As $k$ increases, these solution evaluations decrease and therefore solution evaluations with low frequencies comprise a larger proportion of the entries of the coefficient matrix $\mathbf{A}_{S}$ than solution evaluations with higher frequencies. Therefore, the quadrature rule should be more accurate for low frequencies than for high frequencies.

• 5.

The quadrature rule must be a linear function of the function evaluations and cannot make use of any information about the function being integrated.

As the evaluations, $q_{(j-1)m+i}(x^{(j,p)})$ for $i=1,\ldots,m$, $j=1,\ldots,n$ and $p=1,\ldots,N_{x}$ and $g_{(i-1)n+j}(y^{(i,p)})$ for $i=1,\ldots,m$, $j=1,\ldots,n$ and $p=1,\ldots,N_{y}$, used to approximate the integrals (41)–(44) are determined via the solution of the linear system (53), the quadrature rule must be a linear function of these evaluations and cannot use any other information about the function, such as its derivatives, which are commonly used in approximations of oscillatory integrals [17].

As the integrals appearing in the coefficients (41)–(44) contain an oscillatory function, we investigated Filon [18] and Levin [19] type methods, which are quadrature rules well-suited to highly oscillatory functions. The choice of quadrature rule is further complicated by the integrands being unknown. We investigated a wide variety of different Filon and Levin type quadrature methods, including (but not limited to) the many different methods presented in [20] and found that almost all of them were not suitable for use in our semi-analytical method, as they do not satisfy one or more of the above mentioned constraints. The quadrature methods that were found to be suitable are: the midpoint rule, a Filon-midpoint rule demonstrated in [21] and various Gaussian quadrature methods that do not use abscissas at the endpoints. The Filon-midpoint rule was found to be consistently less accurate than the regular midpoint rule in testing its application to oscillatory integrals. We then tested our semi-analytical method with both the midpoint and Gauss-Legendre quadrature rules for a variety of different problems. For the midpoint rule, we found that the matrix $\mathbf{A}_{\text{S}}$ is of full rank for every problem as long as the restriction $N_{\text{eig}}\geq\text{max}(N_{x},N_{y})-1$ is satisfied. However, for the Gauss-Legendre quadrature rule, we found that for some problems, the matrix $\mathbf{A}_{\text{S}}$ is rank deficient despite this inequality being satisfied. We found that setting the value of $N_{\text{eig}}$ to be sufficiently large yields a matrix $\mathbf{A}_{\text{S}}$ of full rank. However, increasing the value of $N_{\text{eig}}$ is undesirable as the Gauss-Legendre rule can be inaccurate for highly oscillatory functions. To summarise, on the basis of our findings we decided not to use any form of Gaussian quadrature to approximate the integrals (46)–(47) as we cannot be sure that it will be accurate in our approach, bearing in mind that these approximations are key to the success of the semi-analytical method. We therefore choose the midpoint rule to approximate the integrals (46)–(47). We note that more suitable quadrature methods may exist and/or could be developed, and if these methods can be implemented in the form of summations as presented in equations (46)–(47), then they could be easily applied to our semi-analytical method.

The midpoint rule uses evenly spaced abscissas, which leads to a restriction on the maximum number of terms, $N_{\text{eig}}$, taken in the summations appearing in the solution (38). The functions integrated in the definition of the coefficients (41)–(44) are similar to $\widetilde{f}(x)=f(x)\cos\left(k\pi x/L\right)$ with limits of integration of $x=0$ and $x=L$. The period of $\widetilde{f}(x)$ is $2L/k$ and with $N_{x}$ evenly spaced abscissas used in the approximation of the integrals, the spacing between abscissas is $L/(N_{x}-1)$. If $2L/k=L/(N_{x}-1)$ the abscissas align with the peaks of $\widetilde{f}(x)$, which causes each evaluation of $\widetilde{f}(x)$ to occur at an abscissa $x=x^{*}$ at which $\widetilde{f}(x^{*})=f(x^{*})$. This has the effect of returning the value of the integral of $f(x)$ from $x=0$ to $x=L$ instead of $\widetilde{f}(x)$. Hence, to ensure an accurate solution we set an upper limit on the maximum number of terms in the summations in equation (38) to ensure $2L/k<L/(N_{x}-1)$, for all $k=1,\ldots,N_{\text{eig}}$, leading to the restriction $N_{\text{eig}}\leq 2N_{x}-3$. Combining both the lower and upper bounds on $N_{\text{eig}}$ yields the restriction:

[TABLE]

We note that the lower bound in the restriction (62) applies to any quadrature rule, whereas the upper bound applies specifically to the midpoint rule. It also follows that $N_{x}\geq 2$ must hold.

4 Finite volume method

To verify the accuracy of our new semi-analytical method, we compare it to a standard vertex-centered finite volume method [9]. We first define a rectangular mesh over the unit cell $\mathcal{C}$ consisting of uniformly spaced nodes in the $x$ and $y$ directions. The nodes are located at $x=x_{0}+kh_{x}$ for $k=0,\ldots,N_{\text{F}}^{(x)}$ and $y=y_{0}+kh_{y}$ for $k=0,\ldots,N_{\text{F}}^{(y)}$, where $h_{x}=(x_{n}-x_{0})/N_{\text{F}}^{(x)}$ and $h_{y}=(y_{m}-y_{0})/N_{\text{F}}^{(y)}$. Control volumes are formed around each node by connecting the centroid of each rectangular element to the midpoint of its edges. For interior nodes this gives a rectangular control volume comprised of eight “half-sized” edges, each confined to a single element of homogeneous material. The number of nodes in the $x$ and $y$ directions are always chosen to ensure that nodes coincide with interfaces between blocks so that each control volume edge is located within only one block.

To perform the finite volume discretization, the PDE (5) is integrated over each control volume and the integrals along each control volume edge approximated using a midpoint rule. This procedure leads to a semi-discretized finite volume equation of the form:

[TABLE]

where $e_{x}$ and $e_{y}$ are the sets of vertical and horizontal edges of the control volume, $(x_{\text{mid}}^{(s)},y_{\text{mid}}^{(s)})$ are the coordinates of the midpoint of edge $s$ of the control volume and $\delta_{x\xi}$ and $\delta_{y\xi}$ are Kronecker deltas, taking the value of one when the subscripts are equal and zero otherwise. The equation (63) is valid for all nodes, except those located along the boundaries $x=x_{n}$ and $y=y_{m}$. For these nodes, the finite volume equation (63) is modified to accommodate the periodic boundary conditions, as discussed in detail in [4].

Approximating the spatial derivatives appearing in equation (63) using second-order central differences and assembling the resulting finite volume equations for each node yields a linear system of the form:

[TABLE]

where $\mathbf{x_{\text{F}}}$ is a vector of length ${N_{\text{F}}^{(x)}N_{\text{F}}^{(y)}}$ containing the discrete unknown values of $\psi^{(\xi)}$ at each of the nodes. After solving this linear system, we use a bilinear interpolant of $\psi^{(\xi)}$ over each element to approximate the gradient within each element. This allows us to approximate the integrals in equation (3)–(4), yielding the following formulae for the first and second columns of the effective diffusivity:

[TABLE]

where $(x_{c}^{(p)},y_{c}^{(q)})$ is the centroid of the $(p,q)$th element, $\mathbf{e}_{1}=[1,0]^{T}$ and $\mathbf{e}_{2}=[0,1]^{T}$.

5 Results and discussion

In this section we verify the accuracy and efficiency of our new semi-analytical method and demonstrate its applicability to computing effective diffusivities of complex heterogeneous geometries. In all comparisons between the semi-analytical and finite volume methods, we ensure that the spacing between abscissas and the spacing between nodes are identical by setting $N_{\text{F}}^{(x)}=nN_{x}$ and $N_{\text{F}}^{(y)}=mN_{y}$ in the finite volume method, where $N_{x}$ and $N_{y}$ denote the number of abscissas per interface in the semi-analytical method. Both the semi-analytical and finite volume methods were implemented in MATLAB and the linear systems for both methods, equations (53) and (64), were stored in sparse format and solved using the backslash operator.333For both the semi-analytical and finite volume methods, the backslash operator implements the Unsymmetric MultiFrontal PACKage with automatic reordering. All numerical experiments were carried out in MATLAB (version 2017A) running on an early 2015 MacBook Pro with a 2.7GHz dual-core Intel Core i5 processor and 8GB of memory.

The efficiency of the MATLAB code implementing the semi-analytical method can be improved by assuming that $m=n$ (i.e. the same number of blocks in both the horizontal and vertical directions) and that all $mn$ blocks are the same size. The main source of the improvement in efficiency is in being able to compute the summations appearing in the solution (38) only once, instead of once for each abscissa used in the method. All of the test cases we consider and the majority of the problems we encountered in the literature have $m=n$ and identically sized blocks, so all calculations in this paper are computed using the more efficient version of the MATLAB code. The MATLAB codes implementing the semi-analytical method for both $m=n$ and $m\neq n$ are available on our GitHub repository: https://github.com/NathanMarch/Homogenization.

5.1 Accuracy and efficiency of the semi-analytical method

As a first verification of the accuracy of the semi-analytical solution we consider four geometries previously presented in [7]. These four test geometries are depicted in Table 1 and take the form of a square cell of dimension $1$ by $1$ comprised of an $8$ by $8$ array of blocks each of dimension $0.125$ by $0.125$. The dark grey blocks have a diffusivity of $0.1$ and the light grey blocks have a diffusivity of $1$. Case 1 has a square shaped inclusion of dark grey blocks in the centre of the cell, case 2 consists of one larger dark grey inclusion in the centre of the cell and four smaller dark grey inclusions at the corners, case 3 consists of three layers, of which the middle layer is dark grey and case 4 contains an L-shaped dark grey inclusion. We compute the effective diffusivity using both the semi-analytical and finite volume methods for cases 1–4 and compare the results to those presented in [7]. We set $N_{x}=N_{y}=16$ and $N_{\text{eig}}=2N_{x}-3=29$ for the semi-analytical method and $N_{\text{F}}^{(x)}=N_{\text{F}}^{(y)}=4N_{x}=4N_{y}=64$ for the finite volume method. The results in Table 1 demonstrate that both the semi-analytical and finite volume methods are in very good agreement with the results presented by [7], agreeing to at least two significant figures for all entries in each effective diffusivity.

The effective diffusivities for cases 1–3 are diagonal matrices as each geometry is invariant under a rotation of $180^{\circ}$ [22]. Additionally, the diagonal entries of the effective diffusivity for case 1 are identical as the geometry is also invariant under a rotation of $90^{\circ}$. As cases 1–3 have effective diffusivities that are diagonal, the principal directions of diffusion are in the horizontal ($x)$ and vertical ($y$) directions and the diagonal elements represent the diffusivity in these directions (see Table 1). For cases 2 and 3 the diffusivity in the horizontal direction is larger than the diffusivity in the vertical direction. This is because, unlike for flow in the horizontal direction, flow in the vertical direction must pass through or around the low diffusivity (dark grey) blocks. Case 3 is a layered medium, so the diffusivity in the horizontal and vertical directions are the area-weighted arithmetic and harmonic averages of the layer diffusivities, respectively. The eigenvalues of the effective diffusivity for case 4 are $0.527$ and $0.682$ and the corresponding normalised eigenvectors are $[0.982,0.189]^{T}$ and $[0.189,-0.982]^{T}$, which define principal directions of diffusion for case 4 that are rotated $10.9^{\circ}$ anti-clockwise from the standard Cartesian axes (see Table 1).

We now compare the accuracy of the semi-analytical method to the finite volume method by first computing a benchmark effective diffusivity using the finite volume method with a very fine node spacing. We consider the geometries of cases 1 and 4 as presented in [7], represented as both $4$ by $4$ and $16$ by $16$ grids of equal-sized blocks. We label theses cases as cases A–D such that:

Case A is the $4$ by $4$ configuration of Case 1;

• 2.

Case B is the $4$ by $4$ configuration of Case 4;

• 3.

Case C is the $16$ by $16$ configuration of Case 1;

• 4.

Case D is the $16$ by $16$ configuration of Case 4.

For cases A and B each block is of size $0.25$ by $0.25$ while for the $16$ by $16$ configuration each block is of size $0.0625$ by $0.0625$. For cases A and B we compute the effective diffusivity using $N_{x}=N_{y}=4,8,16,32$ and $64$ for the semi-analytical method and $N_{\text{F}}^{(x)}=N_{\text{F}}^{(y)}=4N_{x}=4N_{y}=16,32,64,128$ and $256$ for the finite volume method. For cases C and D we compute the effective diffusivity using $N_{x}=N_{y}=4,8$ and $16$ for the semi-analytical method and $N_{\text{F}}^{(x)}=N_{\text{F}}^{(y)}=16N_{x}=16N_{y}=64,128$ and $256$ for the finite volume method. For all cases, the benchmark effective diffusivity is computed using the finite volume method with $N_{\text{F}}^{(x)}=N_{\text{F}}^{(y)}=1024$. For the semi-analytical method, $N_{\text{eig}}$ is set as the minimum of $2N_{x}-3$ and $100$. This is because, as mentioned in section 3.2, the value $2N_{x}-3$ is the upper limit on the maximum number of terms in the summations in equation (38) and in preliminary testing we found that using more than $100$ terms in the summation had negligible effect on the solution. We compute our relative error $\mathbf{E}=|(\mathbf{D}_{\mathrm{eff}}-\widehat{\mathbf{D}}_{\text{eff}})\varoslash\widehat{\mathbf{D}}_{\text{eff}}|$, where $\varoslash$ represents Hadamard/element-wise division, $\mathbf{D}_{\mathrm{eff}}$ represents the effective diffusivity as calculated from the semi-analytical/finite volume solution under comparison and $\widehat{\mathbf{D}}_{\text{eff}}$ is the benchmark effective diffusivity. We also record the runtime for both the semi-analytical and finite volume methods, which includes the time taken to solve the boundary value problem (21)–(30) and compute the effective diffusivity via equations (58)–(61) and (65)–(66). All runtimes reported in this section were calculated by performing each calculation ten times and taking the median runtime. Results in Table 2 demonstrate that the semi-analytical method is either faster or more accurate (or both) than the finite volume method, indicating that for a specified level of accuracy the semi-analytical method is more efficient.

We now further explore the efficiency of our new semi-analytical method by comparing the size of the linear systems, equations (53) and (64), and the time taken to compute the effective diffusivity for both the semi-analytical and finite volume methods. To do this we consider a square unit cell with $x_{0}=y_{0}=0$ and $x_{n}=y_{m}=1$, consisting of an $m$ by $m$ grid of equal-sized blocks for $m=2,\ldots,40$, i.e., we consider a $2$ by $2$ grid, a $3$ by $3$ grid etc. The unit cell is a checkerboard structure consisting of diffusivities of $0.1$ and $1$ with a diffusivity of $0.1$ in the top left corner. For example, the diffusivities for $m=3$ are:

[TABLE]

As usual, we specify the spacing between abscissas for the semi-analytical method to be the same as the spacing between nodes for the finite volume method. We compute the effective diffusivity using $N_{x}=N_{y}=16$ and $N_{\text{eig}}=2N_{x}-3=29$ for the semi-analytical method and $N_{\text{F}}^{(x)}=N_{\text{F}}^{(y)}=16m$ for the finite volume method. This leads to the linear system (53) in the semi-analytical method having dimension $m^{2}(2N_{x}+1)=33m^{2}$ and the linear system (64) in the finite volume method having dimension $m^{2}N_{x}^{2}=256m^{2}$. Results in Figure 3 demonstrate that the efficiency advantage of the semi-analytical method over the finite volume method becomes more pronounced as $m$ increases. This is because as the number of blocks increases, the size of the linear system increases (see Figure 3) and the time taken to solve these systems becomes the dominant computational cost of calculating the effective diffusivity. This improvement in efficiency that the semi-analytical method has over the finite volume method demonstrates the potential of the semi-analytical method to significantly speed up simulations of transport processes across block, heterogeneous media.

5.2 Accuracy of the solution of the linear system (53)

In this section we demonstrate the accuracy of the solution of the linear system (53) corresponding to cases 1–4 with $N_{x}=N_{y}=16$ and all values of $N_{\text{eig}}$ that satisfy the restriction (62), which are $N_{\text{eig}}=15,\ldots,29$. We stress that it is important that the lower bound of the restriction (62) is satisfied, as otherwise the coefficient matrix $\mathbf{A}_{s}$ will be singular. For each case, the condition number (as calculated using the $2$-norm) of the coefficient matrix $\mathbf{A}_{s}$ is of the order $10^{3}$ and the $2$-norm of the residual of the linear system (53) is of the order $10^{-14}$, which means the continuity conditions (48)–(51) are satisfied with a maximum absolute error of $10^{-14}$. An upper bound on the relative error in the solution of the linear system (53) is $\text{cond}(\mathbf{A}_{s})\left\lVert\mathbf{r}_{s}\right\rVert_{2}/\left\lVert\mathbf{b}_{s}\right\rVert_{2}$, where $\mathbf{r}_{s}$ is the residual of the linear system (53), which yields a maximum relative error of $10^{-11}/\left\lVert\mathbf{b}_{s}\right\rVert_{2}$.

5.3 Results for complex structured heterogeneous geometries

We now demonstrate how the semi-analytical method can be used to homogenize complex heterogeneous geometries. For our demonstration, we consider a square cell of dimension $1$ by $1$ consisting of an $m$ by $m$ array of square blocks and generate aggregated, random geometries of two different media. We first generate $m^{2}$ uniformly distributed, random numbers between [math] and $1$ denoted by $D^{(0)}_{i,j}$ where $i=1,\ldots,m$ and $j=1,\ldots,m$. We then perform the following iteration:

[TABLE]

for all $i=1,\ldots,m$ and $j=1,\ldots,m$, where $D_{i,j}^{(k)}$ represents the diffusivity in the $(i,j)$th block of the unit cell $\mathcal{C}$ (see Figure 1(b)) at the $k$th iteration and the weightings are inspired by those used in lattice Boltzmann methods [24]. For the purposes of the algorithm in equation (68), we assume periodicity to process the blocks along the boundaries, for example, when $i=1$, we set $i-1$ to be equal to $m$. After performing a fixed number of iterations, we are left with an array of diffusivities containing entries between [math] and $1$, where larger values are aggregated together and smaller values are aggregated together. The reason for this aggregation is that during each iteration, the diffusivity of a given block is a weighted average of itself and its nearest neighbours, so in each iteration, the diffusivities are converging to those of its nearest neighbours. To ensure the geometry is comprised of equal parts light grey (diffusivity of 1) and dark grey (diffusivity of 0.1), we then prescribe that the largest $0.5m^{2}$ entries of the array are 1 and the remaining $0.5m^{2}$ entries are 0.1.

Table 3 presents some example geometries generated using the above procedure for $m=10,20,50$ and $100$, which resemble realistic geometries presented in [25] and [22]. These figures demonstrate the highly complex heterogeneous geometries that can be captured using blocks and how the semi-analytical method can be used to calculate the effective diffusivity for such geometries. We note that each iteration in the aggregation procedure could be considered as a change in the underlying geometry of a domain arising throughout the simulation of a diffusive transport problem. In this scenario, the effective diffusivity can be calculated easily and efficiently using our semi-analytical method, as it requires merely changing the diffusivities in some of the blocks. We also note that our semi-analytical method can be applied to any geometry where the diffusivities are represented as a grid of pixels (e.g. CT scans).

5.4 Results for complex unstructured heterogeneous geometries

We now demonstrate the applicability of the semi-analytical method to geometries involving curved edges by first pixelating them, which allows for the geometry to be represented as a grid of blocks. We consider a geometry presented in [9] (see Table 4(a)) and pixelate it several times to generate geometries of different sizes (see Figure 1 (b)). This is a $512$ by $512$ geometry in which the dark grey blocks have a diffusivity of $0.1$ and the light grey blocks have a diffusivity of $1$. We define a pixelation parameter $r$, such that the image is now represented as an $r$ by $r$ grid of square blocks, so decreasing values of $r$ produce greater pixelation. Each block in the pixelated geometry is made up of a $k$ by $k$ grid of blocks, such that $rk=m$ (for the geometry in Figure 1(a), $m=512$). On each $k$ by $k$ grid of blocks, we compute the mean of the diffusivities of the $k^{2}$ blocks. If the mean is greater than or equal to $0.55$, the diffusivity of the entire $k$ by $k$ block is set to be $1$, otherwise it is set to be $0.1$. The threshold value of $0.55$ is chosen as it is halfway between the low diffusivity of $0.1$ and the high diffusivity of $1$ and as the pixellation is performed via a simple averaging process, the computational cost is negligible. On each of the pixelated geometries, we compute the effective diffusivity using the semi-analytical method with $N_{x}=N_{y}=4$ abscissas per interface and compare the result to a benchmark effective diffusivity obtained from the finite volume method of section 4 applied to the original 512 by 512 geometry with a fine mesh ($N^{(x)}_{\text{F}}$ = $N^{(y)}_{\text{F}}$ = $1024$). We consider levels of pixelation corresponding to $r=16,32,64,128$ and for each level of pixelation compute the relative error $\mathbf{E}$ defined in section 5.1.

Results in Table 4 demonstrate that as the pixelation increases (as $r$ decreases), the error in the effective diffusivity generally increases, that is, higher values of $r$ generally yield more accurate effective diffusivities. These results demonstrate the applicability of pixelating complicated geometries to the fast computation of accurate effective diffusivities. For example, pixelating the image to a $32$ by $32$ geometry allows the diagonal entries of the effective diffusivity to be calculated using our semi-analytical method to an accuracy of approximately $1\%$ in under $1$ second. Although a standard finite volume method could also be applied to the pixelated $32$ by $32$ geometry, our semi-analytical method is faster and of equivalent accuracy, as demonstrated in Table 2.

5.5 Extension of semi-analytical method to geometries without aligned interfaces

Although in this paper we consider geometries comprised of square blocks in which all interfaces are aligned (such as in Figure 1(b)), our semi-analytical method can be extended to non-aligned rectangular blocks where interfaces between adjacent blocks of homogeneous material are removed wherever possible. For example, in Figure 4(b), we consider the pixelated geometry corresponding to $r=16$ in Table 4 and divide it into a number of rectangles of different sizes. The rectangles are chosen in such a way as to attempt to minimise the sum of the perimeter of the rectangles. The abscissas used for the numerical integration of the flux functions in equations (46)–(47) are placed on the interfaces between adjacent rectangles. Due to the fewer interfaces in Figure 4(b), the total number of abscissas required is reduced, which would lead to decreased computational time, due to the smaller size of the linear system (53) and increased accuracy, due to the fewer integrals along interfaces (46)–(47) required to be approximated. For example in Figure 4(a), we show the abscissas required for the semi-analytical method for the pixelated $16$ by $16$ geometry from Table 4 with $N_{x}=N_{y}=4$, which requires the solution of a $2304$ by $2304$ linear system in order to calculate the effective diffusivity. The linear system corresponding to the formulation in Figure 4(b), however is only $803$ by $803$, where $764$ of the unknowns in the linear system correspond to the abscissas of the flux function evaluations on the interfaces and the remaining $39$ unknowns correspond to the constants $K_{i,j}$ that appear in equation (38). This extended semi-analytical method reduces the size of the linear system by over $65\%$ for the example in Figure 4 and highlights its potential to greatly reduce the computational time required to calculate effective diffusivities. In this sense, the method could be considered as a form of adaptive approach, as the abscissas used in the calculation of the effective diffusivities are chosen based on the diffusivities of the individual blocks and may change if the underlying geometry changes throughout the course of a simulation. The extended semi-analytical method shares some of the advantages and disadvantages of numerical solution methods that use adaptive meshes that fit the geometry of the unit cell, such as the methods presented in [26] and [27]. All of these methods can be more computationally efficient due to fewer unknowns required to represent the solution $\psi^{(\xi)}$ in sections of the unit cell where the geometry is unchanged. However pre-processing is required to generate an adaptive mesh and remeshing may be required if the geometry changes throughout a simulation.

5.6 Extension of semi-analytical method to discontinuous media

Although in this work we consider positive block diffusivities, our method can be extended to discontinuous media in which some of the blocks have zero diffusivity. To account for this, on any interface between two blocks in which at least one of the blocks has zero diffusivity, the interface and boundary functions (36)–(37) would be set to equal zero along that interface. Additionally, no boundary value problem (10) is imposed on any block in which the diffusivity is zero, as equation (10) would be trivially satisfied. Any $m$ by $n$ geometry with blocks that have zero diffusivity would require fewer unknown interface functions than an $m$ by $n$ geometry without blocks that have zero diffusivity. This would decrease the size of the linear system (53). For example, for the geometry presented in Figure 4(a), with $N_{x}=N_{y}=4$, the linear system has dimension $2304$ by $2304$. However, if the dark grey blocks had a diffusivity of zero instead of $0.1$, the linear system has dimension $1155$ by $1155$. Furthermore, the extended semi-analytical method presented in section 5.5 could also be applied to geometries which contain blocks with zero diffusivity to further reduce the dimension of the linear system (53). For example, for the formulation presented in Figure 4(b), if the dark grey blocks were set to have a diffusivity of zero, removing the interfaces between dark grey and light grey blocks results in a linear system of dimension 327 by 327.

6 Conclusion

In this paper we have presented a semi-analytical method for solving the classical boundary value problem arising from the homogenization of a two-dimensional, pixelated, locally-isotropic, heterogeneous, periodic domain. Our approach involves reformulating the heterogeneous boundary value problem on the unit cell as a family of boundary value problems on the homogeneous blocks by introducing unknown functions representing the flux at the interfaces between adjacent blocks. We then solve each homogeneous boundary value problem via standard techniques to yield an analytical solution depending on the integrals involving unknown interface functions. By applying an appropriate numerical quadrature rule to these integrals and enforcing continuity of the solution across each interface, the solution of the original heterogeneous boundary value problem can be computed, which allows the effective diffusivity to be calculated. Our numerical experiments demonstrated our semi-analytical method yields accurate results when applied to standard test problems from the literature. For all of these test problems considered in this work, we found that our semi-analytical method is faster and/or more accurate than a standard finite volume method. The primary reason for this is that the linear system for the semi-analytical method is much smaller than the linear system for the finite volume method as unknowns are located only along interfaces between adjacent blocks and not in the interior of the blocks (see Figure 1(d)–(e)). We also demonstrated that our semi-analytical method can be applied to complex heterogeneous geometries by pixellating and representing the geometry as a grid of blocks, as demonstrated in Table 4.

While all of the heterogeneous media considered in this paper consist of only two different materials, our semi-analytical method is not restricted to such media as the diffusivity in each block can be different to all other blocks. Furthermore, we only considered periodic boundary conditions (15)–(18), however, our semi-analytical method can be modified to accomodate additional forms of boundary conditions used in other homogenization techniques [7] and other steady-state diffusion problems on two-dimensional, block heterogeneous media. There is also room for improvement in our semi-analytical method that may further improve its efficiency and accuracy. For example: (i) allowing for blocks of different sizes where not all interfaces are aligned would reduce the number of interface functions and thus reduce the size of the linear system (53), as discussed in section 5.5; (ii) applying a quadrature rule to approximate the integrals involving the interface functions that allows for more than $2N_{x}-3$ terms to be taken in the summations would likely improve the accuracy; (iii) applying an iterative technique to the solution of the linear system, as demonstrated in [28], may allow for improvement in the efficiency of our semi-analytical method. In our future work, we plan to investigate the effect of the accuracy of an effective diffusivity tensor on the accuracy of the coarse-scale solution of diffusion problems in complex heterogeneous media and the use of our semi-analytical method for alleviating bottlenecks in coarse-grained models of non-periodic heterogeneous media, where the domain is decomposed into several sub-domains, which are individually homogenized.

Acknowledgments

The first author acknowledges discussions with Nadja Ray and Andreas Rupp from Friedrich-Alexander University of Erlangen-Nürnberg on randomly generating aggregated heterogeneous domains. All authors acknowledge the helpful comments of the anonymous reviewers and editors that helped improve the quality of the manuscript. The second and third authors acknowledge funding from the Australian Research Council (DE150101137, DP150103675).

Appendix A Boundary conditions for blocks in bottom row and right column

The boundary conditions (36)–(37) are valid for all blocks except those in either the bottom row ($i=m$) or right column $(j=n)$. For the blocks in the bottom row and right column, the boundary conditions are:

[TABLE]

for $j=1,\ldots,n-1$,

[TABLE]

for $i=1,\ldots,m-1$, and

[TABLE]

Appendix B Coefficients for blocks in bottom row and right column

The coefficients (41)–(44) are valid for all blocks except those in either the bottom row ($i=m$) or right column $(j=n)$. For the blocks in the bottom row and right column, the coefficients are defined as:

[TABLE]

for $j=1,\ldots,n-1$,

[TABLE]

for $i=1,\ldots,m-1$, and

[TABLE]

Appendix C Solvability conditions for blocks in bottom row and right column

The solvability conditions (45) are valid for all blocks except those in either the bottom row ($i=m$) or right column $(j=n)$. For the blocks in the bottom row and right column, the solvability conditions are:

[TABLE]

for $j=1,\ldots,n-1$,

[TABLE]

for $i=1,\ldots,m-1$, and

[TABLE]