3.1. Two spatial scales

Applying the MMS as in Chapman & McBurnie (Reference Chapman and McBurnie2011) and Richardson & Chapman (Reference Richardson and Chapman2011), we consider the spatial domain on two distinct length scales: the macroscale $\boldsymbol {x}:= \hat {\boldsymbol {x}}$, relative to which the porous medium is of unit length, and a microscale coordinate $\boldsymbol {y}$, in which we impose strict periodicity. The latter can be achieved via a mapping that transforms each cell (comprising the porous material) to a tessellating periodic cell. Here, we choose to transform to a square cell of unit area (see figure 2a). As such, our mapping will stretch the obstacles comprising the macroscale filter by a factor of $1/(\epsilon a(x_1))$ in the longitudinal direction and by $1/\epsilon$ in the transverse direction, i.e.

(3.3a,b)\begin{equation} \frac{\mathrm{d} y_1}{\mathrm{d}\kern0.06em x_1} = \frac{1}{\epsilon a(x_1)} \quad \text{and} \quad \frac{\mathrm{d} y_2}{\mathrm{d}\kern0.06em x_2} = \frac{1}{\epsilon}, \end{equation}

thus motivating the definitions

(3.4a,b)\begin{equation} {y_1}:= \frac{1}{\epsilon} \int^{x_1} \! \frac{\mathrm{d}s}{a(s)} \quad \text{and} \quad y_2 := \frac{x_2}{\epsilon}. \end{equation}

Figure 2.

An arbitrary cell within the porous medium (orange rectangle in figure 1) represented in (a) transformed microscale coordinates and (b) naive microscale coordinates. The transformed microscale coordinates $y_1$ and $y_2$ and the naive microscale coordinates $Y_1$ and $Y_2$ are related via (3.4a,b) and (3.5). The transformed microscale coordinates allow the slow variation in cell width $a$ to be scaled out of the cell problem, such that each naive rectangular cell is transformed into a square.

Note that any arbitrary distribution of obstacles in the longitudinal directions (i.e. arbitrary longitudinal heterogeneity) can be imposed via the functions $a(x_1)$ and $\lambda (x_1)$, while in the transverse direction the porous medium is exactly periodic. Additional heterogeneity in the transverse direction can be considered through a more general mapping (Richardson & Chapman Reference Richardson and Chapman2011).

To understand the implications of the mapping (3.4a,b) on a single cell, we note that (on the macroscale) the domain of the porous medium, $\boldsymbol {x}\in \varOmega$, comprises a fluid domain $\varOmega _f$ and a complementary solid domain $\varOmega _s$. A single cell can be obtained by discretising the porous medium into rectangular cells of height $\epsilon$ and width $\epsilon a(x_1)$. For any single rectangular cell in the domain, the transformation (3.4a,b) yields a transformed microscale, $\boldsymbol {y}$. On the transformed microscale, each cell $\omega$ comprises a fluid phase $\omega _f(x_1)$ and solid obstacles, the union of which is denoted $\omega _s(x_1):= \omega \setminus \omega _f(x_1)$. The fluid–solid interface $\partial \omega _s(x_1)$ is the union of the boundary of the obstacles. Each cell has four additional boundaries that separate it from neighbouring cells. We denote the top and bottom boundaries $\partial \omega _{=}$ and the left and right boundaries $\partial \omega _{\|}$ with the union of these being the unit cell boundary $\partial \omega$.

Following Chapman & McBurnie (Reference Chapman and McBurnie2011) and Richardson & Chapman (Reference Richardson and Chapman2011), we choose the microscale mapping such that the cell size is the same throughout the domain. Hence, for example, circles will approximately map to ellipses. Since the untransformed cell size varies spatially through the domain, this microscale mapping will lead to obstacles that vary by an $O(\epsilon )$ amount between neighbouring cells, but by an $O(1)$ amount over the macroscale. We systematically account for these variations using the methodology presented in Bruna & Chapman (Reference Bruna and Chapman2015) and Dalwadi et al. (Reference Dalwadi, Griffiths and Bruna2015).

Finally, we note that the transformed microscale variable can be difficult to interpret physically. As such, it will be helpful to define a ‘naive’ microscale coordinate

(3.5)\begin{equation} \boldsymbol{Y}:=\boldsymbol{x}/\epsilon, \end{equation}

in which each cell is of unit transverse height but of longitudinal width $a(x_1)$ (figure 2b). After completing the homogenisation procedure in the transformed microscale (3.4a,b), we will transform the relevant cell problems to the naive microscale (3.5), in order to present them more intuitively and subsequently solve them numerically. Note that the domains and boundaries in the (naive) rectangular $\boldsymbol {Y}$-cell will be denoted as in the square $\boldsymbol {y}$-cell, but with the addition of a superscript $\star$. Further, we emphasise that the microscale (cell) problems we derive and solve are not physical flow or transport problems, but rather mathematical constructs that enable us to invoke the MMS.

We now perform the homogenisation. Following the MMS, we take $\boldsymbol {x}$ and $\boldsymbol {y}$ to be independent spatial parameters. We therefore re-write all functions of $\hat {\boldsymbol {x}}$ as functions of $\boldsymbol {x}$ and $\boldsymbol {y}$: $\hat {\boldsymbol {v}}(\hat {\boldsymbol {x}}):=\boldsymbol {v}(\boldsymbol {x},\boldsymbol {y})$, $\hat {p}(\hat {\boldsymbol {x}}) := p(\boldsymbol {x},\boldsymbol {y})$, and $\hat {c}(\hat {\boldsymbol {x}},t) := c(\boldsymbol {x},\boldsymbol {y},t)$. We denote functions of $\boldsymbol {Y}$ (rather than of $\boldsymbol {y}$) with a superscript $\star$. Spatial derivatives then become

(3.6a)\begin{equation} \frac{\partial}{\partial \hat{x}_i} = \frac{\partial}{\partial x_i} + \frac{\sigma_{ij}}{\epsilon} \frac{\partial}{\partial y_j}, \end{equation}

for $i,j = 1, 2$, and where $\sigma _{ij}=(\boldsymbol {\sigma })_{ij}$ and

(3.6b)\begin{equation} \boldsymbol{\sigma}=\begin{pmatrix} \displaystyle\frac{1}{a(x_1)} & 0\\ 0 & 1 \end{pmatrix}. \end{equation}

Alternatively, in vector form, the spatial derivatives become

(3.6c)\begin{equation} \hat{\boldsymbol{\nabla}} := \boldsymbol{\nabla}_x+\frac{1}{\epsilon}\boldsymbol{\nabla}_y^a, \end{equation}

where $\boldsymbol {\nabla }_x$ is the gradient operator with respect to the coordinate $\boldsymbol {x}$ and where

(3.6d)\begin{equation} \boldsymbol{\nabla}_y^a:= \left(\frac{1}{a}\frac{\partial}{\partial y_1},\frac{\partial}{\partial y_2}\right)^\intercal, \end{equation}

is the gradient operator associated with the $\boldsymbol {y}$-coordinate transform. For a given quantity $Z(\boldsymbol {x},\boldsymbol {y},t) = Z^\star (\boldsymbol {x},\boldsymbol {Y},t)$, there are two different averages of interest: the intrinsic (fluid) average

(3.7)\begin{equation} \langle Z\rangle(\boldsymbol{x},t) := \frac{1}{|\omega_f(x_1)|} \! \int_{\omega_f(x_1)} Z(\boldsymbol{x},\boldsymbol{y},t) \, \mathrm{d}S_y \equiv \frac{1}{|\omega_f^{{\star}} (x_1)|}\! \int_{\omega_f^\star (x_1)} Z^\star(\boldsymbol{x},\boldsymbol{Y},t) \, \mathrm{d}S_Y, \end{equation}

where the total fluid area in the transformed cell $|\omega _f|$ (or naive cell $|\omega _f^{\star }|$) is a function of $a(x_1)$ and $\lambda (x_1)$, and the volumetric average

(3.8)\begin{equation} \frac{1}{|\omega(x_1)|}\int_{\omega(x_1)} \! Z(\boldsymbol{x},\boldsymbol{y},t)\, \mathrm{d} S_y \equiv \frac{1}{|\omega^{{\star}}(x_1)|}\int_{\omega^{{\star}}(x_1)} \! Z^\star(\boldsymbol{x},\boldsymbol{Y},t) \, \mathrm{d}S_Y, \end{equation}

where $|\omega |=1$ and $|\omega ^{\star }|=a$. Here, $\mathrm {d}S_y:=\mathrm {d} y_1\,\mathrm {d} y_2$ is an area element of the transformed microscale fluid region, $\mathrm {d}S_Y:=\mathrm {d}Y_1\,\mathrm {d}Y_2$ is an area element of the naive microscale fluid region and the porosity $\phi$ is

(3.9)\begin{equation} \phi(x_1) = \frac{|\omega_f(x_1)|}{|\omega(x_1)|}\equiv |\omega_f(x_1)| \left(= \frac{|\omega_f^{{\star}}(x_1)|}{|\omega^{{\star}}(x_1)|} \right). \end{equation}

Thus, $\langle c\rangle$ is the amount of solute per unit fluid area within the porous medium, while $\phi \langle c\rangle$, the volumetric average of the concentration, is the amount of solute per unit total area.

We define the average velocity, pressure and concentration as

(3.10a–c)\begin{align} \boldsymbol{V}(\hat{\boldsymbol{x}}) \equiv \boldsymbol{V}(\boldsymbol{x}) := \langle \boldsymbol{v} \rangle , \quad P(\hat{\boldsymbol{x}}) \equiv P(\boldsymbol{x}) := \langle p \rangle, \quad \text{and} \quad C(\hat{\boldsymbol{x}},t) \equiv C(\boldsymbol{x},t) := \langle c \rangle, \end{align}

respectively. Note that $\phi \boldsymbol {V}$ is the standard Darcy flux.

3.2. Flow problem

For a passive tracer, the flow problem (2.6a) does not depend on $c$. Using (3.6), (2.6) in an arbitrary cell become

