Numerical investigation into coarse-scale models of
diffusion in complex heterogeneous media

According to the homogenization literature (Hornung, 1997; Szymkiewicz, 2005; Davit et al., 2013; Carr et al., 2013a; Carr and Turner, 2014; Carr et al., 2016, 2017), the effective diffusivity $\mathbf{D}_{\text{eff}}(x,y)$ used in the coarse-scale model (2) can be calculated using the solution of an appropriate boundary value problem (or homogenization problem). In this work, we consider three common formulations for the homogenization problem involving periodic, confined and uniform boundary conditions (see Figure 2) (Szymkiewicz, 2013; Renard and de Marsily, 1997). In this section, we discuss these three types of boundary conditions, their associated boundary value problems and how they are used to calculate the effective diffusivities for an arbitrary homogenization cell $\Omega_{C}=[x_{a},x_{b}]\times[y_{a},y_{b}]$ and in section 3.4 we discuss our numerical solution strategy for solving these boundary value problems.

We define a square mesh over the domain $\Omega=[x_{0},x_{m}]\times[y_{0},y_{m}]$ that consists of $N_{f}+1$ uniformly spaced nodes in both the $x$ and $y$ directions. The $x$-coordinates of the nodes are $x^{p}=x_{0}+ph_{f}$ for $p=0,\ldots,N_{f}$ and the $y$-coordinates are $y^{q}=y_{0}+qh_{f}$ for $q=0,\ldots,N_{f}$, where $h_{f}=(x_{m}-x_{0})/N_{f}=(y_{m}-y_{0})/N_{f}$. The spacing $h_{f}$ is chosen to ensure that nodes coincide with interfaces between blocks. We define the solution approximation at the node located at $(x^{p},y^{q})$ as $u^{p,q}$. We decompose the domain into square elements by connecting the four nodes located at $(x^{p-1},y^{q-1}),(x^{p},y^{q-1}),(x^{p-1},y^{q})$ and $(x^{p},y^{q})$ for $p=1,\ldots,N_{f}$ and $q=1,\ldots,N_{f}$. We define a set of edges $\mathcal{E}_{f}^{p,q}$ (see Figure 3) around the node located at $(x^{p},y^{q})$ and connect them to form a control volume. The faces are constructed by connecting the centroid of an element to the midpoint of its edges (see Figure 3).

We integrate the governing partial differential equation (1) over the control volume and apply the divergence theorem and midpoint rule to yield a finite volume equation. For a node located at $(x^{p},y^{q})$, except for those located on a boundary in which Dirichlet conditions are imposed, the finite volume equation takes the form:

where $\Phi_{s}^{p,q}$ is the flux through the edge $s$ and is calculated as:

where $(x_{s}^{p},y_{s}^{q})$ is the midpoint of edge $s$ and $\mathbf{n}^{p,q}_{s}$ is the outward facing normal vector to edge $s$ with length $h_{f}/2$. We note that the flux corresponding to the edges located on the boundary is zero as the Neumann boundary conditions that we consider in this paper are homogeneous. For any node located on a boundary with Dirichlet conditions, the solution approximations corresponding to the nodes located along the boundary can be calculated directly, so no finite volume equation is required for these nodes. Any node adjacent to one of these boundary nodes is modified to account for these known solution values.

On each element, we approximate the solution with a bilinear interpolant $\widetilde{u}(x,y)$. For an element $[x^{p-1},x^{p}]\times[y^{q-1},y^{q}]$, the bilinear interpolant takes the form:

The coefficients $a_{1},a_{2},a_{3}$ and $a_{4}$ can be computed by noting that within each element, the solution values at the corners are known and by enforcing that the solution value as computed by the bilinear interpolant is equal to these values on each of the four corners, the following linear system can be generated:

This allows for the coefficients $a_{1},a_{2},a_{3}$ and $a_{4}$ to be expressed in terms of the solution approximations $u^{p-1,q-1},u^{p,q-1},u^{p-1,q}$ and $u^{p,q}$ as:

where $\alpha_{1},\ldots,\alpha_{16}$ can be computed in terms of the coordinates of the four corners of the element. This bilinear interpolant allows $\nabla\widetilde{u}(x,y)$ to be computed as:

This gradient is then used to approximate the flux term appearing in the finite volume equation (17).

The finite volume scheme for the coarse scale model is similar to that of the fine-scale model except a coarser mesh can be used as $\mathbf{D_{\text{eff}}}(x,y)$ varies at a coarser scale compared to $D(x,y)$.

We define a square mesh for the coarse-scale model over the domain $\Omega=[x_{0},x_{m}]\times[y_{0},y_{m}]$ that consists of $N_{c}+1$ uniformly spaced nodes in both the $x$ and $y$ directions. The $x$-coordinates of the nodes are $x^{p}=x_{0}+ph_{c}$ for $p=0,\ldots,N_{c}$ and the $y$-coordinates are $y^{q}=y_{0}+qh_{c}$ for $q=0,\ldots,N_{c}$, where $h_{c}=(x_{m}-x_{0})/N_{c}=(y_{m}-y_{0})/N_{c}$. The spacing $h_{c}$ is chosen to ensure that each control volume edge is located entirely within one homogenization cell. We define the solution at the node located at $(x^{p},y^{q})$ as $U^{p,q}$. This solution is approximated on each element with a bilinear interpolant $\widetilde{U}(x,y)$. The remainder of the finite volume method for the coarse-scale model is implemented identically to the finite volume method for the fine-scale model.

The finite volume equations for each of the nodes for the fine-scale and coarse-scale models are collected to form a system of linear equations. For the fine-scale model, the linear system takes the form:

where $\mathbf{u}$ is a vector containing the discrete unknown values $u^{p,q}$ at each node, except for those corresponding to Dirichlet boundary conditions (as the value of the solution is known at these nodes), $\mathbf{A}_{f}$ is a matrix whose entries are determined using the finite volume equation (17) and the boundary conditions and $\mathbf{b}_{f}$ is a vector whose entries are determined by the boundary conditions. The solution of this linear system is computed using the backslash operator in MATLAB. Similarly, the linear system for the coarse-scale model takes the form:

where $\mathbf{U}$ is a vector containing the discrete unknown values $U^{p,q}$ at each node, except for those corresponding to Dirichlet boundary conditions and $\mathbf{A}_{c}$ and $\mathbf{b}_{c}$ are formed similarly to the fine-scale model. The solution of equation (27) is also calculated using the backslash operator in MATLAB.

We define a square mesh over the homogenization cell $\Omega_{C}=[x_{a},x_{b}]\times[y_{a},y_{b}]$ that consists of $N_{e}+1$ uniformly spaced nodes in both the $x$ and $y$ directions. The $x$-coordinates of the nodes are $x^{p}=x_{a}+ph_{e}$ for $p=0,\ldots,N_{e}$ and the $y$-coordinates are $y^{q}=y_{0}+qh_{e}$ for $q=0,\ldots,N_{e}$, where $h=(x_{b}-x_{a})/N_{e}=(y_{b}-y_{a})/N_{e}$. The spacing $h_{e}$ is chosen to ensure that nodes coincide with interfaces between blocks. The governing equation for the effective diffusivity, (3), (9) or (14), depending on the chosen boundary conditions, is solved using the finite volume scheme for the fine-scale model with slight modifications. For periodic conditions, the finite volume equation for a node located at $(x^{p},y^{q})$ takes the form:

where $(x_{s}^{p},y_{s}^{q})$ is the midpoint of edge $s$, $\mathcal{E}_{e}^{p,q}$ is the set of edges surrounding the node, $\mathbf{n}^{p,q}_{s}$ is the outward facing normal vector of edge $s$ with length $h_{e}/2$ and $\nabla\xi=[1,0]^{T}$ if $\xi=x$ and $\nabla\xi=[0,1]^{T}$ if $\xi=y$. The nodes along the boundaries are treated by modifying the finite volume equation (28) to account for the periodic boundary conditions, as detailed by Carr et al. (2013a). The zero mean condition (6) is accommodated by approximating the integrals with the midpoint rule, leading to the summation:

where $\psi_{p,q}^{(\xi)}=\psi^{(\xi)}(x^{p},y^{q})$. For uniform and confined conditions, the finite volume equation for a node located at $(x^{p},y^{q})$, except for those located on a boundary in which Dirichlet conditions are imposed, takes the form:

where:

We note that the flux corresponding to the edges located on the boundary is zero as the Neumann boundary conditions that occur in the confined conditions are homogeneous. For any node located on a boundary with Dirichlet conditions, the solution approximations corresponding to the nodes located along the boundary can be calculated directly, so no finite volume equation is required for these nodes. Any node adjacent to one of these boundary nodes is modified to account for these known solution values.

The finite volume equations for each of the nodes for the effective diffusivity are collected to form a system of linear equations of the form:

where $\boldsymbol{\psi}^{(\xi)}$ is a vector containing the discrete unknown values $\psi^{(\xi)}_{p,q}$ at each node, except for those corresponding to Dirichlet boundary conditions, $\mathbf{A}_{e}$ is a matrix whose entries are determined using the finite volume equations (28)–(30) and the boundary conditions and $\mathbf{b}_{e}$ is a vector whose entries are determined by the boundary conditions. The solution of this linear system is calculated using the backslash operator in MATLAB. The effective diffusivities are then calculated by approximating the integrals appearing in equations (7)–(8) and (12)–(13) with a midpoint rule, as was used previously by March et al. (2021), leading to the following approximations when periodic boundary conditions are imposed:

where $(x_{\text{cen}}^{p},y_{\text{cen}}^{q})$ is the centroid of the $(p,q)$th element, $\mathbf{e}_{1}=[1,0]^{T}$ and $\mathbf{e}_{2}=[0,1]^{T}$. The corresponding approximations when uniform and confined boundary conditions are imposed are:

In this test we consider the previously mentioned 100 different geometries with volume fractions $\varepsilon_{1}$ of $25\%,50\%$ and $75\%$ and diffusivity ratios $\varepsilon_{2}$ of $10$ and $100$. We consider the homogenization parameters $k=1,2,3,4,5,6,10,12,15,20,30$ and all three types of boundary conditions on the homogenization cell. The mesh size for the coarse-scale model is set to be $N_{c}=N_{f}$. Preliminary results did not show significant differences in errors across different volume fractions $\varepsilon_{1}$, so only the results for $\varepsilon_{1}=50\%$ are shown in Figure 4. For the diffusivity ratio $\varepsilon_{2}=10$ and homogenization parameters $k\leq 20$, the periodic and confined boundary conditions performed similarly and were more accurate than the uniform boundary conditions. For the homogenization parameters $k>20$, all three boundary conditions yielded similar results. In the case of diffusivity ratio $\varepsilon_{2}=100$ and homogenization parameters $k\leq 15$, the periodic and confined boundary conditions performed approximately equally and were significantly more accurate than the uniform boundary conditions. For $15\leq k\leq 30$, the confined boundary conditions performed slightly better than the periodic boundary conditions, which performed better than the uniform boundary conditions. As both plots in Figure 4 look similar, the different coarse-scale boundary conditions (39)–(40) appear to have little effect on the performance of the coarse-scale model.

We also compare the average flux calculated using the coarse-scale model to the average flux calculated using the fine-scale model. To calculate the average flux for the fine-scale model, we first calculate the gradient $\nabla\widetilde{u}(x,y)$ at the centroid of the element $[x^{p-1},x^{p}]\times[y^{q-1},y^{q}]$ using the finite volume scheme discussed in section 3 by evaluating equation (25) at the centroid of the element $(x_{\text{cen}}^{p},y_{\text{cen}}^{q})$ for $p=1,\ldots,N_{f}$ and $q=1,\ldots,N_{f}$. We then calculate the flux vector $\mathbf{J}_{f}^{p,q}$ as:

for $p=1,\ldots,N_{f}$ and $q=1,\ldots,N_{f}$. The flux vector for the coarse-scale model, $\mathbf{J}_{c}^{p,q}$, is calculated similarly:

for $p=1,\ldots,N_{c}$ and $q=1,\ldots,N_{c}$. We note that we use the fine-scale diffusivity $D(x,y)$ in the calculation of the flux corresponding to the coarse-scale model, rather than the effective diffusivity $\mathbf{D}_{\text{eff}}(x,y)$. In preliminary testing, we used $\mathbf{D}_{\text{eff}}(x,y)$ in equation (42), rather than the fine-scale diffusivity $D(x,y)$, however we found that the estimates of the flux vector were more accurate (as compared to the flux calculated using equation (41)) when $D(x,y)$ is used. We note that using $D(x,y)$ in the calculation of the flux for the coarse-scale model is reasonable, as $D(x,y)$ is known everywhere in the domain $\Omega$ and is used to calculate the effective diffusivity $\mathbf{D}_{\text{eff}}(x,y)$. We define the matrices $\mathbf{J}_{f}\in\mathbb{R}^{2\times N_{f}^{2}}$ and $\mathbf{J}_{c}\in\mathbb{R}^{2\times N_{c}^{2}}$, which contain the components of the flux vectors $\mathbf{J}_{f}^{p,q}$ and $\mathbf{J}_{c}^{p,q}$, respectively, on each element. To calculate the average flux vectors, we compute the mean of the $N_{f}^{2}$ and $N_{c}^{2}$ columns of $\mathbf{J}_{f}$ and $\mathbf{J}_{c}$, denoted as $\mathbf{\overline{J}}_{f}\in\mathbb{R}^{2\times 1}$ and $\mathbf{\overline{J}}_{c}\in\mathbb{R}^{2\times 1}$, respectively. The second error metric that we use is the flux-based error, defined as the error in the average flux between the coarse-scale and fine-scale solutions:

where $\varoslash$ represents Hadamard/element-wise division and $\vline\,\cdot\,\vline$ indicates absolute value. We then average this error over the $100$ different geometries:

We note that $\mathbf{E}^{J}\in\mathbb{R}^{2\times 1}$ and that the first component of $\mathbf{E}^{J}$ represents the average error of the flux in the $x$-direction and that the second component represents the average error of the flux in the $y$-direction. However, we consider only the first component $\mathbf{E}^{J}_{[1]}$ in the analysis for the BC1 (39) boundary conditions and second component $\mathbf{E}^{J}_{[2]}$ in the analysis for the BC2 (40) boundary conditions, as the flux in the $y$-direction for BC1 boundary conditions and the flux in the $x$-direction for BC2 boundary conditions is very small, which can cause large relative errors in the components of $\mathbf{E}^{J}$ calculated using these fluxes.

As seen in Figure 5, for both BC1 and BC2 conditions, all three boundary conditions performed comparably in terms of accuracy for the diffusivity ratio $\varepsilon_{2}=10$ and that the periodic and confined conditions performed similarly for $\varepsilon_{2}=100$, with uniform conditions performing slightly better for some values of $k$. Similarly to Figure 4, the flux-based error increases as $k$ increases, from approximately $0.2$ to $0.75$ for $\varepsilon_{2}=10$ and $0.3$ to $1$ for $\varepsilon_{2}=100$. We also note that the flux-based error shown in Figure 5 is much higher than the solution-based error shown in Figure 4. We believe that the larger flux-based error is due to the the tendency of coarse-scale solutions to “smoothen” over fine-scale details (see Figure 7), which can cause the flux calculated using the coarse-scale solution to differ significantly from the flux calculated using the fine-scale solution.

In this test, we consider the same 100 different geometries from Test 1. We consider only the homogenization parameters $k=1,2,3,4,5,6$ and all three types of boundary conditions on the homogenization cell. The mesh size for the coarse-scale model is set to be $N_{c}=r$, which is the coarsest mesh that allows for each element to be contained entirely within homogeneous material. This differs from Test 1, as the simulations in Test 2 use a coarser mesh. As the mesh size varies with $k$, to calculate the solution-based error (37), bilinear interpolation is used to interpolate the coarse-scale solution onto the fine-scale mesh. The choices of $k$ are considered only up to $k=6$, as opposed to the $k=30$ used in Test 1 as the larger $k$ values correspond to very coarse meshes. To compare the results for Test 2 against those for Test 1, we calculate the error ratio $R=E_{\text{var}}/E_{\text{con}}$, where $E_{\text{var}}$ is the mean error (44) as calculated using the varying mesh and $E_{\text{con}}$ is the mean error (44) as calculated using the constant mesh.