(3.11a)$$\begin{gather} -\left(\boldsymbol{\nabla}_x+\frac{1}{\epsilon}\boldsymbol{\nabla}_y^a\right){p}+(\epsilon^2\nabla^2_x+\epsilon\boldsymbol{\nabla}_x\boldsymbol{\cdot}\boldsymbol{\nabla}_y^a +\epsilon\boldsymbol{\nabla}_y^a\boldsymbol{\cdot}\boldsymbol{\nabla}_x +(\nabla_y^a)^2)\boldsymbol{v} = \boldsymbol{0}, \quad {\boldsymbol{y}}\in\omega_f(x_1), \end{gather}$$

(3.11b)$$\begin{gather}(\epsilon\boldsymbol{\nabla}_x+\boldsymbol{\nabla}_y^a)\boldsymbol{\cdot}{{\boldsymbol{v}}} = 0, \quad {\boldsymbol{y}}\in\omega_f(x_1), \end{gather}$$

(3.11c)$$\begin{gather}\boldsymbol{v} = \boldsymbol{0}, \quad {\boldsymbol{y}}\in \partial\omega_s(x_1), \end{gather}$$

where (3.11b) has been multiplied by epsilon.

To proceed using the MMS, we must also impose periodicity of $\boldsymbol {v}$, $p$ and $c$ over a single microscale cell (i.e. local periodicity). Enforcing periodicity of all quantities at both the top and bottom, $\partial \omega _{=}$, and left and right, $\partial \omega _{\|}$, cell boundaries leads to

(3.11d)\begin{equation} \boldsymbol{v} \quad \text{and}\quad p \quad \text{periodic on} \quad \boldsymbol{y}\in\partial\omega_{=}\quad \text{and}\quad \partial\omega_{\|}. \end{equation}

We now seek an asymptotic solution to (3.11) by expanding $\boldsymbol {v}$ and $p$ in powers of $\epsilon$

(3.12a)$$\begin{gather} \boldsymbol{v}(\boldsymbol{x},\boldsymbol{y}) = \boldsymbol{v}^{(0)}(\boldsymbol{x},\boldsymbol{y})+\epsilon\boldsymbol{v}^{(1)}(\boldsymbol{x},\boldsymbol{y})+\cdots \quad \text{as}\ \epsilon\to0, \end{gather}$$

(3.12b)$$\begin{gather}p(\boldsymbol{x},\boldsymbol{y} ) = p^{(0)}(\boldsymbol{x},\boldsymbol{y})+\epsilon p^{(1)}(\boldsymbol{x},\boldsymbol{y})+\cdots \quad \text{as}\ \epsilon\to0. \end{gather}$$

Considering terms of $O(1/\epsilon )$ in (3.11a) gives

(3.13)\begin{equation} \boldsymbol{\nabla}_y^ap^{(0)}= \boldsymbol{0}, \end{equation}

from which we conclude the standard result that, at leading order, the pressure is uniform on the microscale: $p^{(0)} = p^{(0)}(\boldsymbol {x})$.

Considering terms of $O(1)$ in (3.11) gives

(3.14a)$$\begin{gather} -\boldsymbol{\nabla}_x p^{(0)}-\boldsymbol{\nabla}_y^a{p^{(1)}}+(\nabla_y^a)^2\boldsymbol{v}^{(0)} = \boldsymbol{0}, \quad {\boldsymbol{y}}\in\omega_f(x_1), \end{gather}$$

(3.14b)$$\begin{gather}\boldsymbol{\nabla}_y^a\boldsymbol{\cdot}{{\boldsymbol{v}^{(0)}}} = 0, \quad {\boldsymbol{y}}\in\omega_f(x_1), \end{gather}$$

(3.14c)$$\begin{gather}{\boldsymbol{v}^{(0)}} = \boldsymbol{0}, \quad {\boldsymbol{y}}\in \partial\omega_s(x_1), \end{gather}$$

with

(3.14d)\begin{equation} \boldsymbol{v}^{(0)}\quad \text{and}\quad p^{(1)} \quad \text{periodic on} \quad \boldsymbol{y}\in\partial\omega_{=}\quad \text{and}\quad \partial\omega_{\|}. \end{equation}

The form of (3.14) suggests that we can scale $\boldsymbol {\nabla }_x p^{(0)}$ out of the problem via the substitutions

(3.15a)$$\begin{gather} \boldsymbol{v}^{(0)} ={-}\boldsymbol{{\mathcal{K}}}(\boldsymbol{x}, \boldsymbol{y}) \boldsymbol{\cdot} \boldsymbol{\nabla}_x p^{(0)}, \end{gather}$$

(3.15b)$$\begin{gather}p^{(1)} ={-}\boldsymbol{\varPi}(\boldsymbol{x},\boldsymbol{y})\boldsymbol{\cdot}\boldsymbol{\nabla}_x p^{(0)} +\breve{p}(\boldsymbol{x}), \end{gather}$$

where $\breve {p}(\boldsymbol {x})$ is a scalar function, $\boldsymbol {{\mathcal {K}}}(\boldsymbol {x}, \boldsymbol {y})$ is a tensor function, and $\boldsymbol {\varPi }(\boldsymbol {x},\boldsymbol {y})$ is a vector function. Using (3.15a,b) and the fact that $p^{(0)}$ is independent of $\boldsymbol {y}$, (3.14a–c) becomes

(3.16a)$$\begin{gather} (\boldsymbol{\mathsf{I}}-\boldsymbol{\nabla}_y^a\otimes\boldsymbol{\varPi}+(\nabla^a_y)^2\boldsymbol{{\mathcal{K}}})\boldsymbol{\cdot}\boldsymbol{\nabla}_x p^{(0)} = \boldsymbol{0}, \quad {\boldsymbol{y}}\in\omega_f(x_1), \end{gather}$$

(3.16b)$$\begin{gather}(\boldsymbol{\nabla}_y^a\boldsymbol{\cdot}{\boldsymbol{{\mathcal{K}}}})\boldsymbol{\cdot}\boldsymbol{\nabla}_x p^{(0)} = 0, \quad {\boldsymbol{y}}\in\omega_f(x_1), \end{gather}$$

(3.16c)$$\begin{gather}\boldsymbol{{\mathcal{K}}}\boldsymbol{\cdot}\boldsymbol{\nabla}_x p^{(0)} = \boldsymbol{0}, \quad {\boldsymbol{y}}\in \partial\omega_s(x_1), \end{gather}$$

with

(3.16d)\begin{equation} \mathcal{K}_{ij}:=\left(\boldsymbol{\mathcal{K}}\right)_{ij} \quad \text{and} \quad \varPi_{i}:=\left(\boldsymbol{\varPi}\right)_{i} \quad \text{periodic on} \quad \boldsymbol{y}\in\partial\omega_{=} \quad \text{and}\quad \partial\omega_{\|}, \end{equation}

where $\boldsymbol{\mathsf{I}}$ is the identity tensor and where

(3.16e)\begin{equation} \left(\boldsymbol{\nabla}_y^a\otimes\boldsymbol{\varPi}\right)_{ij} = \sigma_{ik} \frac{\partial \varPi_j}{\partial y_k} \quad \text{and} \quad(\boldsymbol{\nabla}_y^a\boldsymbol{\cdot}{\boldsymbol{{\mathcal{K}}}})_{i} = \sigma_{jk} \frac{\partial {\mathcal{K}}_{ji}}{\partial y_k}. \end{equation}

Note that, in the above, we have adopted the summation convention; we will adopt the summation convention throughout this manuscript. Equations (3.16) must hold for arbitrary $\boldsymbol {\nabla }_x p^{(0)}$, hence $\boldsymbol {{\mathcal {K}}}(\boldsymbol {x}, \boldsymbol {y})$ and $\boldsymbol {\varPi }(\boldsymbol {x},\boldsymbol {y})$ must satisfy the system

(3.17a)$$\begin{gather} \boldsymbol{\mathsf{I}}-\boldsymbol{\nabla}_y^a\otimes\boldsymbol{\varPi}+(\nabla^a_y)^2\boldsymbol{{\mathcal{K}}} = \boldsymbol{0}, \quad {\boldsymbol{y}}\in\omega_f(x_1), \end{gather}$$

(3.17b)$$\begin{gather}\boldsymbol{\nabla}_y^a\boldsymbol{\cdot}{\boldsymbol{{\mathcal{K}}}} = \boldsymbol{0}, \quad {\boldsymbol{y}}\in\omega_f(x_1), \end{gather}$$

(3.17c)$$\begin{gather}\boldsymbol{{\mathcal{K}}} = \boldsymbol{0}, \quad {\boldsymbol{y}}\in \partial\omega_s(x_1), \end{gather}$$

with

(3.17d)\begin{equation} {\mathcal{K}}_{ij} \quad \text{and} \quad \varPi_{i} \quad \text{periodic on} \ \boldsymbol{y}\in\partial\omega_{=} \quad \text{and}\quad \partial\omega_{\|}. \end{equation}

In general, (3.17) must be solved numerically for each desired cell geometry (i.e. pairs of $a$ and $\lambda$). Note that (3.17) are independent of $\boldsymbol {\nabla }_xp^{(0)}$, justifying our scalings in (3.15).

To derive a macroscale relationship between velocity and pressure from (3.15a), we expand the averaged quantities defined in (3.10a–c) in powers of $\epsilon$

(3.18a)$$\begin{gather} \boldsymbol{V}(\hat{\boldsymbol{x}}) = \boldsymbol{V}^{(0)}(\hat{\boldsymbol{x}})+\epsilon\boldsymbol{V}^{(1)}(\hat{\boldsymbol{x}})+\cdots\quad \text{as}\ \epsilon\to0, \end{gather}$$

(3.18b)$$\begin{gather}P(\hat{\boldsymbol{x}}) = P^{(0)}(\hat{\boldsymbol{x}})+\epsilon P^{(1)}(\hat{\boldsymbol{x}})+\cdots\quad\ \text{as}\ \epsilon\to0. \end{gather}$$

Note that

(3.19)\begin{equation} P^{(0)}(\hat{\boldsymbol{x}}) = \langle p^{(0)}(\boldsymbol{x})\rangle \equiv p^{(0)}(\boldsymbol{x}), \end{equation}

since $p^{(0)}$ is independent of $\boldsymbol {y}$. We then take the intrinsic average of (3.15a) to determine that the leading-order macroscale velocity depends on gradients in the leading-order macroscale pressure according to Darcy's law

(3.20a)\begin{equation} \phi\boldsymbol{V}^{(0)} ={-}\boldsymbol{\mathsf{K}}(\phi, a)\boldsymbol{\cdot}\hat{\boldsymbol{\nabla}}P^{(0)}, \end{equation}

where we have introduced the macroscale permeability tensor

(3.20b)\begin{equation} \boldsymbol{\mathsf{K}}(\phi, a) := \phi\langle{\boldsymbol{\mathcal{K}}}\rangle, \end{equation}

and where $\phi$ and $a$ are known functions of $\hat {{x}}_1$. Prescribing both $\phi$ and $a$ determines $\lambda$ via a simple geometric relation, specific to the chosen geometry of the porous material.

Being averaged in $\boldsymbol {y}$, (3.20a) depends on $\hat {\boldsymbol {x}} = \boldsymbol {x}$ only and we have therefore replaced $\boldsymbol {\nabla }_x$ with $\hat {\boldsymbol {\nabla }}$. If the cell geometry has symmetric reflectional symmetry along both the $y_1$ and $y_2$ axes, the symmetry of the boundary conditions implies that $\boldsymbol{\mathsf{K}}$ is diagonal. If it is additionally the case that $a=1$, then $\boldsymbol{\mathsf{K}}$ reduces to a scalar multiple of $\boldsymbol{\mathsf{I}}$.

Equation (3.20a) provides two equations for three unknowns. To develop another constraint in terms of $\boldsymbol {V}^{(0)}$ and $P^{(0)}$, consider the $O(\epsilon )$ terms from (3.11b) and (3.11c)

(3.21a)$$\begin{gather} \boldsymbol{\nabla}_x\boldsymbol{\cdot}\boldsymbol{v}^{(0)}={-}\boldsymbol{\nabla}_y^a\boldsymbol{\cdot}\boldsymbol{v}^{(1)},\quad \boldsymbol{y}\in\omega_f \end{gather}$$

(3.21b)$$\begin{gather}\boldsymbol{v}^{(1)} = \boldsymbol{0}, \quad \boldsymbol{y}\in\partial\omega_s. \end{gather}$$

We take the intrinsic average of (3.21a) and apply the divergence theorem to the right-hand side which vanishes by (3.21b). Then, applying the transport theorem (A12), derived in Appendix A, to the left-hand side of the intrinsic average of (3.21a) yields

(3.22)\begin{equation} \boldsymbol{\nabla}_x\boldsymbol{\cdot}\int_{\omega_f}\boldsymbol{v}^{(0)}\,\mathrm{d}S_y = 0, \end{equation}

where we have used (3.14c). Expressing (3.22) in terms of the averaged quantity $\boldsymbol {V}^{(0)}$ gives

(3.23)\begin{equation} \hat{\boldsymbol{\nabla}}\boldsymbol{\cdot}(\phi\boldsymbol{V}^{(0)}) = 0, \end{equation}

which closes the system defined in (3.20a). Similarly to (3.20a), there is no $\boldsymbol {y}$-dependence in (3.23), so we have replaced $\boldsymbol {\nabla }_x$ with $\hat {\boldsymbol {\nabla }}$.

To evaluate $\boldsymbol{\mathsf{K}}$, we find it convenient to map the system (3.17) to the naive microscale coordinate $\boldsymbol {Y}$, defined in (3.5). In the naive microscale coordinate, different values of $a$ manifest as physical changes to the domain rather than as changes to the governing equations, yielding more intuitive cell problems. This mapping gives

(3.24a)$$\begin{gather} \boldsymbol{\mathsf{I}}-\boldsymbol{\nabla}_Y\otimes\boldsymbol{\varPi}^\star{+}\nabla_Y^2\boldsymbol{{\mathcal{K}}}^\star{=} \boldsymbol{0}, \quad {\boldsymbol{Y}}\in\omega_f^{{\star}}(x_1), \end{gather}$$

(3.24b)$$\begin{gather}\boldsymbol{\nabla}_Y\boldsymbol{\cdot}{\boldsymbol{{\mathcal{K}}}}^\star{=} \boldsymbol{0}, \quad {\boldsymbol{Y}}\in\omega_f^{{\star}}(x_1), \end{gather}$$

(3.24c)$$\begin{gather}\boldsymbol{{\mathcal{K}}}^\star{=} \boldsymbol{0}, \quad {\boldsymbol{Y}}\in \partial\omega_s^{{\star}}(x_1), \end{gather}$$

with

(3.24d)\begin{equation} \mathcal{K}^\star_{ij}:=\left(\boldsymbol{\mathcal{K}}^\star\right)_{ij} \quad \text{and} \quad \varPi_{i}^\star:=\left(\boldsymbol{\varPi}^\star\right)_{i} \quad \text{periodic on} \quad \boldsymbol{Y}\in\partial\omega_{=}^{{\star}} \quad \text{and}\quad \partial\omega_{\|}^{{\star}}, \end{equation}

where $\boldsymbol {\nabla }_Y$ is the gradient operator with respect to the coordinate $\boldsymbol {Y}$. Note that

(3.24e)\begin{equation} \left(\boldsymbol{\nabla}_Y\otimes\boldsymbol{\varPi}^\star\right)_{ij} = \frac{\partial \varPi_j^\star}{\partial Y_i} \quad \text{and} \quad(\boldsymbol{\nabla}_Y\boldsymbol{\cdot}{\boldsymbol{{\mathcal{K}}}}^\star)_{i} = \frac{\partial {\mathcal{K}}^\star_{ji}}{\partial Y_j}. \end{equation}

In § 4 we consider a porous medium with a simple, prescribed microstructure. In that section, we solve (3.24) using COMSOL Multiphysics, graphically present $\boldsymbol{\mathsf{K}}(\phi, a)$ and discuss its implications.

3.3. Transport problem

We now perform a similar homogenisation procedure for the solute-transport problem (2.7). The main difference between the classic homogenisation procedure and our procedure here is that we use the transformed microscale $\boldsymbol {y}$ to convert a locally periodic tessellating cell structure into a strictly periodic tessellating cell structure. As such, we proceed following the framework of Chapman & McBurnie (Reference Chapman and McBurnie2011) and Richardson & Chapman (Reference Richardson and Chapman2011). A key step is to consider the unit normal $\hat {\boldsymbol {n}}_s$ that appears in (2.7b). In general, under the microscale transformation (3.4a,b), $\hat {\boldsymbol {n}}_s$ will not be transformed to the geometric normal of the transformed cell. Hence, we must take care when transforming the normal into multiple-scales form.

Under the multiple-scales framework, the unit normal to the solid interface is written as a function of both the macro- and microscales $\hat {\boldsymbol {n}}_s(\hat {\boldsymbol {x}}) = \boldsymbol {n}_s(\boldsymbol {x},\boldsymbol {y})$, and similarly for the function $\hat {f}_s(\hat {\boldsymbol {x}}) = f_s(\boldsymbol {x},\boldsymbol {y})$, which vanishes on the solid interface. The consistent transformation of $\boldsymbol {n}_s \equiv n_s^{i} \boldsymbol {e}_i$ requires the consistent application of the MMS derivative transformation (3.6) to the definition of $\boldsymbol {n}_s$ in terms of $f_s$ given by (2.8), to obtain the transformed unit normal

(3.25a)\begin{equation} \boldsymbol{n}_s = \frac{\left(\boldsymbol{\nabla}_y^a + \epsilon \boldsymbol{\nabla}_x\right)f_s}{\left|\left(\boldsymbol{\nabla}_y^a + \epsilon \boldsymbol{\nabla}_x \right)f_s\right|} = \frac{\displaystyle \left(\sigma_{ij} \frac{\partial f_s}{\partial y_j} + \epsilon \frac{\partial f_s}{\partial x_i}\right) \boldsymbol{e}_i }{\left[\displaystyle \sigma_{kl}\sigma_{km} \frac{\partial f_s}{\partial y_l}\frac{\partial f_s}{\partial y_m} \right]^{1/2} + O(\epsilon)}. \end{equation}

It will also be helpful to define the leading-order transformed unit normal $\boldsymbol {n}^{Y} = n^Y_i \boldsymbol {e}_i$ as follows:

(3.25b)\begin{equation} \boldsymbol{n}^{Y} = \frac{\displaystyle \sigma_{ij} \frac{\partial f_s}{\partial y_j} \boldsymbol{e}_i}{\left[\displaystyle \sigma_{kl}\sigma_{km} \frac{\partial f_s}{\partial y_l}\frac{\partial f_s}{\partial y_m} \right]^{1/2}}, \end{equation}

such that $\boldsymbol {n}_s \sim \boldsymbol {n}^{Y}$ as $\epsilon \to 0$. However, we also note that the geometric unit normal $\boldsymbol {n}^y = n^y_i \boldsymbol {e}_i$ is defined as

(3.25c)\begin{equation} \boldsymbol{n}^y := \frac{\boldsymbol{\nabla}_y\, f_s}{\left|\boldsymbol{\nabla}_y\, f_s \right|} = \frac{\displaystyle \frac{\partial f_s}{\partial y_i} \boldsymbol{e}_i}{\left[\displaystyle \frac{\partial f_s}{\partial y_j} \frac{\partial f_s}{\partial y_j} \right]^{1/2}}, \end{equation}

where $\boldsymbol {\nabla }_y$ is the gradient operator with respect to the coordinate $\boldsymbol {y}$. Importantly, the transformed normal (3.25a) and geometric normal (3.25c) are not equal. Moreover, comparing (3.25b) and (3.25c) reveals that they are not even equal to leading order in $\epsilon$ unless $a \equiv 1$.

To facilitate our subsequent manipulation of the transformed problem, we write the transformed normal $\boldsymbol {n}_s$ in terms of the geometric normal $\boldsymbol {n}^y$. Since (3.25c) can be rearranged to obtain $\partial f_s /\partial y_i = |\boldsymbol {\nabla }_y\, f_s | n^y_i$, we can re-write the transformed normal (3.25a) as