Similarly to Test 1, preliminary results did not show significant differences in errors across different volume fractions $\varepsilon_{1}$, so only the results for $\varepsilon_{1}=50\%$ are shown in Figure 6. The results varied significantly for different combinations of coarse-scale boundary conditions. The purpose of this test is to determine what effect the use of a coarser mesh has on the relative error. We expect that the error would increase with a coarser mesh, as the spatial discretisation error associated with finite volume methods would generally increase as the number of nodes decreases. For both coarse-scale boundary conditions, the relative errors for diffusivity ratios $\varepsilon_{2}=10$ were generally lower than the corresponding relative errors for diffusivity ratios $\varepsilon_{2}=100$.

For BC1 (Figure 6(a)), the error ratio is between approximately $1.02$ and $1.08$ for any diffusivity ratio or volume fraction and for BC2 (Figures 6(b)), the error ratio is between approximately $1.02$ and $1.1$. There is no clear relationship between the homogenization parameter $k$ and the error ratio for any choices of coarse-scale boundary conditions, boundary conditions on the homogenization cell, volume fraction $\varepsilon_{1}$ or diffusivity ratio $\varepsilon_{2}$. However as the error ratio is between $1.02$ and $1.1$, the spatial discretisation error introduced through the use of a coarser mesh is between $2\%$ and $10\%$ of the homogenization error and thus the homogenization error is approximately $10$ to $50$ times larger than the spatial discretisation error associated with our implemented finite volume method.

Based on the results generated and presented in section 4 and other properties of the different boundary conditions for the homogenization cell, we make some recommendations for homogenizing block heterogeneous media. We assume that for a given problem, the coarse-scale boundary conditions, initial condition, $m$ and $D_{i,j}$ for all $i=1,\ldots,m$ and $j=1,\ldots,m$ are all known. We also assume that there exists a value $N_{\text{max}}$, such that a mesh comprised of $N_{\text{max}}+1$ uniformly spaced nodes in both the $x$ and $y$ directions is the most refined mesh that is computationally feasible for a given problem. We also assume that any mesh comprised of $N_{\text{fea}}+1$ uniformly spaced nodes in both the $x$ and $y$ directions, where $N_{\text{fea}}<N_{\text{max}}$ is computationally feasible. Given these assumptions, we can make the following recommendations:

• •

  For $N_{\text{max}}>m$, the fine-scale solution model should be used with $N_{f}=N_{\text{max}}$;

• •

  For $N_{\text{max}}\leq m$, the coarse-scale model should be used with homogenization parameter $k=\lceil{m/N_{\text{max}}}\rceil$, periodic boundary conditions on the homogenization cells; and a mesh constructed by setting $N_{c}=m/k$.

We recommend that the fine-scale model be used if $N_{\text{max}}>m$ as the fine-scale solution does not require any homogenization. We recommend the use of periodic conditions in the case of $N_{\text{max}}\leq m$ for the following reasons:

• •

  The periodic and confined boundary conditions have comparable solution-based errors that are significantly lower than the error for the uniform conditions across all tests, as seen in Figure 4.

• •

  The periodic and confined boundary conditions compute the correct effective diffusivities for layered media, whereas the uniform boundary conditions give an incorrect value for the effective diffusivity in the direction perpendicular to the layers (Szymkiewicz, 2013).

• •

  The implementation of the finite volume methods for the calculation of the effective diffusivities corresponding to homogenization cells with either periodic or uniform conditions requires the inversion of only one coefficient matrix, as the coefficient matrix $\mathbf{A}_{e}$ appearing in the linear system (32) is identical for both $\xi=x$ and $\xi=y$, whereas the linear system that arises from confined conditions uses a different coefficient matrix for both $\xi=x$ and $\xi=y$ and thus requires two matrix inversions.

• •

  The periodic and uniform conditions are guaranteed to return symmetric positive definite effective diffusivities, which ensures that the principal axes of diffusion and the diffusivities in those directions are physical, whereas the confined conditions do not (Szymkiewicz, 2013).