(3.25d)\begin{equation} \boldsymbol{n}_s = \frac{\displaystyle \left(\sigma_{ij} n^y_j + \epsilon N_i\right) \boldsymbol{e}_i}{\left[\displaystyle \sigma_{kl}\sigma_{km} n^y_ln^y_m \right]^{1/2} + O(\epsilon)}, \end{equation}

where the macroscale perturbation to the normal $\boldsymbol {N} = N_i \boldsymbol {e}_i$ is defined as

(3.25e)\begin{equation} \boldsymbol{N}:=\frac{\boldsymbol{\nabla}_x\,f_s}{|\boldsymbol{\nabla}_y\, f_s|}. \end{equation}

The macroscale perturbation to the normal $\boldsymbol {N}$ formally quantifies the effect of the transformed microscale structure varying over the macroscale within the MMS framework.

Having defined the transformed normal in terms of the geometric normal, we are now in a position to proceed with the homogenisation. Under the spatial transformations (3.6), (2.7) become

(3.26a)\begin{align} &\epsilon\frac{\partial {{c}}}{\partial {t}} = \left(\epsilon \frac{\partial}{\partial x_i}+\sigma_{ij} \frac{\partial}{\partial y_j}\right) \left[\frac{\partial c}{\partial x_i}+\dfrac{\sigma_{ik}}{\epsilon}\frac{\partial c}{\partial y_k}- Pe \ {v_i}{c}\right], \quad {{\boldsymbol{y}}}\in\omega_f(x_1), \end{align}

(3.26b)\begin{align} &-\epsilon\gamma{c} \left[\displaystyle \sigma_{kl}\sigma_{km} n^y_ln^y_m \right]^{1/2} +O(\epsilon^2)\notag\\ &\quad = (\sigma_{ij} n^y_j + \epsilon N_i) \left[\frac{\partial c}{\partial x_i}+\dfrac{\sigma_{ik}}{\epsilon}\frac{\partial c}{\partial y_k}- Pe \ {v_i}{c}\right], \quad {{\boldsymbol{y}}}\in \partial\omega_s(x_1), \end{align}

with

(3.26c)\begin{equation} v_i, \quad c, \quad \text{periodic on} \ \boldsymbol{y}\in\partial\omega_{=}\quad \text{and}\quad \partial\omega_{\|}, \end{equation}

where $\boldsymbol {v} = v_i \boldsymbol {e}_i$. Note that, for clarity of presentation in what follows, we have multiplied by $\epsilon$ when deriving (3.26a) from (2.7a). We now consider an expansion of the concentration field of the form

(3.27)\begin{equation} c(\boldsymbol{x}, \boldsymbol{y},t) = c^{(0)}(\boldsymbol{x}, \boldsymbol{y},t) + \epsilon c^{(1)}(\boldsymbol{x}, \boldsymbol{y},t) + \epsilon^2 c^{(2)}(\boldsymbol{x}, \boldsymbol{y},t) +\cdots \quad \text{as} \ \epsilon \to 0. \end{equation}

Note that we take $Pe, \gamma = O(1)$ to be constants independent of $\epsilon$, corresponding to a distinguished limit where all transport mechanisms balance over the macroscale (cf. (3.45c)). Considering (3.26) at leading order – that is, $O(1/\epsilon )$ – we obtain

(3.28a)$$\begin{gather} \sigma_{ij} \sigma_{ik} \dfrac{\partial^2 c^{(0)}}{\partial y_j\partial y_k} = 0, \quad \boldsymbol{y} \in \omega_f(x_1), \end{gather}$$

(3.28b)$$\begin{gather}\sigma_{ij} \sigma_{ik} n_j^y \dfrac{\partial c^{(0)}}{\partial y_k} = 0, \quad \boldsymbol{y} \in \partial \omega_s({x_1}), \end{gather}$$

(3.28c)$$\begin{gather}{v}^{(0)}_i, \quad c^{(0)}\quad \text{periodic on} \quad \boldsymbol{y}\in\partial\omega_{=} \quad \text{and}\quad \partial\omega_{\|}. \end{gather}$$

By inspection, we find that $c^{(0)} = c^{(0)}(\boldsymbol {x},t)$ is a non-trivial solution to this system. By linearity, this solution is unique and therefore the leading-order concentration is independent of $\boldsymbol {y}$.

Considering (3.26) at $O(1)$, we obtain

(3.29a)$$\begin{gather} \sigma_{ij} \sigma_{ik} \dfrac{\partial^2 c^{(1)}}{\partial y_j\partial y_k} = 0, \quad \boldsymbol{y} \in \omega_f(x_1), \end{gather}$$

(3.29b)$$\begin{gather}\sigma_{ij} \sigma_{ik} n_j^y \dfrac{\partial c^{(1)}}{\partial y_k} ={-}\sigma_{ij}n_j^y \dfrac{\partial c^{(0)}}{\partial x_i}, \quad \boldsymbol{y} \in \partial \omega_s({x_1}), \end{gather}$$

(3.29c)$$\begin{gather}{v}^{(1)}_i, \quad c^{(1)} \quad \text{periodic on} \quad \boldsymbol{y}\in\partial\omega_{=}\quad\text{and}\quad \partial\omega_{\|}, \end{gather}$$

where we have used $c^{(0)} = c^{(0)}(\boldsymbol {x},t)$, microscale incompressibility (3.14b) and the no-slip and no-penetration conditions (3.14c) on the solid surface. The form of (3.29) suggests that we can scale $\boldsymbol {\nabla }_x c^{(0)}$ out of the problem via the substitution

(3.30)\begin{equation} c^{(1)}(\boldsymbol{x}, \boldsymbol{y},t) ={-} \dfrac{\partial c^{(0)}}{\partial x_n} \varGamma_n (\boldsymbol{x}, \boldsymbol{y}) + \breve{c}(\boldsymbol{x},t), \end{equation}

where $\breve {c}$ is a scalar function and the functions $\varGamma _n$ satisfy the following cell problems:

(3.31a)$$\begin{gather} \sigma_{ij} \sigma_{ik} \dfrac{\partial^2 \varGamma_n}{\partial y_j\partial y_k} = 0, \quad \boldsymbol{y} \in \omega_f(x_1), \end{gather}$$

(3.31b)$$\begin{gather}\sigma_{ij} \sigma_{ik} n_j^y \dfrac{\partial \varGamma_n}{\partial y_k} = \sigma_{ni}n_i^y, \quad \boldsymbol{y} \in \partial \omega_s({x_1}), \end{gather}$$

(3.31c)$$\begin{gather}\varGamma_n \quad \text{periodic on} \quad \boldsymbol{y}\in\partial\omega_{=}\quad \text{and}\quad \partial\omega_{\|}. \end{gather}$$

Note that we enforce

(3.31d)\begin{equation} \langle\varGamma_n\rangle = 0, \end{equation}

which uniquely defines $\varGamma _n$. Equations (3.31) are obtained by substituting (3.30) into (3.29). Equations (3.31) must then be solved numerically for $n\in \{1,2\}$ and each desired cell geometry (i.e. pairs of $a$ and $\lambda$). Note that (3.31) are independent of $\boldsymbol {\nabla }_x c^{(0)}$, justifying our scalings in (3.30).

The goal of this analysis remains to determine a macroscale equation for the concentration. Since there are no macroscopic transport mechanisms present at this order, there is not enough information to determine a macroscale governing equation for the concentration. Hence, we must proceed to the next order in (3.26), which yields

(3.32a)$$\begin{gather} \frac{\partial c^{(0)}}{\partial t} = \sigma_{ij}\dfrac{\partial \mathcal{A}_i}{\partial y_j} + \dfrac{\partial \mathcal{B}_i}{\partial x_i}, \quad \boldsymbol{y} \in \omega_f(x_1), \end{gather}$$

(3.32b)$$\begin{gather}-\gamma c^{(0)} \left[\sigma_{kl}\sigma_{km} n^y_ln^y_m \right]^{1/2} = \sigma_{ij}n_j^y \mathcal{A}_i + N_i \mathcal{B}_i \quad \boldsymbol{y} \in \partial \omega_s({x_1}), \end{gather}$$

(3.32c)$$\begin{gather}{v}^{(2)}_i,\quad c^{(2)} \quad\text{periodic on} \quad \boldsymbol{y}\in\partial\omega_{=}\quad \text{and}\quad \partial\omega_{\|}, \end{gather}$$

where

(3.32d)$$\begin{gather} \mathcal{A}_i := \sigma_{ij} \dfrac{\partial c^{(2)}}{\partial y_j} + \dfrac{\partial c^{(1)}}{\partial x_i} - Pe \ \left( v^{(0)}_i c^{(1)} + v^{(1)}_i c^{(0)} \right), \end{gather}$$

(3.32e)$$\begin{gather}\mathcal{B}_i := \sigma_{ij} \dfrac{\partial c^{(1)}}{\partial y_j} + \dfrac{\partial c^{(0)}}{\partial x_i} - Pe \ v^{(0)}_i c^{(0)}. \end{gather}$$

Integrating (3.32a) over the transformed microscale fluid domain $\omega _f$ gives

(3.33)\begin{equation} |\omega_f| \frac{\partial c^{(0)}}{\partial t} = \int_{\omega_f} \, \sigma_{ij}\dfrac{\partial \mathcal{A}_i}{\partial y_j} \, \mathrm{d}S_y + \int_{\omega_f} \, \dfrac{\partial \mathcal{B}_i}{\partial x_i} \, \mathrm{d}S_y. \end{equation}

Applying the divergence theorem to the first integral on the right-hand side of (3.33) yields

(3.34)\begin{equation} \int_{\omega_f} \, \sigma_{ij}\dfrac{\partial \mathcal{A}_i}{\partial y_j} \, \mathrm{d}S_y = \int_{\partial\omega_s} \, \sigma_{ij}n^y_j \mathcal{A}_i\, \mathrm{d}s_y+\int_{\partial\omega} \, \sigma_{ij}n^{\square}_j \mathcal{A}_i\, \mathrm{d}s_y, \end{equation}

where $\mathrm {d}s_y$ signifies an element of a scalar line integral and $\boldsymbol {n}^{\square } = n^{\square }_j \boldsymbol {e}_j$ is the outward-facing unit normal to the external square boundary $\partial \omega$. Since $\mathcal {A}_i$ is periodic on $\partial \omega$, the last term on the right-hand side of (3.34) vanishes. Then, using (3.32b), we may re-write (3.34) as

(3.35)\begin{equation} \int_{\omega_f} \, \sigma_{ij}\dfrac{\partial \mathcal{A}_i}{\partial y_j} \, \mathrm{d}S_y ={-}\int_{\partial\omega_s} \, N_i \mathcal{B}_i \, \mathrm{d}s_y -\int_{\partial\omega_s} \, \gamma c^{(0)}\left[\sigma_{kl}\sigma_{km} n^y_ln^y_m \right]^{1/2} \, \mathrm{d}s_y. \end{equation}

To manipulate the final integral on the right-hand side of (3.33), we apply the transport theorem (A12)

(3.36)\begin{equation} \int_{\omega_f} \, \dfrac{\partial \mathcal{B}_i}{\partial x_i} \, \mathrm{d}S_y = \dfrac{\partial}{\partial x_i} \int_{\omega_f} \, \mathcal{B}_i \, \mathrm{d}S_y + \int_{\partial \omega_s}\, N_i \mathcal{B}_i \,\mathrm{d}s_y. \end{equation}

Thus, combining (3.33), (3.35) and (3.36) we obtain

(3.37)\begin{align} |\omega_f| \frac{\partial c^{(0)}}{\partial t} &= \dfrac{\partial}{\partial x_i} \int_{\omega_f} \, \left[\sigma_{ij} \dfrac{\partial c^{(1)}}{\partial y_j} + \dfrac{\partial c^{(0)}}{\partial x_i} - Pe \ v^{(0)}_i c^{(0)}\right] \, \mathrm{d}S_y\notag\\ &\quad - \gamma c^{(0)} \int_{\partial\omega_s} \,\left[\sigma_{kl}\sigma_{km} n^y_ln^y_m \right]^{1/2} \, \mathrm{d}s_y. \end{align}

Using the definitions of $c^{(1)}$ (3.30) and $\boldsymbol {V}^{(0)}$ (3.18a) and dividing through by $|\omega _f| = \phi$, we can re-write (3.37) as

(3.38a)\begin{equation} \frac{\partial C^{(0)}}{\partial t}= \dfrac{1}{|\omega_f| }\dfrac{\partial}{\partial \hat{x}_i}\left[|\omega_f| {\mathsf{D}}_{ij}(\phi,a) \dfrac{\partial C^{(0)}}{\partial \hat{x}_j} - Pe \ |\omega_f| V^{(0)}_i C^{(0)} \right] - \gamma F(\phi,a) C^{(0)}, \end{equation}

where we have expanded the intrinsic concentration $C$ in powers of $\epsilon$

(3.38b)\begin{equation} C(\boldsymbol{x},t) = C^{(0)}(\boldsymbol{x},t) + \epsilon C^{(1)}(\boldsymbol{x}, t) + \epsilon^2 C^{(2)}(\boldsymbol{x},t) +\cdots \quad \text{as} \ \epsilon \to 0, \end{equation}

and have noted that $C^{(0)}= c^{(0)}$. Note also that we have replaced $x_i$ with $\hat {x}_i$ in (3.38a) since it has been averaged over the microscale and is thus independent of $\boldsymbol {y}$. In (3.38a), the components of the effective diffusivity tensor ${\mathsf{D}}_{ij}(\phi, a)$ are defined as

(3.38c)\begin{equation} {\mathsf{D}}_{ij}(\phi, a) := \delta_{ij} - \frac{1}{|\omega_f|} \int_{\omega_f} \, \sigma_{ij} \dfrac{\partial \varGamma_j}{\partial y_j} \, \mathrm{d}S_y = \delta_{ij} - \frac{1}{\phi} \int_{\omega_f} \, \sigma_{ij} \dfrac{\partial \varGamma_j}{\partial y_j} \, \mathrm{d}S_y, \end{equation}

where $\delta _{ij}$ is the Kronecker delta. The effective adsorption strength $F(\phi, a)$ is defined as

(3.38d)\begin{equation} F(\phi,a) := \dfrac{1}{|\omega_f|} \int_{\partial\omega_s} \!\left[\sigma_{kl}\sigma_{km} n^y_ln^y_m \right]^{1/2} \, \mathrm{d}s_y. \end{equation}

Hence, our homogenised transport equation is given by (3.38a), (3.38c) and (3.38d). In order to interpret the coefficients in this equation physically and evaluate them numerically, we now transform our coefficients into the naive microscale.

3.3.1. Transforming into the naive microscale coordinate

To interpret the rate $F(\phi,a)$ physically, it is helpful to map its definition (3.38d) to the naive microscale coordinate $\boldsymbol {Y}$, defined in (3.5), in a similar way to Richardson & Chapman (Reference Richardson and Chapman2011). Firstly, consider an arbitrary vector function $\boldsymbol {z}(\boldsymbol {x},\boldsymbol {y},t) = \boldsymbol {z}^\star (\boldsymbol {x},\boldsymbol {Y},t)$, such that $\boldsymbol {z} = z_i\boldsymbol {e}_i$ and $\boldsymbol {z}^\star = z_i^\star \boldsymbol {e}_i$ then (3.8), with $Z = \boldsymbol {\nabla }_y^a\boldsymbol {\cdot }\boldsymbol {z} \equiv \boldsymbol {\nabla }_Y\boldsymbol {\cdot }\boldsymbol {z}^\star = Z^\star$ gives

(3.39)\begin{equation} \frac{1}{|\omega(x_1)|}\int_{\omega(x_1)} \, \sigma_{ij}\frac{\partial z_i}{\partial y_j}, \mathrm{d} S_y \equiv \frac{1}{|\omega^{{\star}}(x_1)|}\int_{\omega^{{\star}}(x_1)} \ \frac{\partial z^\star_i}{\partial Y_i} \, \mathrm{d}S_Y. \end{equation}

Applying the divergence theorem to both sides of (3.39) leads to the relation

(3.40)\begin{equation} \dfrac{1}{|\omega|}\int_{\partial\omega_s} \,\sigma_{ij} z_j n^y_i \, \mathrm{d}s_y = \dfrac{1}{|\omega^{{\star}}|} \int_{\partial\omega_s^{{\star}}} \!{n^Y_i}^\star z^\star_i \, \mathrm{d}s_Y, \end{equation}

where ${n}_i^Y(\boldsymbol {y}) = {n_i^{Y}}^\star (\boldsymbol {Y})$. Note that as $\boldsymbol {\sigma }$ is diagonal $\sigma _{ij} z_j n^y_i = \sigma _{ij} n^y_jz_i$. Additionally, (3.25) lead to the relation

(3.41)\begin{equation} \sigma_{ij} n^y_j = n^{Y}_{i} \left[\sigma_{kl} \sigma_{km} n_l^y n_m^y \right]^{1/2}. \end{equation}

Thus, setting $z_i=n_i^Y$ gives

(3.42)\begin{equation} \int_{\partial\omega_s} \left[\displaystyle \sigma_{kl}\sigma_{km} n^y_ln^y_m \right]^{1/2}\,\mathrm{d}s_y =\int_{\partial\omega_s} \sigma_{ij} n^y_jn_i^Y\,\mathrm{d}s_y = \frac{|\omega|}{|\omega^\star|}\int_{\partial\omega_s^\star}\, \mathrm{d}s_Y = \frac{|\partial\omega_s^\star|}{|\omega^\star|}, \end{equation}

since $\boldsymbol {n}^Y\boldsymbol {\cdot }\boldsymbol {n}^Y = \boldsymbol {n}^{Y^\star }\boldsymbol {\cdot }\boldsymbol {n}^{Y^\star }= 1$ and $|\omega |=1$. Hence,

(3.43)\begin{equation} F(\phi,a) = \frac{|\partial\omega_s^\star|}{|\omega_f\|\omega^\star|}\equiv \frac{|\partial\omega_s^\star|}{|\omega_f^\star|}, \end{equation}

and we deduce that $F$ represents the obstacle perimeter within a cell normalised by the fluid area within a cell.

Additionally, in order to evaluate ${\mathsf{D}}_{ij}$ we find it convenient to map the system (3.31) to the naive microscale coordinate $\boldsymbol {Y}$, defined in (3.5). This mapping transforms the cell problems (3.31) to

(3.44a)$$\begin{gather} \nabla^2_{Y} \varGamma_n^\star{=} 0, \quad \boldsymbol{Y} \in \omega_f^{{\star}}(x_1), \end{gather}$$

(3.44b)$$\begin{gather}{\boldsymbol{n}^Y}^\star \boldsymbol{\cdot} \nabla_{Y} \varGamma_n^\star{=} {n_k^Y}^\star, \quad \boldsymbol{Y} \in \partial \omega_s^{{\star}}({x_1}), \end{gather}$$

(3.44c)$$\begin{gather}\varGamma_n^\star \quad \text{periodic on} \ \boldsymbol{Y}\in\partial\omega_{=}^{{\star}}\quad \text{and}\quad \partial\omega_{\|}^{{\star}}, \end{gather}$$

with

(3.44d)\begin{equation} \langle\varGamma_n^\star\rangle=0. \end{equation}

The components of the effective diffusivity tensor (3.38c) become

(3.45a)\begin{equation} {\mathsf{D}}_{ij}(\phi, a) = \delta_{ij} - \frac{1}{|\omega_f^{{\star}}|} \int_{\omega_f^{{\star}}} \, \dfrac{\partial \varGamma_j^\star}{\partial Y_i} \, \mathrm{d}S_Y = \delta_{ij} - \frac{1}{a \phi} \int_{\omega_f^{{\star}}} \! \dfrac{\partial \varGamma_j^\star}{\partial Y_i} \, \mathrm{d}S_Y. \end{equation}

We can also write (3.45a) in tensor form $\boldsymbol{\mathsf{D}}$ as