We recommend $k=\lceil{m/N_{\text{max}}}\rceil$ as smaller values of $k$ are more accurate and it is the smallest value of $k$ giving a computationally-feasible mesh.

We consider the solution of a steady-state diffusion problem on binary media with volume fraction $\varepsilon_{1}=50\%$, diffusivity ratio $\varepsilon_{2}=100$, $m=60$, $N_{c}=m$ and $k=1,3,5,10,15,30$. Figure 7 shows these solutions as well as the original geometry with each principal axes of diffusion overlaid on each homogenization cell. For low levels of homogenization, $k=3$ and $k=5$, the solutions compare well to the fine-scale solution $k=1$, although the contours representing $U=0.1,0.2$ and $0.3$ are significantly different in shape. For $k=10$, the majority of the solution field is signficantly different than that of the fine-scale solution, with the possible exception of the contours representing $U=0.7,0.8$ and $0.9$. For $k=15$, the solution field bears only a vague resemblance to the fine-scale solution as the homogenized domain is comprised of only $16$ homogenization cells, indicating that a smaller value of $k$ is required to adequately capture the fine-scale behaviour. For $k=30$ the solution is almost linear as there are only four different effective diffusivities across the entire domain.

We repeat the experiments presented in section 5.1 for $k=2,3,6$ with a mesh constructed using $N_{c}=r$ as discussed in section 4.2. We consider only these homogenization parameters as we found that for larger values of $k$ that the corresponding mesh is too coarse to generate an accurate solution. We present these solutions in Figure 8 and note that the appropriate comparison in this figure is between pairs of figures with the same values of $k$, as each solution is computed with the same homogenization parameter and the only difference is in the mesh used in the finite volume method to solve the coarse-scale model. For $k=2$, the solutions are very similar, with the only noticeable difference in the solutions occurring in the contour representing $U=0.4$ and $U=0.3$. Similar results hold for $k=3$, with the only noticeable difference in the contours of the solution representing $U=0.4$ and $U=0.3$. However for $k=6$, the solutions differ significantly, indicating that the mesh with $N_{c}=10$ that is used in this simulation is too coarse to generate an accurate solution.

In this section, motivated by fracturing in shale gas reservoirs (Zhou et al., 2017; Kong et al., 2018), chalk engineering (Welch et al., 2015; Wattier et al., 2018) and cell wall collapse in wood drying (Carr et al., 2013a), we compare coarse-scale models of diffusion on two different geometries, in which the second geometry is a slight modification of the first geometry. The purpose of this section is to investigate what effects these modifications have on the solution of steady-state diffusion problems and to what extent these changes can be seen in the coarse-scale solution. We consider a simulation in the media presented in Figure 9. In this geometry there are three sections that have been changed to generate the modified geometry (highlighted using red ellipses). In two of these changes a number of low diffusivity dark grey blocks are replaced with high diffusivity light grey blocks and in the other change one light grey block has been replaced with a dark grey block. These first two changes allow for two continuous paths of light grey material joining the left and right boundaries of the domain while the other change has yielded a section of light grey material that is entirely surrounded by dark grey material.

In the comparison of the fine-scale solutions represented in Figure 9(d) and Figure 9(j), the modifications to the geometry have caused the size of the regions between the contours corresponding to $u=1$ and $u=0.9$ and between $u=0.1$ and $u=0$ to decrease and the size of the other regions to increase. For the solutions corresponding to $k=3$ (Figures 9 (e) and (f)), they are both similar to the respective fine-scale solutions, with the largest difference occurring in the regions between the contours $U=0.3$ and $U=0.2$. For $k=5$, there are significant differences in the region between the contours $U=1$ and $U=0.9$ for the coarse-scale solution Figure 9(l) compared to the fine-scale solution Figure 9(j), whereas differences corresponding to the comparison of the coarse-scale solution Figure 9(f) to the fine-scale solution Figure 9(d) are less pronounced. However, the reverse situation is true for the region between the contours $U=0$ and $U=0.1$. Overall, the process of homogenization has yielded similar qualitative differences in the coarse-scale solutions for both the original geometry Figure 9(a) and the modified geometry Figure 9(g).