(3.45b)\begin{equation} \boldsymbol{\mathsf{D}}(\phi, a) = \boldsymbol{\mathsf{I}} - \frac{1}{|\omega_f^{{\star}}|} \int_{\omega_f^{{\star}}} \boldsymbol{\nabla}_{Y}\otimes{\boldsymbol{\varGamma}}^\star \, \mathrm{d}S_Y. \end{equation}

To evaluate $\boldsymbol{\mathsf{D}}$, we solve the transformed cell problems (3.44) numerically in COMSOL Multiphysics.

Using the results from this subsection, we may re-write (3.38a) in vector/tensor form as

(3.45c)\begin{equation} \frac{\partial C^{(0)}}{\partial t}= \dfrac{1}{\phi} \hat{\nabla} \boldsymbol{\cdot} \left(\phi \boldsymbol{\mathsf{D}} \hat{\nabla} C^{(0)} - Pe \ \phi \boldsymbol{V}^{(0)} C^{(0)}\right) - \gamma \dfrac{|\partial\omega_s^{{\star}}|}{|\omega_f^{{\star}}|} C^{(0)}. \end{equation}

Again, since there is no $\boldsymbol {y}$-dependence in (3.45c), we have replaced $\boldsymbol {\nabla }_x$ with $\hat {\boldsymbol {\nabla }}$. Equation (3.45c) describes macroscopic transport by advection and diffusion in a porous medium with chemical sorption, where $\phi \boldsymbol{\mathsf{D}}\boldsymbol {\cdot }\hat {\boldsymbol {\nabla }}C^{(0)}$ is the diffusive flux per unit area of porous medium and $\boldsymbol{\mathsf{D}}\boldsymbol {\cdot }\hat {\boldsymbol {\nabla }}C^{(0)}$ is the diffusive flux per unit area of fluid. The form of the effective macroscale transport equation (3.45c) is similar to that obtained in Dalwadi et al. (Reference Dalwadi, Griffiths and Bruna2015), where a simpler problem with a constant cell size is considered, resulting in a more straightforward upscaling procedure. Here, we have formally accounted for a slowly varying cell size and a slowly varying microscale geometry. The most significant difference between structure of the macroscale (3.45c) and the equivalent result obtained from the classic homogenisation of a strictly periodic problem arises from the slowly varying microscale geometry, which manifests through the explicit (non-trivial) porosity dependence of the diffusive term (cf. Bruna & Chapman Reference Bruna and Chapman2015). We find that the effect of the slowly varying cell size is less important to the structural form of the derived macroscale equations, which was not known at the outset.

Here, we have formally accounted for slowly varying cell size and slowly varying obstacle size. The resulting model has the same form as that derived in previous work for uniform cell size and slowly varying obstacle size (cf. Bruna & Chapman Reference Bruna and Chapman2015), in the sense that there are three macroscopic coefficients $\boldsymbol{\mathsf{K}}$, $\boldsymbol{\mathsf{D}}$ and $F$ that vary with the local microstructure via the values of $\phi$, $R$ and here also $a$. Allowing for slowly varying cell size has not otherwise altered the mathematical structure of the macroscale problem at leading order, suggesting that different types of microscale heterogeneity can lead to a similar mathematical structure on the macroscale. The key difference between these heterogeneous results and the classical result for a uniform microstructure is that factors of porosity appear in front of the time derivative and within the divergence, the latter multiplying the macroscopic solute flux.

In § 4, we calculate and discuss the permeability, effective diffusivity and effective adsorption strength for a porous medium with a simple, prescribed microstructure. Note that although our model problem of a one-dimensional filter in § 4.2 features flow in the longitudinal direction, the macroscopic flow and transport equations (3.20), (3.23) and (3.45c) and the results in § 4.1 are valid for any arbitrary flow direction.

4.1. Macroscale flow and transport properties

For this specific geometry we explore the impact of microstructure on macroscopic flow and transport by analysing the permeability and effective net diffusivity tensors, $\boldsymbol{\mathsf{K}}$ and $\phi \boldsymbol{\mathsf{D}}$, respectively. To determine $\boldsymbol{\mathsf{K}}$ we solve (3.24) in COMSOL Multiphysics using the ‘Laminar Flow (spf)’ interface (‘Fluid Flow’ $\to$ ‘Single Phase Flow’ $\to$ ‘Laminar flow (spf)’). The domain is discretised using the ‘Physics-controlled mesh’ with the element size set to ‘Extremely fine’. Similarly, to evaluate $\boldsymbol{\mathsf{D}}$, we solve (3.44), in COMSOL Multiphysics using the ‘Laplace Equation (lpeq)’ interface (‘Classical PDEs’ $\to$ ‘Mathematics branch’ $\to$ ‘Laplace Equation (lpeq)’). For the flow problem, the domain is discretised using the ‘Physics-controlled mesh’ with the element size set to ‘Extremely fine.’

The tensors $\boldsymbol{\mathsf{K}}$ and $\boldsymbol{\mathsf{D}}$ depend on microstructure via $a$, $\phi$ and $R$, any two of which are independent and the third prescribed by (4.1) (figure 4). We therefore have one additional degree of microstructural freedom relative to Dalwadi et al. (Reference Dalwadi, Griffiths and Bruna2015) and this allows us to explore the anisotropy in the system. We explore the effect of the $a$, $R$ and $\phi$ parameter space on $\boldsymbol{\mathsf{K}}$ and $\phi \boldsymbol{\mathsf{D}}$ in figures 5 and 6, respectively. The effective diffusivity $\boldsymbol{\mathsf{D}}$ is shown for reference in figure 10a–d (Appendix B).

Figure 4.

The porosity $\phi$ of a rectangular cell increases with aspect ratio $a$ and decreases with obstacle radius $R$ according to (4.1). (a) Aspect ratio $a$ vs $\phi$ for $R\in \{0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5\}$ (dark to light). The attainable region of the $\phi$-$a$ plane (shaded grey) is bounded below by $\phi _{min}$ given by (3.9) with $a=2R$ (dot-dashed line). The cell geometry for two distinct points with $R= 0.25$ is shown pictorially. (b) Aspect ratio $a$ vs $R$ for $\phi \in \{0.35, 0.45, 0.55, 0.65, 0.75, 0.85, 0.925, 0.975, 0.995, 0.999\}$ (dark to light). The attainable region of the $R$–$a$ plane (shaded grey) is also bounded below by $a_{min}:=2R$ (dot-dashed line). Note that the smallest attainable $\phi$ for any $R$, $a$ combination is $\phi _{min}(R=1/2)=1-{\rm \pi} /4$.

Figure 5.

The longitudinal permeability ${{\mathsf{K}}}_{11}$ (a,b), transverse permeability ${{\mathsf{K}}}_{22}$ (c,d) and the permeability–anisotropy ratio ${{\mathsf{K}}}_{22}/{{\mathsf{K}}}_{11}$ (e, f) depend strongly on microstructure. Panels (a,c,e) show ${{\mathsf{K}}}_{11}$, ${{\mathsf{K}}}_{22}$ and ${{\mathsf{K}}}_{22}/{{\mathsf{K}}}_{11}$ against $\phi$ for fixed values of $R\in \{0.1, 0.2, 0.3, 0.4, 0.49\}$, with $a$ varying according to (4.1). (b,d, f) Values of the same quantities against $R$ for fixed values of $\phi \in \{0.65, 0.7, 0.8, 0.9, 0.95\}$, with $a$ varying according to (4.1). For a given value of $R$, the minimum porosity $\phi _{min}(R)$ is given by (4.1) with $a=2R$. Note that ${{\mathsf{K}}}_{11}$ is non-zero at $\phi _{min}(R)$ for all values of $\phi$ (dot-dashed curve, (a,b)), whereas ${{\mathsf{K}}}_{22}$ vanishes at $\phi _{min}$ ((c,e) red vertical asymptotes). There exists a smallest possible $R$ for any given $\phi$ ((d, f) red vertical asymptotes). In all cases, ${{\mathsf{K}}}_{11}$ and ${{\mathsf{K}}}_{22}$ are as defined in (3.20b) and calculated using COMSOL Multiphysics. The permeability is isotropic when $a\equiv 1$ (solid black curves; Dalwadi et al. Reference Dalwadi, Griffiths and Bruna2015). The limit $R\to 0$ and $a \to 0$ corresponds to a set of parallel but disconnected channels with unit transverse width, for which $\phi \to 1$, ${{\mathsf{K}}}_{22}\to 0$ and ${{\mathsf{K}}}_{11}\to 1/12$ (black diamonds in a,b).

Figure 6.

The net longitudinal diffusivity $\phi {{\mathsf{D}}}_{11}$ (a,b), net transverse diffusivity $\phi {{\mathsf{D}}}_{22}$ (c,d) and the diffusivity–anisotropy ratio ${{\mathsf{D}}}_{22}/{{\mathsf{D}}}_{11}$ (e, f) depend strongly on microstructure. Panels (a,c,e) show $\phi {{\mathsf{D}}}_{11}$, $\phi {{\mathsf{D}}}_{22}$ and ${{\mathsf{D}}}_{22}/{{\mathsf{D}}}_{11}$ against $\phi$ for fixed values of $R\in \{0.1, 0.2, 0.3, 0.4, 0.49\}$, with $a$ varying according to (4.1). Panels (b,d, f) show the same quantities against $R$ for fixed values of $\phi \in \{0.65, 0.7, 0.8, 0.9, 0.95\}$, with $a$ varying according to (4.1). For a given value of $R$, the minimum porosity $\phi _{min}(R)$ is given by (4.1) with $a=2R$. Note that $\phi {{\mathsf{D}}}_{11}$ is non-zero at $\phi _{min}(R)$ for all values of $\phi$ (dot-dashed curve, a,b), whereas $\phi {{\mathsf{D}}}_{22}$ vanishes at $\phi _{min}$. In all cases, ${{\mathsf{D}}}_{11}$ and ${{\mathsf{D}}}_{22}$ are as defined in (3.45a) and calculated using COMSOL Multiphysics. Both ${{\mathsf{D}}}_{11}$ and ${{\mathsf{D}}}_{22}$ tend to 1 as $\phi \to 1$, which corresponds to the limit of free-space diffusion. The net effective diffusivity, $\phi \boldsymbol{\mathsf{D}}$, is isotropic when $a=1$ (solid black curves), in agreement with the results presented in Dalwadi et al. (Reference Dalwadi, Griffiths and Bruna2015).

We have validated our analysis for this geometry in a number of ways. Firstly, we have compared our results with those in Dalwadi et al. (Reference Dalwadi, Griffiths and Bruna2015) for the special case of $a\equiv 1$, confirming both the final homogenised equations ((3.20), (3.45a), (3.45c) and (4.2)) and the detailed numerical results (black lines; figures 5–7). Secondly, we have confirmed our results in the Hele-Shaw limit of parallel, disconnected channels where the longitudinal permeability is $1/12$ (black diamond; figure 5) and the transverse permeability vanishes. Finally, we have confirmed that both the transverse permeability and the transverse effective diffusivity vanish when the transverse connectivity vanishes ($a\to {}2R$ or $2R\to {}a$; red lines in figures 5 and 6), and that the permeability diverges and the effective diffusivity tends to unity as the obstacles vanish ($\phi \to 1$).

Figure 7.

The effective adsorption rate $F$ depends strongly on microstructure. (a) Effective adsorption strength $F$ against $\phi$ for fixed values of $R\in \{0.1, 0.2, 0.3, 0.4, 0.49\}$, with $a$ varying according to (4.1). (b) Effective adsorption strength $F$ against $R$ for fixed values of $\phi \in \{0.65, 0.7, 0.8, 0.9, 0.95\}$, with $a$ varying according to (4.1). In all cases, $F$ is as defined in (4.2). For a given value of $R$, the minimum porosity $\phi _{min}(R)$ is given by (4.1) with $a=2R$. The results of Dalwadi et al. (Reference Dalwadi, Griffiths and Bruna2015) are again reproduced when $a=1$ (solid black).

Increasing $\phi$ at fixed $R$ is achieved by increasing $a$ (figure 4a), such that the obstacles move further apart in the longitudinal direction only; as a result, ${{\mathsf{K}}}_{11}$, ${{\mathsf{K}}}_{22}$, $\phi {{\mathsf{D}}}_{11}$ and $\phi {{\mathsf{D}}}_{22}$ all increase (figures 5a,c and 6a,c). As $\phi \to 1$ ($a\to \infty$), both ${{\mathsf{K}}}_{11}$ and ${{\mathsf{K}}}_{22}$ diverge as the resistance to flow vanishes (figure 5a,c), and both $\phi {{\mathsf{D}}}_{11}$ and $\phi {{\mathsf{D}}}_{22}$ tend to 1 as molecular diffusion becomes unobstructed (figure 6a,c). As $\phi \to \phi _{min}(R)$ ($a\to 2R$) at fixed $R$, the obstacles move closer together in the longitudinal direction and the pore space becomes disconnected in the transverse direction, so that ${{\mathsf{K}}}_{22}$ and $\phi {{\mathsf{D}}}_{22}$ vanish; ${{\mathsf{K}}}_{11}$ and $\phi {{\mathsf{D}}}_{11}$ are minimised but do not vanish. Further taking $R\to 0$, the longitudinal problem reduces to a set of disconnected parallel channels of unit transverse width, for which ${{\mathsf{K}}}_{11}=1/12$ (figure 5a,b, black diamond).

Increasing $R$ at fixed $\phi$ is similarly achieved by increasing $a$ (figure 4b), in which case the transverse channels between obstacles grow wider while the longitudinal channels between obstacles grow narrower. As a result, ${{\mathsf{K}}}_{11}$ and $\phi {{\mathsf{D}}}_{11}$ decrease while ${{\mathsf{K}}}_{22}$ and $\phi {{\mathsf{D}}}_{22}$ increase. As $R \to 1/2$ at fixed $\phi$, the longitudinal channels close and ${{\mathsf{K}}}_{11}$ and $\phi {{\mathsf{D}}}_{11}$ vanish, but the transverse channels become wider and ${{\mathsf{K}}}_{22}$ and $\phi {{\mathsf{D}}}_{22}$ are maximised. The longitudinal permeability, ${{\mathsf{K}}}_{11}$, is weakly non-monotonic in $R$ for larger values of $\phi$ (figure 5b), which means that the longitudinal permeability of a high-porosity porous medium can be maximised for a given $\phi$ by appropriately varying $R$ and $a$.

When $a\equiv 1$, equivalent to the case considered in Dalwadi et al. (Reference Dalwadi, Griffiths and Bruna2015), $\boldsymbol{\mathsf{K}}$ and $\phi \boldsymbol{\mathsf{D}}$ become isotropic. Increasing $\phi$ corresponds to decreasing $R$, in which case both the longitudinal and transverse spacing between obstacles decreases (figure 4) which decreases ${{\mathsf{K}}}_{11}={{\mathsf{K}}}_{22}$ and $\phi {{\mathsf{D}}}_{11}=\phi {{\mathsf{D}}}_{22}$ (figures 5 and 6, solid black lines). For $a\not =1$, our microstructure is inherently anisotropic (${{\mathsf{K}}}_{11}\neq {{\mathsf{K}}}_{22}$, $\phi {{\mathsf{D}}}_{11}\neq \phi {{\mathsf{D}}}_{22}$). For $a<1$, the longitudinal channels are wider than the transverse channels, such that ${{\mathsf{K}}}_{22}/{{\mathsf{K}}}_{11}<1$ and ${{\mathsf{D}}}_{22}/{{\mathsf{D}}}_{11}<1$ and both ratios vanish as $a\to 2R$ ($\phi \to \phi _{min}$ where $\phi _{min}$ is given by (4.1); figures 5(e, f) and 6(e, f), respectively). For $a>1$, the longitudinal channels are narrower than the transverse channels, such that ${{\mathsf{K}}}_{22}/{{\mathsf{K}}}_{11}> 1$ and ${{\mathsf{D}}}_{22}/{{\mathsf{D}}}_{11}>1$. The permeability–anisotropy ratio, ${{\mathsf{K}}}_{22}/{{\mathsf{K}}}_{11}$, increases monotonically with both $\phi$ and $R$, and diverges as $\phi \to 1$ at fixed $R$ (${{\mathsf{K}}}_{22}$ diverges faster than ${{\mathsf{K}}}_{11}$ because the obstacles never get further apart in the transverse direction) and as $R\to 1/2$ at fixed $\phi$ (${{\mathsf{K}}}_{11}$ vanishes; figure 5e, f). The diffusivity–anisotropy ratio, ${{\mathsf{D}}}_{22}/{{\mathsf{D}}}_{11}$, increases monotonically with $R$ for all $\phi \in (\phi _{min},1)$, diverging as $R\to 1/2$ (figure 6f). For fixed $R$, this ratio increases monotonically with $\phi$ for $a\leq 1$, is equal to unity for $a=1$ (isotropic geometry), and must approach unity as $\phi \to 1$ ($a\to \infty$; unobstructed molecular diffusion; figure 6e). These bounds require that ${{\mathsf{D}}}_{22}/{{\mathsf{D}}}_{11}$ has an intermediate maximum in $\phi$ (or in $a$) at fixed $R$, the amplitude of which diverges as $R\to 1/2$. Specifically, the non-monotonicity in the ratio ${{\mathsf{D}}}_{22}/{{\mathsf{D}}}_{11}$ occurs due to the relative rates of increase of $\phi {{\mathsf{D}}}_{11}$ and $\phi {{\mathsf{D}}}_{22}$. For $\phi =\phi _{min}$ ($a = 2R$) there is no transverse connectivity thus separating the obstacles slightly (a small increase in $a$) leads to a sharp increase in $\phi {{\mathsf{D}}}_{22}$ but only a slight increase in $\phi {{\mathsf{D}}}_{11}$ since the longitudinal connectivity is unchanged and most longitudinal mixing occurs in the longitudinal channels. Conversely, as $\phi \to 1$ ($a\to \infty$) the longitudinal spacing between obstacles diverges which means that $\phi {{\mathsf{D}}}_{11}$ is very sensitive to changes in $a$ as longitudinal mixing occurs predominantly between longitudinally adjacent obstacles (in the transverse channels), thus $\phi {{\mathsf{D}}}_{11}$ rapidly approaches unity. However, significant transverse connectivity is preserved for large $a$ so increasing $a$ further has minimal effect on $\phi {{\mathsf{D}}}_{22}$ since most transverse mixing occurs in the transverse channels in this limit. Note that the longitudinal diffusivity ${{\mathsf{D}}}_{11}$ is non-monotonic in $\phi$ for each $R$ as $a$ varies (Appendix B, figure 10).

The partially absorbing boundary condition on the microscale, whose strength is measured by the parameter ${\gamma }$ in (2.7b), leads to an effective sink term in the macroscale transport problem, whose strength is measured by $\gamma F$, where $F$ is given in (4.2). Here, $F$ is the ratio of the perimeter of an obstacle to the fluid area within a cell, which are $2{\rm \pi} {}R$ and $a \phi$, respectively, for a rectangular array of circular obstacles. We consider the impact of microstructure on the removal of solute in more detail in § 4.2. Note that $F$ decreases as $\phi$ increases at fixed $R$, as should be expected, but also as $R$ increases at fixed $\phi$; the latter occurs because an increase in obstacle size requires a correspondingly larger increase in cell size to keep $\phi$ constant. We consider the impact of microstructure on the removal of solute in more detail in § 4.2.

4.2. Simple one-dimensional filter

We now use the homogenised model to understand the effect of microstructure and Péclet number on filter efficiency in the context of a simple one-dimensional steady-state filtration problem. We identify the performance of the filter with the rate at which it removes solute, and thus use the leading-order outlet concentration $C_{{out}}^{(0)}$ as a measure of filtration efficiency. Specifically, we consider (3.20) and (3.38) at steady state, with imposed flux and concentration at the inlet,

(4.3a)$$\begin{gather} \phi{\boldsymbol{V}^{(0)}} = \boldsymbol{e}_1 \quad \text{at} \ \hat{x}_1 = 0, \end{gather}$$

(4.3b)$$\begin{gather}C^{(0)} = 1 \quad \text{at} \ \hat{x}_1 = 0, \end{gather}$$

and passive outflow at the outlet,

(4.3c)\begin{equation} \frac{\partial C^{(0)}}{\partial \hat{x}_1} = 0 \quad \text{at} \ \hat{x}_1 = 1. \end{equation}

Since these boundary conditions (4.3a) are compatible with unidirectional flow, we take $\boldsymbol {V}^{(0)}(\hat {\boldsymbol {x}})=V_1^{(0)}(\hat {x}_1)\boldsymbol {e}_1$ and $C^{(0)}(\hat {\boldsymbol {x}},t) = C^{(0)}(\hat {{x}}_1)$. Thus, (3.23) leads to

(4.4)\begin{equation} \frac{\mathrm{d}}{\mathrm{d}\hat{x}_1} \left( \phi V_1^{(0)} \right) = 0, \end{equation}

which, on application of the inlet condition (4.3a), gives the macroscale flux $\phi {V}_1^{(0)} \equiv 1$ for all $\hat {x}_1$. The associated pressure drop across the entire filter, $\Delta P^{(0)}$, is obtained by integrating (3.20) and using the fact that $\phi {V}_1^{(0)} \equiv 1$, which gives

(4.5)\begin{equation} \Delta P^{(0)} = \int_0^1 \frac{1}{{{\mathsf{K}}}_{11}(x)}\,\mathrm{d}x. \end{equation}

Note that the right-hand side of (4.5) is a measure of the total flow resistance of the entire filter; the inverse of this quantity can be thought of as an effective permeability for the entire filter.

Hence, the homogenised governing equation for the steady concentration distribution $C^{(0)}(\hat {x}_1)$ (3.45c) becomes

(4.6)\begin{equation} \dfrac{1}{\phi}\frac{\mathrm{d}}{\mathrm{d} \hat{x}_1} \left[ \phi {{\mathsf{D}}}_{11}(\phi, R) \frac{\mathrm{d} C^{(0)}}{\mathrm{d} \hat{x}_1} - Pe \ C^{(0)} \right] = \gamma F(\phi, R) C^{(0)} \quad \text{for} \ \hat{x}_1 \in (0, 1), \end{equation}

where $F(\phi, R)$ is defined in (3.43). We solve (4.6) subject to (4.3b) and (4.3c) numerically using a finite-difference scheme. We discretise the interval [0, 1] using a uniform mesh of size $\Delta x = 1/N$ and we approximate derivatives using a second-order-accurate central-difference formula. The results presented below were obtained using $N = 500$. Below, we consider filters with varying porosity and filters with uniform porosity, but with heterogeneous microstructure in both cases.

4.2.1. Porosity gradients

We first consider filters with varying porosity. Recall that a given porosity $\phi$ may be achieved in two different ways: by fixing $a(\hat {x}_1)$ and varying $R(\hat {x}_1)$, as considered in Dalwadi et al. (Reference Dalwadi, Griffiths and Bruna2015) for $a\equiv 1$; or by fixing $R(\hat {x}_1)$ and varying $a(\hat {x}_1)$. We consider these options in figure 8, for the same three porosity fields in both cases: linearly increasing with $\hat {x}_1$, uniform in $\hat {x}_1$, or linearly decreasing with $\hat {x}_1$. Specifically, we take $\phi (\hat {x}_1) = \phi _0 + m_\phi (\hat {x}_1-0.5)$, where $\phi _0$ is the average porosity (also the mid-point porosity) and $m_\phi$ is the porosity gradient. We take $\phi _0=0.8$ and $m_\phi =0.3$ (increasing in $\hat {x}_1$), $m_\phi =0$ (uniform in $\hat {x}_1$), or $m_\phi =-0.3$ (decreasing in $\hat {x}_1$). Note that $|m_\phi |$ defines the filter microstructure and $\mathrm {sgn}(m_\phi )$ is simply the orientation of the filter.

Figure 8.

Steady-state concentration field $C^{(0)}(\hat {x}_1)$ for $Pe \in \{ 0, 1, 5\}$ (a,c) and outlet concentration $C_{out}^{(0)}:= C^{(0)}(1)$ as a function of $Pe$ (b,d). Panels (a,b) show $a\equiv 1$ and $\phi (\hat {x}_1) = 0.8 + m_\phi (\hat {x}_1 - 0.5)$. Panels (c,d) show $R\equiv 0.4$ and $\phi (\hat {x}_1) = 0.8 + m_\phi (\hat {x}_1 - 0.5)$. In both cases, $m_\phi = 0.3$ (solid), $m_{\phi } = 0$ (dotted) and $m_\phi = -0.3$ (dashed).

When varying $\phi$ by varying $R$ at fixed $a\equiv 1$, as considered by Dalwadi et al. (Reference Dalwadi, Griffiths and Bruna2015), the sign of the porosity gradient has a modest impact on the concentration distribution within the filter: $m_\phi >0$ leads to a steeper gradient in $C^{(0)}$ near the inlet and a shallower gradient in $C^{(0)}$ near the outlet, whereas $m_\phi <0$ leads to a more uniform gradient in $C^{(0)}$ throughout the filter (figure 8a). However, the outlet concentration $C_{out}^{(0)}:= C^{(0)}(1)$ is remarkably insensitive to $m_\phi$. The outlet concentration is slightly lower for $m_\phi <0$, and this slight difference decreases as $Pe$ increases (figure 8b). As $Pe$ increases, advection becomes stronger causing more solute to be swept through the filter; as a result, $C^{(0)}(\hat {x}_1)$ increases with $Pe$ for all $\hat {x}_1$ and $C^{(0)}_{out}$ more than doubles as $Pe$ increases from 0 to 10. The case when $a\equiv 1$ is considered in more detail in Dalwadi et al. (Reference Dalwadi, Griffiths and Bruna2015). Varying $\phi$ by varying $a$ at fixed $R\equiv 0.4$ leads to qualitatively similar results, but $C^{(0)}(x)$ and $C^{(0)}_{out}$ are more sensitive to $m_\phi$ (figure 8c,d). For all $Pe$, attaining a desired porosity gradient via varying $R$ leads to a more efficient filter than varying $a$, in the sense that $C_{out}^{(0)}$ is lower for the same $\phi (x_1)$.

4.2.2. Microstructural gradients with uniform porosity

We now consider filters with uniform porosity but gradients in microstructure. We therefore fix $\phi$ and simultaneously vary $R$ and $a$ with $\hat {x}_1$, recalling that $a$ is related to $R$ via (4.1). We consider two types of variation: an imposed gradient in $R$ with $a$ varying to maintain constant $\phi$ (via (4.1)) or an imposed gradient in $a$ with $R$ varying to maintain constant $\phi$ (via (4.1)). We consider these options in figure 9.

Figure 9.

Steady-state concentration field $C^{(0)}(\hat {x}_1)$ for $\phi \in \{0.7, 0.8, 0.9\}$ (a,c), and outlet concentration $C_{out}^{(0)}:= C^{(0)}(1)$ as a function of $\phi$ (b,d). Panels (a,b) show $Pe = 1$ and $R(\hat {x}_1) = 0.36 + m_R(\hat {x}_1-0.5)$ for $m_R = 0.14$ (solid), $m_{R} = 0$ (dotted), $m_R = -0.14$ (dashed). The aspect ratio $a$ varies following (4.1). For $R_0=0.36$ and $m_R=\pm 0.14$ the minimum attainable porosity is $1-0.145{\rm \pi} \approx 0.54$ ((b); black dots). Panels (c,d) show $Pe = 1$ and $a(\hat {x}_1) = 1.34 + m_a(\hat {x}_1-0.5)$ for $m_a = 1.8$ (solid), $m_{a} = 0$ (dotted) and $m_a = -1.8$ (dashed). The radius $R$ varies following (4.1). For $a_0=1.34$ and $m_a=\pm 1.8$ the minimum attainable porosity is $1-0.11{\rm \pi} \approx 0.66$ ((d); black dots).

We first consider $R(\hat {x}_1) = R_0 + m_R(\hat {x}_1-0.5)$, where $R_0$ is the average obstacle radius (also the mid-point radius) and $m_R$ is the gradient (figure 9a,b). We take $R_0 = 0.36$ and $m_R=0.14$ (increasing in $\hat {x}_1$), $m_R=0$ (uniform in $\hat {x}_1$) or $m_R=-0.14$ (decreasing in $\hat {x}_1$).

When varying $R$ linearly, the sign of $m_R$ has a modest impact on the concentration distribution: for any uniform porosity, $m_R>0$ leads to a shallower gradient in $C^{(0)}$ and a higher concentration at every point within the filter, including the outlet. Thus, for any uniform porosity, $m_R<0$ is always a more efficient filter than $m_R>0$.

We next consider $a(\hat {x}_1) = a_0 + m_a(\hat {x}_1-0.5)$, where $a_0$ is the average cell width (also the cell width at the mid-point) and $m_a$ is the gradient of $a$ over $\hat {x}_1$ (figure 9c,d). We take $a_0 = 1.34$ and $m_a=1.8$ (increasing in $\hat {x}_1$), $m_a=0$ (constant in $\hat {x}_1$) or $m_a=-1.8$ (decreasing in $\hat {x}_1$). From (4.1) for fixed $\phi$, it can be seen that $a \propto R^2$, thus, a decrease in $a$ must be mirrored by a decrease in $R$ to maintain a uniform $\phi$. Thus, we expect the same qualitative behaviour for a linear gradient in $R$ (figure 9a,b) as for a linear gradient in $a$ (figure 9c,d). We find that $m_a<0$ leads to a more efficient filter for all $\phi$.

Note that, for any $\phi \in (1-0.11{\rm \pi},1)$ – that is, the range of porosities attainable for imposed linear gradients in both $R$ and $a$ (see figure 9 caption) – prescribing $a$ and taking $m_a<0$ predicts most efficient filter considered here (comparing figures 9b and 9d). Similarly, for large $\phi \gtrsim 0.725$ prescribing $a$ and taking $m_a>0$ predicts the least efficient filter, whereas, for small $\phi \lesssim 0.725$ prescribing $R$ and taking $m_R>0$ predicts the least efficient filter.