2.1.1. Fluid incompressibility

Let us now inspect the fluid incompressibility assumption and use it in the governing equations.

We first define a volume fraction $\varphi _k$ , which simply represents the fraction of the volume occupied by the $k$ th component with respect to the total volume. From the definition, the sum of all volume fractions satisfies the condition

(2.9) \begin{align} \sum _{k=1}^{n}\varphi _k=1. \end{align}

We can now rewrite the density $\rho _k$ as

(2.10) \begin{align} \rho _k=\varphi _k\rho _k^T, \end{align}

where $\rho _k^T$ is the so-called true density of the $k$ th component of the mixture. The true density is defined as a density of the stand-alone component, that is, if it were not included in a mixture. Then, the incompressibility of a given component $k$ means that $\rho _k^T$ constant. Crucially, considering both the solvent and solute in our filtration problem to be incompressible does not result into the mixture being divergence free. That is because individual components in a mixture can rearrange in space; hence, their volume fraction volumes can change (Klika et al. Reference Klika, Whiteley, Brown and Gaffney2019). Consequently, rewriting the partial mass balances in terms of volume fractions, we obtain

(2.11) \begin{align} \frac {\partial {\varphi _k}}{\partial {t}}=-{\boldsymbol{\nabla }} \,{\boldsymbol{\cdot }}\,\left (\varphi _k\boldsymbol{u}_k\right ). \end{align}

The incompressibility of all the components of the mixture entails

(2.12) \begin{align} 0={\boldsymbol{\nabla }} \,{\boldsymbol{\cdot }}\,\left (\sum _{k=1}^{n}\varphi _k\boldsymbol{u}_k\right )\!, \end{align}

by virtue of (2.9), but which is different from vanishing divergence of (barycentric) velocity.

Rewriting partial velocities $\boldsymbol{u}_k$ in terms of diffusion fluxes $\boldsymbol{J}_k$ , we get after some algebra that the incompressibility condition (2.12) can be expressed as

(2.13) \begin{align} {\boldsymbol{\nabla }} \,{\boldsymbol{\cdot }}\,\boldsymbol{u}=-\sum _{k=1}^{n}\frac {{\boldsymbol{\nabla }} \,{\boldsymbol{\cdot }}\,\boldsymbol{J}_k}{\rho _k^T}. \end{align}

With the description of a mixture of incompressible components ready, we move to the description of the filtration problem.

2.2.1. Creeping flow of the mixture

We assume steady creeping flow of a Newtonian binary mixture, so inertia is neglected. This assumption concerns the momentum balance and does not impose $\tilde {{\boldsymbol{\nabla }}}\boldsymbol{\cdot }\tilde {\boldsymbol{u}}=0$ . Indeed, for constant shear and dilatational viscosities, the Newtonian stress may be written as

(2.15) \begin{align} \boldsymbol{\sigma }=-\pi \unicode{x1D644}_d+2\mu \textbf{D}(\tilde {\boldsymbol{u}}) +\lambda (\tilde {{\boldsymbol{\nabla }}}\boldsymbol{\cdot }\tilde {\boldsymbol{u}})\unicode{x1D644}_d. \end{align}

Its divergence contains the gradient $(\lambda +\mu )\tilde {{\boldsymbol{\nabla }}}(\tilde {{\boldsymbol{\nabla }}}\boldsymbol{\cdot }\tilde {\boldsymbol{u}})$ , which can be absorbed into a modified pressure $\tilde p=\pi -(\lambda +\mu )\tilde {{\boldsymbol{\nabla }}}\boldsymbol{\cdot }\tilde {\boldsymbol{u}}$ . The steady momentum balance therefore has the familiar Stokes form even though the barycentric velocity need not be solenoidal:

(2.16a) \begin{align} -\tilde {{\boldsymbol{\nabla }}}\tilde {p}+\mu \tilde {{{\nabla }}}^2\tilde {\boldsymbol{u}}=0, \hspace {10pt} \tilde {\boldsymbol{x}}\in \varOmega _{\kern-1pt f}. \end{align}

Here, $\tilde {p}(\tilde {\boldsymbol{x}},\tilde {t})$ stands for the pressure of the fluid and the constant $\mu$ is its viscosity. We replace the fluid mass balance equation with the derived incompressibility condition (2.13). Rewriting it for a two-component mixture using our new notation, we obtain

(2.16b) \begin{align} \tilde {{\boldsymbol{\nabla }}}\boldsymbol{\cdot }\tilde {\boldsymbol{u}}=\left (\frac {1}{\rho _v^T}-\frac {1}{\rho _u^T}\right )\tilde {{\boldsymbol{\nabla }}}\boldsymbol{\cdot }\boldsymbol{J}_u, \end{align}

where we used $\boldsymbol{J}_v=-\boldsymbol{J}_u$ as a result of (2.5) for $n=2$ . The index $v$ denotes the fluid (solvent) and the index $u$ denotes the solute. By substituting the assumed diffusion flux closure (2.7), we obtain the final form of the incompressibility condition as

(2.16c) \begin{align} \tilde {{\boldsymbol{\nabla }}}\boldsymbol{\cdot }\tilde {\boldsymbol{u}}=\left (\frac {1}{\rho _u^T}-\frac {1}{\rho _v^T}\right )\tilde {{\boldsymbol{\nabla }}}\boldsymbol{\cdot } ({\rm D}\tilde {{\boldsymbol{\nabla }}}\tilde {c}), \hspace {10pt} \tilde {\boldsymbol{x}}\in \varOmega _{\kern-1pt f}. \end{align}

It is important to note that the right-hand side of (2.16c ) is generally non-zero because we assume $\rho _u^T\neq \rho _v^T$ for two different substances. As pointed out earlier, this non-zero term results from the observation that for an incompressible fluid, only the true densities $\rho _k^T$ are constant. However, if we assumed incompressibility simply as the constancy of the overall densities $\rho _k$ , the right-hand side of (2.16c ) would be zero, derived directly from the diffusion equation (2.8). Thus, our approach can be viewed as a generalisation of a special case with $\tilde {{\boldsymbol{\nabla }}}\boldsymbol{\cdot }\tilde {\boldsymbol{u}}=0$ , corresponding to $\rho _u^T=\rho _v^T$ , used for example by Dalwadi et al. (Reference Dalwadi, Griffiths and Bruna2015). We aim to explore in depth the effects of this non-solenoidal mixture velocity field. The final missing element is the boundary condition at the stationary, rigid, impermeable solid skeleton. We impose both no slip and no penetration of the barycentric advective velocity,

(2.16d) \begin{align} \tilde {\boldsymbol{u}}=0, \hspace {10pt} \tilde {\boldsymbol{x}}\in \varGamma _s. \end{align}

This does not suppress adsorption: solute mass crosses the interface through the diffusive flux in (2.14b ), while the mixture has no advective transport through the stationary solid.

2.2.2. Microscopic governing equations

The final forms of the detailed governing equations are obtained after non-dimensionalisation, see Appendix A for more details. Here, $\delta$ is the ratio of the pore spacing to the macroscopic filter length, introduced at the beginning of § 2; throughout, $0\lt \delta \ll 1$ . The equations describing diffusion of the solute are

(2.17a) \begin{align} \frac {\partial {c}}{\partial {t}}={\boldsymbol{\nabla }} \,{\boldsymbol{\cdot }}\,\left (\frac {1}{{Pe}}{\boldsymbol{\nabla }} c-\boldsymbol{u} c\right )&, \hspace {10pt} \boldsymbol{x}\in \varOmega _{\kern-1pt f}, \\[-12pt] \nonumber \end{align}

(2.17b) \begin{align} -k\delta c=\boldsymbol{n}\boldsymbol{\cdot }\frac {1}{{Pe}}{\boldsymbol{\nabla }} c&, \hspace {10pt} \boldsymbol{x}\in \varGamma _s, \\[10pt] \nonumber \end{align}

and flow of the solvent

(2.18a) \begin{align} -{\boldsymbol{\nabla }} p+\delta ^2{{\nabla }}^2\boldsymbol{u}=0&, \hspace {10pt} \boldsymbol{x}\in \varOmega _{\kern-1pt f}, \\[-12pt] \nonumber \end{align}

(2.18b) \begin{align} {\boldsymbol{\nabla }} \,{\boldsymbol{\cdot }}\,\boldsymbol{u}=\frac {A}{{Pe}}{{\nabla }}^2 c&, \hspace {10pt} \boldsymbol{x}\in \varOmega _{\kern-1pt f}, \\[-12pt] \nonumber \end{align}

(2.18c) \begin{align} \boldsymbol{u}=0&, \hspace {10pt} \boldsymbol{x}\in \varGamma _s. \\[10pt] \nonumber \end{align}

We introduced a dimensionless parameter, $A= [(1/\rho _u^T)- (1/\rho _v^T )]c_\infty$ , with $c_\infty$ being a characteristic solute concentration outside of the filter. Note that the ratio $A/Pe$ measures the deviation from the classical incompressibility condition for a single component fluid, which would lead to a zero right-hand side of (A11b ), as described in § 2.2.1.

3.1.1. Transformation of derivatives

As is standard, we replace the spatial gradients as

(3.5) \begin{align} \boldsymbol{\nabla }\longmapsto {\boldsymbol{\nabla }}_{x}+\frac {1}{\delta }{\boldsymbol{\nabla }}_{y}. \end{align}

Hence, for the fluid system (2.18), we obtain after a multiplication by a suitable power of $\delta$ (so that the lowest power of $\delta$ is equal to $\delta ^{0}$ ),

(3.6a) \begin{align} -\left (\delta {\boldsymbol{\nabla }}_{x}+{\boldsymbol{\nabla }}_{y}\right ) p+\delta \left (\delta {\boldsymbol{\nabla }}_{x}+{\boldsymbol{\nabla }}_{y}\right )^2\boldsymbol{u}=0&, \hspace {10pt} \boldsymbol{y}\in \omega _{\kern-1pt f}(\boldsymbol{x}), \\[-12pt] \nonumber \end{align}

(3.6b) \begin{align} \big (\delta ^2{\boldsymbol{\nabla }}_{x}+\delta {\boldsymbol{\nabla }}_{y}\big )\boldsymbol{\cdot }\boldsymbol{u}=\frac {A}{{Pe}}\left (\delta {\boldsymbol{\nabla }}_{x}+{\boldsymbol{\nabla }}_{y}\right )^2 c&, \hspace {10pt} \boldsymbol{y}\in \omega _{\kern-1pt f}(\boldsymbol{x}), \\[-12pt] \nonumber \end{align}

(3.6c) \begin{align} \boldsymbol{u}=0&, \hspace {10pt} \boldsymbol{y}\in \partial \omega _s(\boldsymbol{x}). \\[10pt] \nonumber \end{align}

For the diffusion system (2.17), we proceed analogously, however, the inhomogeneous boundary condition requires further care. We need to first describe the normal vector $\boldsymbol{n}$ , perpendicular to the interface, and depending on both $\boldsymbol{x}$ and $\boldsymbol{y}$ . We now introduce a level set function describing the location of the interface, $\chi (\boldsymbol{x},\boldsymbol{y})=R(\boldsymbol{x})-\|\boldsymbol{y}\|$ (hence, a geometry more general than a ball is admissible), while the interface is simply $\chi (\boldsymbol{x},\boldsymbol{y})=0$ . Now, we just need to realise that the normal vector to the boundary points in the direction as follows (Klika & Gaffney Reference Klika and Gaffney2023):

(3.7) \begin{align} \boldsymbol{n}\propto \boldsymbol{\nabla }\chi =\left ({\boldsymbol{\nabla }}_{x}+\frac {1}{\delta }{\boldsymbol{\nabla }}_{y}\right )\chi (\boldsymbol{x},\boldsymbol{y})={\boldsymbol{\nabla }}_{x} R(\boldsymbol{x})-\frac {1}{\delta }\frac {\boldsymbol{y}}{\|\boldsymbol{y}\|}={\boldsymbol{\nabla }}_{x} R(\boldsymbol{x})+\frac {1}{\delta }\boldsymbol{n}_y, \end{align}

where we defined a new unit vector on the microscale $\boldsymbol{n}_y=-({\boldsymbol{y}}/{\|\boldsymbol{y}\|})$ . The normal vector $\boldsymbol{n}$ can be then in line with previous reasoning redefined as

(3.8) \begin{align} \boldsymbol{n}(\boldsymbol{x},\boldsymbol{y})=\frac {\boldsymbol{n}_y+\delta {\boldsymbol{\nabla }}_{x} R(\boldsymbol{x})}{\|\boldsymbol{n}_y+\delta {\boldsymbol{\nabla }}_{x} R(\boldsymbol{x})\|}. \end{align}

This normal follows from the total gradient of the particular physical level set $\chi (\boldsymbol{x},\boldsymbol{x}/\delta )$ . It is therefore not an additional assumption that $\boldsymbol n$ decouples into unrelated slow and fast fields. Alternative parametrised-cell approaches retain a general $\boldsymbol n(\boldsymbol{x},\boldsymbol{y})$ and generate extra geometric terms when transport identities are applied (Burridge & Keller Reference Burridge and Keller1981; Holmes Reference Holmes2013; Penta et al. Reference Penta, Ambrosi and Shipley2014, Reference Penta, Ambrosi and Quarteroni2015). Those formulations address a broader geometry, while here we restrict to the level-set construction of Bruna & Chapman (Reference Bruna and Chapman2015) and Klika & Gaffney (Reference Klika and Gaffney2023). A rigorous comparison of the two formal limits is beyond the scope of the present fluid-mechanical model.

We may now transform the diffusion problem into the multiple-scales form, which reads (after a multiplication by a suitable power of $\delta$ )

(3.9a) \begin{align} \delta ^2\frac {\partial {c}}{\partial {t}}=\left (\delta {\boldsymbol{\nabla }}_{x}+{\boldsymbol{\nabla }}_{y}\right )\boldsymbol{\cdot }\left (\frac {1}{{Pe}}\left (\delta {\boldsymbol{\nabla }}_{x}+{\boldsymbol{\nabla }}_{y}\right ) c-\delta \boldsymbol{u} c\right )&, \hspace {10pt} \boldsymbol{y}\in \omega _{\kern-1pt f}(\boldsymbol{x}), \\[-12pt] \nonumber \end{align}

(3.9b) \begin{align} -k\delta ^2 c=\left (\boldsymbol{n}_y+\delta {\boldsymbol{\nabla }}_{x} R(\boldsymbol{x})\right )\boldsymbol{\cdot }\frac {1}{{Pe}}\left (\delta {\boldsymbol{\nabla }}_{x}+{\boldsymbol{\nabla }}_{y}\right ) c+\textit {O} (\delta ^3)&, \hspace {10pt} \boldsymbol{y}\in \partial \omega _s(\boldsymbol{x}), \\[10pt] \nonumber \end{align}

where we also multiplied (3.9b ) by the term $\|\boldsymbol{n}_y+\delta {\boldsymbol{\nabla }}_{x} R(\boldsymbol{x})\|$ .

3.1.2. Multiple-scale expansion

We now expand the state variables in $\delta$ into regular asymptotic expansions:

(3.10) \begin{align} c(\boldsymbol{x},\boldsymbol{y},t)&=c_0(\boldsymbol{x},\boldsymbol{y},t)+\delta c_1(\boldsymbol{x},\boldsymbol{y},t)+\delta ^2 c_2(\boldsymbol{x},\boldsymbol{y},t)+{\cdots} , \\[-12pt] \nonumber \end{align}

(3.11) \begin{align} \boldsymbol{u}(\boldsymbol{x},\boldsymbol{y},t)&=\boldsymbol{u}_0(\boldsymbol{x},\boldsymbol{y},t)+\delta \boldsymbol{u}_1(\boldsymbol{x},\boldsymbol{y},t)+\delta ^2 \boldsymbol{u}_2(\boldsymbol{x},\boldsymbol{y},t)+{\cdots} , \\[-12pt] \nonumber \end{align}

(3.12) \begin{align} p(\boldsymbol{x},\boldsymbol{y},t)&=p_0(\boldsymbol{x},\boldsymbol{y},t)+\delta p_1(\boldsymbol{x},\boldsymbol{y},t)+\delta ^2 p_2(\boldsymbol{x},\boldsymbol{y},t)+{\cdots} , \\[10pt] \nonumber \end{align}

which finally gives the following form of the governing equations (3.6) and (3.9):

(3.13a) \begin{align} -\left (\delta {\boldsymbol{\nabla }}_{x}+{\boldsymbol{\nabla }}_{y}\right )\left (p_0+\delta p_1+\delta ^2 p_2\right )+\delta \left (\delta {\boldsymbol{\nabla }}_{x}+{\boldsymbol{\nabla }}_{y}\right )^2\left (\boldsymbol{u}_0+\delta \boldsymbol{u}_1\right )=0&, \hspace {10pt} \boldsymbol{y}\in \omega _{\kern-1pt f}(\boldsymbol{x}), \\[-12pt] \nonumber \end{align}

(3.13b) \begin{align} \big (\delta ^2{\boldsymbol{\nabla }}_{x}+\delta {\boldsymbol{\nabla }}_{y}\big )\boldsymbol{\cdot }\left (\boldsymbol{u}_0+\delta \boldsymbol{u}_1\right )=\frac {A}{{Pe}}\left (\delta {\boldsymbol{\nabla }}_{x}+{\boldsymbol{\nabla }}_{y}\right )^2 \big (c_0+\delta c_1+\delta ^2 c_2\big )&, \hspace {10pt} \boldsymbol{y}\in \omega _{\kern-1pt f}(\boldsymbol{x}), \\[-12pt] \nonumber \end{align}

(3.13c) \begin{align} \boldsymbol{u}_0+\delta \boldsymbol{u}_1+\delta ^2 \boldsymbol{u}_2=0&, \hspace {10pt} \boldsymbol{y}\in \partial \omega _s(\boldsymbol{x}) \\[10pt] \nonumber \end{align}

and

(3.14a) \begin{align} \delta ^2\frac {\partial {c_0}}{\partial {t}} & =\left (\delta {\boldsymbol{\nabla }}_{x}+{\boldsymbol{\nabla }}_{y}\right )\boldsymbol{\cdot }\left (\frac {1}{{Pe}}\left (\delta {\boldsymbol{\nabla }}_{x}+{\boldsymbol{\nabla }}_{y}\right )\big (c_0+\delta c_1+\delta ^2 c_2\big ) \right. \nonumber \\& \quad \left. -\delta \left (\boldsymbol{u}_0+\delta \boldsymbol{u}_1\right ) \vphantom{\frac{1}{2}}\left (c_0+\delta c_1\right )\right )\!, \hspace {10pt} \boldsymbol{y}\in \omega _{\kern-1pt f}(\boldsymbol{x}), \\[-10pt] \nonumber \end{align}

(3.14b) \begin{align} -k\delta ^2 c_0 & =\left (\boldsymbol{n}_y+\delta {\boldsymbol{\nabla }}_{x} R(\boldsymbol{x})\right )\boldsymbol{\cdot }\frac {1}{{Pe}}\left (\delta {\boldsymbol{\nabla }}_{x}+{\boldsymbol{\nabla }}_{y}\right )\big (c_0+\delta c_1+\delta ^2 c_2\big ), \hspace {10pt} \boldsymbol{y}\in \partial \omega _s(\boldsymbol{x}), \\[10pt] \nonumber \end{align}

where we kept only terms up to the power of $\delta ^2$ . Now, we put together all terms with the same power of $\delta$ , which gets us equations for the leading order solutions $c_0,\boldsymbol{u}_0,p_0$ and corrections $c_i,\boldsymbol{u}_i,p_i, \hspace {3pt} i\in \mathbb{N}$ .

3.1.3. Order-by-order solutions

First, let us assume that $A/{Pe}=\textit {O}(\delta ^{0})$ . Further, note that the no-slip fluid flow condition translates for all $i\in \mathbb{N}_0$ into the same condition

(3.15) \begin{align} \boldsymbol{u}_i=0, \hspace {10pt} \boldsymbol{y}\in \partial \omega _s. \end{align}

Collecting the leading order problem $\textit {O}(\delta ^{0})$ in the equations in (3.13) and (3.14), we obtain

(3.16a) \begin{align} {\boldsymbol{\nabla }}_{y}^2 c_0=0&, \hspace {10pt} \boldsymbol{y}\in \omega _{\kern-1pt f}(\boldsymbol{x}), \\[-12pt] \nonumber \end{align}

(3.16b) \begin{align} \boldsymbol{n}_y\boldsymbol{\cdot }{\boldsymbol{\nabla }}_{y} c_0=0&, \hspace {10pt} \boldsymbol{y}\in \partial \omega _s(\boldsymbol{x}). \\[10pt] \nonumber \end{align}

The solution of (3.16) is $c_0=c_0(\boldsymbol{x},t)$ , meaning that $c_0$ is independent of $\boldsymbol{y}$ . The homogeneous Neumann problem determines only that the fast-variable dependence is constant; that constant is the macroscopic field $c_0(\boldsymbol{x},t)$ .

In the same manner, we obtain the leading-order equations for (3.13) as

(3.17a) \begin{align} {\boldsymbol{\nabla }}_{y} p_0=0&, \hspace {10pt} \boldsymbol{y}\in \omega _{\kern-1pt f}(\boldsymbol{x}), \\[-12pt] \nonumber \end{align}

(3.17b) \begin{align} {\boldsymbol{\nabla }}_{y}^2 c_0=0&, \hspace {10pt} \boldsymbol{y}\in \omega _{\kern-1pt f}(\boldsymbol{x}). \\[10pt] \nonumber \end{align}

Equation (3.17b ) is already satisfied by the solution $c_0=c_0(\boldsymbol{x},t)$ of (3.16) and (3.17a ) is fulfilled by $p_0=p_0(\boldsymbol{x},t)$ being independent of $\boldsymbol{y}$ . Both $c$ and $p$ are therefore independent of the microscale $\boldsymbol{y}$ in the leading order.

We now proceed to the first-order subleading problem, that is, $\textit {O}(\delta )$ . Collecting terms with $\delta ^1$ in (3.13), and utilising $c_0=c_0(\boldsymbol{x},t)$ , yields

(3.18a) \begin{align} -{\boldsymbol{\nabla }}_{y} p_1+{{\nabla }}_{y}^2\boldsymbol{u}_0={\boldsymbol{\nabla }}_{x} p_0&, \hspace {10pt} \boldsymbol{y}\in \omega _{\kern-1pt f}(\boldsymbol{x}), \\[-12pt] \nonumber \end{align}

(3.18b) \begin{align} {\boldsymbol{\nabla }}_{y}\boldsymbol{\cdot }\boldsymbol{u}_0=\frac {A}{{Pe}}{{\nabla }}_{y}^2 c_1&, \hspace {10pt} \boldsymbol{y}\in \omega _{\kern-1pt f}(\boldsymbol{x}). \\[10pt] \nonumber \end{align}

Similarly, the first subleading problem for solute transport reads

(3.19a) \begin{align} 0=\frac {1}{{Pe}}{\boldsymbol{\nabla }}_{y}^2 c_1-c_0{\boldsymbol{\nabla }}_{y}\boldsymbol{\cdot }\boldsymbol{u}_0&, \hspace {10pt} \boldsymbol{y}\in \omega _{\kern-1pt f}(\boldsymbol{x}), \\[-12pt] \nonumber \end{align}

(3.19b) \begin{align} 0=\boldsymbol{n}_y\boldsymbol{\cdot }\left (\frac {1}{{Pe}}{\boldsymbol{\nabla }}_{x} c_0 +\frac {1}{{Pe}}{\boldsymbol{\nabla }}_{y} c_1\right )&, \hspace {10pt} \boldsymbol{y}\in \partial \omega _s(\boldsymbol{x}), \\[10pt] \nonumber \end{align}

where we used the boundary condition (3.15) and the independence of $c_{0}$ on $y$ .

We substitute ${\boldsymbol{\nabla }}_{y}\boldsymbol{\cdot }\boldsymbol{u}_0= (A/{Pe}){\boldsymbol{\nabla }}_{y}^2 c_1$ from (3.18b ) into (3.19a ) and get

(3.20) \begin{align} \left (1-A c_0\right ){\boldsymbol{\nabla }}_{y}^2 c_1=0&, \hspace {10pt} \boldsymbol{y}\in \omega _{\kern-1pt f}(\boldsymbol{x}), \nonumber \\ \boldsymbol{n}_y\boldsymbol{\cdot }{\boldsymbol{\nabla }}_{y} c_1=-\boldsymbol{n}_y\boldsymbol{\cdot }{\boldsymbol{\nabla }}_{x} c_0&, \hspace {10pt} \boldsymbol{y}\in \partial \omega _s(\boldsymbol{x}). \end{align}

Note that the factor $ (1-Ac_0 )$ is non-zero. Recalling definitions $A= [(1/\rho _2^T)- (1/\rho _1^T )]c_\infty$ and $\tilde {c}_0=c_\infty c_0$ entails

(3.21) \begin{align} 1-Ac_0=1-\left (\frac {1}{\rho _2^T}-\frac {1}{\rho _1^T}\right )c_\infty c_0=1-\frac {\tilde {c}_0}{\rho _2^T}+\frac {\tilde {c}_0}{\rho _1^T}. \end{align}

If $\tilde {c}_0=0$ , then also $c_0=0$ and we simply obtain $1-Ac_0=1$ . Now, let us assume that $\tilde {c}_0\gt 0$ . From the definition of the true density, we have that $\rho _1^T\gt 0$ and $\rho _2^T\geq \tilde {c}_0$ , and thus, $1- (\tilde {c}_0/\rho _2^T)\geq 0$ . All together, then

(3.22) \begin{align} 1-Ac_0=1-\frac {\tilde {c}_0}{\rho _2^T}+\frac {\tilde {c}_0}{\rho _1^T}\gt 0. \end{align}

Hence, we may simplify the first subleading diffusion problem into

(3.23a) \begin{align} {\boldsymbol{\nabla }}_{y}^2 c_1=0&, \hspace {10pt} \boldsymbol{y}\in \omega _{\kern-1pt f}(\boldsymbol{x}), \\[-12pt] \nonumber \end{align}

(3.23b) \begin{align} \boldsymbol{n}_y\boldsymbol{\cdot }{\boldsymbol{\nabla }}_{y} c_1=-\boldsymbol{n}_y\boldsymbol{\cdot }{\boldsymbol{\nabla }}_{x} c_0&, \hspace {10pt} \boldsymbol{y}\in \partial \omega _s(\boldsymbol{x}). \\[10pt] \nonumber \end{align}

Equations (3.23) are linear in $c_1$ and uncoupled from fluid flow equations (3.18). Their general solution contains a $\boldsymbol{y}$ -independent first-order mean,

(3.24) \begin{align} c_1(\boldsymbol{x},\boldsymbol{y},t)=-{\boldsymbol{\nabla }}_{x} c_0(\boldsymbol{x},t)\boldsymbol{\cdot }{\varGamma }(\boldsymbol{x},\boldsymbol{y})+\bar c_1(\boldsymbol{x},t), \end{align}

where ${\varGamma }(\boldsymbol{x},\boldsymbol{y})$ is a vector satisfying

(3.25a) \begin{align} {\boldsymbol{\nabla }}_{y}^2\varGamma _i=0&, \hspace {10pt} \boldsymbol{y}\in \omega _{\kern-1pt f}(\boldsymbol{x}), \\[-12pt] \nonumber \end{align}

(3.25b) \begin{align} \boldsymbol{n}_y\boldsymbol{\cdot }{\boldsymbol{\nabla }}_{y}\varGamma _i=n_{y,i}&, \hspace {10pt} \boldsymbol{y}\in \partial \omega _s(\boldsymbol{x}), \\[-12pt] \nonumber \end{align}

(3.25c) \begin{align} \varGamma _i \hspace {5pt} \text{periodic}&, \hspace {10pt} \boldsymbol{y}\in \partial \omega (\boldsymbol{x}), \\[-12pt] \nonumber \end{align}

(3.25d) \begin{align} \int _{\omega _{\kern-1pt f}(\boldsymbol{x})}\varGamma _i\,\,\mathrm{d}\boldsymbol{y}=0.& \\[10pt] \nonumber \end{align}

Here, notation $n_{y,i}$ denotes the ith component of the vector $\boldsymbol{n}_y$ .

The last condition fixes the additive constant in the Neumann cell problem. Correspondingly, we choose the standard corrector gauge $\int _{\omega _{\kern-1pt f}(\boldsymbol{x})}c_1\,\,\mathrm{d}\boldsymbol{y}=0$ , which sets $\bar c_1=0$ . This gauge leaves the leading-order macroscopic equations unchanged; retaining $\bar c_1$ would only alter the $O(\delta )$ reconstruction. Hence, the formula used in the following is (3.24) with $\bar c_1=0$ .

The above-mentioned set of equations is called the cell problem and its solution ${\varGamma }$ is given directly by the geometry of our model. Also, let us note that in (3.25), the coordinate $\boldsymbol{x}$ appears only through the choice of the unit cell $\omega (\boldsymbol{x})$ . This means we can find the function ${\varGamma }$ for all cells simultaneously, where $x$ acts as a parameter. For further theoretical evaluation, we consider the function ${\varGamma }$ as already known.

Now, returning to the flow problem, the right-hand side of (3.18b ) vanishes as ${\boldsymbol{\nabla }}_{y}^2 c_1=0$ . Thus, together with the boundary condition (3.15), we obtain

(3.26a) \begin{align} -{\boldsymbol{\nabla }}_{y} p_1+{\boldsymbol{\nabla }}_{y}^2\boldsymbol{u}_0={\boldsymbol{\nabla }}_{x} p_0&, \hspace {10pt} \boldsymbol{y}\in \omega _{\kern-1pt f}(\boldsymbol{x}), \\[-12pt] \nonumber \end{align}

(3.26b) \begin{align} {\boldsymbol{\nabla }}_{y}\boldsymbol{\cdot }\boldsymbol{u}_0=0&, \hspace {10pt} \boldsymbol{y}\in \omega _{\kern-1pt f}(\boldsymbol{x}), \\[-12pt] \nonumber \end{align}

(3.26c) \begin{align} \boldsymbol{u}_0=0&, \hspace {10pt} \boldsymbol{y}\in \partial \omega _s(\boldsymbol{x}), \\[-12pt] \nonumber \end{align}

(3.26d) \begin{align} \boldsymbol{u}_0,p_1\ \text{periodic}&, \hspace {10pt} \boldsymbol{y}\in \partial \omega (\boldsymbol{x}). \\[10pt] \nonumber \end{align}

Equations (3.26) are again linear in $p_1$ and $\boldsymbol{u}_0$ , therefore, we can also factorise the solution as

(3.27a) \begin{align} \boldsymbol{u}_0(\boldsymbol{x},\boldsymbol{y},t)&=-\unicode{x1D646}(\boldsymbol{x},\boldsymbol{y},t)\boldsymbol{\cdot }{\boldsymbol{\nabla }}_{x} p_0(\boldsymbol{x},t), \\[-12pt] \nonumber \end{align}

(3.27b) \begin{align} p_1(\boldsymbol{x},\boldsymbol{y},t)&=-{\varPi }(\boldsymbol{x},\boldsymbol{y},t)\boldsymbol{\cdot }{\boldsymbol{\nabla }}_{x} p_0(\boldsymbol{x},t)+\bar {p}(\boldsymbol{x},t), \\[10pt] \nonumber \end{align}

where $\unicode{x1D646}(\boldsymbol{x},\boldsymbol{y},t)$ is a matrix and ${\varPi }(\boldsymbol{x},\boldsymbol{y},t)$ a vector satisfying

(3.28a) \begin{align} \unicode{x1D644}_d-{\boldsymbol{\nabla }}_{y}{\varPi }+{\boldsymbol{\nabla }}_{y}^2\unicode{x1D646}=0&, \hspace {10pt} \boldsymbol{y}\in \omega _{\kern-1pt f}(\boldsymbol{x}), \\[-12pt] \nonumber \end{align}

(3.28b) \begin{align} {\boldsymbol{\nabla }}_{y}\boldsymbol{\cdot }\unicode{x1D646}=0&, \hspace {10pt} \boldsymbol{y}\in \omega _{\kern-1pt f}(\boldsymbol{x}), \\[-12pt] \nonumber \end{align}

(3.28c) \begin{align} \unicode{x1D646}=0&, \hspace {10pt} \boldsymbol{y}\in \partial \omega _s(\boldsymbol{x}), \\[-12pt] \nonumber \end{align}

(3.28d) \begin{align} \unicode{x1D646},{\varPi } \hspace {5pt} \text{periodic}&, \hspace {10pt} \boldsymbol{y}\in \partial \omega (\boldsymbol{x}), \\[-12pt] \nonumber \end{align}

(3.28e) \begin{align} \int _{\omega _{\kern-1pt f}(\boldsymbol{x})}{\varPi }\,\,\mathrm{d}\boldsymbol{y}=\boldsymbol 0.& \\[10pt] \nonumber \end{align}

The final condition fixes the arbitrary additive constant in the pressure corrector ${\varPi }$ ; the velocity corrector $\unicode{x1D646}$ is already fixed by its no-slip boundary condition. Here, $\unicode{x1D644}_d$ denotes the identity matrix for the corresponding spatial dimension $d$ , and the terms ${\boldsymbol{\nabla }}_{y}{\varPi },{\boldsymbol{\nabla }}_{y}^2\unicode{x1D646}$ represent matrices and ${\boldsymbol{\nabla }}_{y}\boldsymbol{\cdot }\unicode{x1D646}$ a vector defined for $i,j\in \{1,{\cdots} \hspace {-0.7pt},d\}$ as follows:

(3.29) \begin{align} \left ({\boldsymbol{\nabla }}_{y}{\varPi }\right )_{\bullet j}&={\boldsymbol{\nabla }}_{y}{\varPi }_j, \nonumber \\ \left ({\boldsymbol{\nabla }}_{y}^2\unicode{x1D646}\right )_{ij}&={\boldsymbol{\nabla }}_{y}^2\unicode{x1D646}_{ij}, \nonumber \\ \left ({\boldsymbol{\nabla }}_{y}\boldsymbol{\cdot }\unicode{x1D646}\right )_i&={\boldsymbol{\nabla }}_{y}\boldsymbol{\cdot }\unicode{x1D646}_{\bullet i}. \end{align}

As in the diffusion problem, by the factorisation, we managed to split the dependence on $\boldsymbol{x}$ and $\boldsymbol{y}$ . Functions $\unicode{x1D646}$ and ${\varPi }$ are again dependent on $\boldsymbol{x}$ only through the change of the boundary $\partial \omega _s(\boldsymbol{x})$ and we can evaluate them for all cells $\omega (\boldsymbol{x})$ at once. As with the function ${\varGamma }$ , we, for further progress, assume that the functions $\unicode{x1D646}$ and ${\varPi }$ are known solutions to (3.28).

3.2.1. Diffusion

For the diffusion system (3.14), we continue with the asymptotic expansion to the second order in $\delta$ :

(3.30a) \begin{align} \frac {\partial {c_0}}{\partial {t}} &={\boldsymbol{\nabla }}_{y}\boldsymbol{\cdot }\left (\frac {1}{{Pe}}\left ({\boldsymbol{\nabla }}_{x} c_1+{\boldsymbol{\nabla }}_{y} c_2\right )-\boldsymbol{u}_0 c_1-\boldsymbol{u}_1 c_0\right )\qquad \nonumber \\& \quad +{\boldsymbol{\nabla }}_{x}\boldsymbol{\cdot }\left (\frac {1}{{Pe}}\left ({\boldsymbol{\nabla }}_{x} c_0+{\boldsymbol{\nabla }}_{y} c_1\right )-\boldsymbol{u}_0 c_0\right ), \hspace {10pt} \boldsymbol{y}\in \omega _{\kern-1pt f}(\boldsymbol{x}), \\[-12pt] \nonumber \end{align}

(3.30b) \begin{align} -kc_0 & =\boldsymbol{n}_y\boldsymbol{\cdot }\frac {1}{{Pe}}\left ({\boldsymbol{\nabla }}_{x} c_1+{\boldsymbol{\nabla }}_{y} c_2\right )+\qquad \nonumber \\& \quad +{\boldsymbol{\nabla }}_{x} R(\boldsymbol{x})\boldsymbol{\cdot }\frac {1}{{Pe}}\left ({\boldsymbol{\nabla }}_{x} c_0+{\boldsymbol{\nabla }}_{y} c_1\right ), \hspace {10pt} \boldsymbol{y}\in \partial \omega _s(\boldsymbol{x}), \\[10pt] \nonumber \end{align}

Next, we integrate (3.30a ) over the fluid volume $\omega _{\kern-1pt f}(\boldsymbol{x})$ as

(3.31) \begin{align} \int _{\omega _{\kern-1pt f}(\boldsymbol{x})}\frac {\partial {c_0}}{\partial {t}}\,\mathrm{d}\boldsymbol{y} &= \nonumber \\& ={\boldsymbol{\nabla }}_{x}\boldsymbol{\cdot }\int _{\omega _{\kern-1pt f}(\boldsymbol{x})}\left (\frac {1}{{Pe}}\left ({\boldsymbol{\nabla }}_{x} c_0+{\boldsymbol{\nabla }}_{y} c_1\right )-\boldsymbol{u}_0 c_0\right )\,\mathrm{d}\boldsymbol{y}-\int _{\partial \omega _s(\boldsymbol{x})}\hspace {-12pt}kc_0\,{\mathrm ds}, \end{align}

where we used the divergence theorem, the boundary condition (3.30b ) and the Reynolds transport theorem (Dalwadi et al. Reference Dalwadi, Griffiths and Bruna2015; Klika & Gaffney Reference Klika and Gaffney2023)

(3.32) \begin{align} {\boldsymbol{\nabla }}_{x}\boldsymbol{\cdot }\int _{\omega _{\kern-1pt f}(\boldsymbol{x})}\boldsymbol{v}\,\mathrm{d}\boldsymbol{y}=\int _{\omega _{\kern-1pt f}(\boldsymbol{x})}{\boldsymbol{\nabla }}_{x}\boldsymbol{\cdot }\boldsymbol{v}\,\mathrm{d}\boldsymbol{y}-\int _{\partial \omega _s(\boldsymbol{x})}{\boldsymbol{\nabla }}_{x} R(\boldsymbol{x})\boldsymbol{\cdot }\boldsymbol{v}\,{\mathrm ds}. \end{align}

When the $\boldsymbol{y}$ -divergence theorem is applied on $\omega _{\kern-1pt f}(\boldsymbol{x})$ , its boundary consists of the solid interface and the outer faces of the reference cell. The integrals on each pair of opposite outer faces cancel because the correctors and their fluxes have matching periodic boundary values in $\boldsymbol{y}$ . The remaining solid-interface integral is combined with the boundary condition and with the transport identity (3.32). The latter is associated with the level-set cell family chosen previously.

Note that (3.31) give the same result as the Fredholm alternative theorem (Allaire et al. Reference Allaire, Brizzi, Mikelić and Piatnitski2010), which states that the system (3.30) for the unknown $c_2$ has a solution only if a condition given by (3.31) is satisfied.

Now, given $c_0(\boldsymbol{x},t)$ is independent of $\boldsymbol{y}$ and that $\int _{\omega _{\kern-1pt f}(\boldsymbol{x})}\,\mathrm{d}\boldsymbol{y}=|\omega _{\kern-1pt f}(\boldsymbol{x})|=\phi (\boldsymbol{x})$ , (3.31) then simplifies to

(3.33) \begin{align} \phi (\boldsymbol{x})\frac {\partial {c_0}}{\partial {t}}={\boldsymbol{\nabla }}_{x}\boldsymbol{\cdot }\left (\frac {1}{{Pe}}\int _{\omega _{\kern-1pt f}(\boldsymbol{x})}\left ({\boldsymbol{\nabla }}_{x} c_0+{\boldsymbol{\nabla }}_{y} c_1\right )\,\mathrm{d}\boldsymbol{y}-c_0\int _{\omega _{\kern-1pt f}(\boldsymbol{x})}\boldsymbol{u}_0\,\mathrm{d}\boldsymbol{y}\right )-kc_0|\partial \omega _s(\boldsymbol{x})|. \end{align}

If we define the Jacobi matrix $ (\unicode{x1D645}_\varGamma )_{ij}= \partial {\varGamma _i}/\partial {y_j}, \hspace {3pt} i,j\in \{1,{\cdots} \hspace {-0.7pt},d\}$ , we can, using (3.24), write

(3.34) \begin{align} {\boldsymbol{\nabla }}_{y} c_1={\boldsymbol{\nabla }}_{y}\left (-{\boldsymbol{\nabla }}_{x} c_0\boldsymbol{\cdot }{\varGamma }\right )={\boldsymbol{\nabla }}_{y}\left (-\sum _{i}\varGamma _i\frac {\partial {c_0}}{\partial {x_i}}\right )=-\sum _{i}\left ({\boldsymbol{\nabla }}_{y}\varGamma _i\right )\frac {\partial {c_0}}{\partial {x_i}}=-\unicode{x1D645}_\varGamma ^T\boldsymbol{\cdot }{\boldsymbol{\nabla }}_{x} c_0. \end{align}

Therefore, we have another form of (3.33) given only in terms of $c_{0}$ :

(3.35) \begin{align} \phi (\boldsymbol{x})\frac {\partial {c_0}}{\partial {t}}&={\boldsymbol{\nabla }}_{x}\boldsymbol{\cdot }\left (\frac {1}{{Pe}}\int _{\omega _{\kern-1pt f}(\boldsymbol{x})}\left (\unicode{x1D644}_d-\unicode{x1D645}_\varGamma ^T\right )\boldsymbol{\cdot }{\boldsymbol{\nabla }}_{x} c_0\,\mathrm{d}\boldsymbol{y}-c_0\int _{\omega _{\kern-1pt f}(\boldsymbol{x})}\boldsymbol{u}_0\,\mathrm{d}\boldsymbol{y}\right )-kc_0|\partial \omega _s(\boldsymbol{x})| \nonumber \\[5pt]&={\boldsymbol{\nabla }}_{x}\boldsymbol{\cdot }\left (\frac {1}{{Pe}}\left (\phi (\boldsymbol{x})\unicode{x1D644}_d-\int _{\omega _{\kern-1pt f}(\boldsymbol{x})}\unicode{x1D645}_\varGamma ^T\,\mathrm{d}\boldsymbol{y}\right )\boldsymbol{\cdot }{\boldsymbol{\nabla }}_{x} c_0-c_0\int _{\omega _{\kern-1pt f}(\boldsymbol{x})}\boldsymbol{u}_0\,\mathrm{d}\boldsymbol{y}\right ) \nonumber \\[4pt]& \quad -kc_0|\partial \omega _s(\boldsymbol{x})|. \end{align}

Noting that $C(\boldsymbol{x},t)\sim \phi (\boldsymbol{x})c_0(\boldsymbol{x},t)$ at the leading order, we can rewrite (3.35) into

(3.36) \begin{align} \frac {\partial {C}}{\partial {t}}={\boldsymbol{\nabla }}_{x}\boldsymbol{\cdot }\left (\tilde {D}\boldsymbol{\cdot }{\boldsymbol{\nabla }}_{x} C-\frac {C}{\phi }\left (\boldsymbol{U}+\tilde {D}\boldsymbol{\cdot }{\boldsymbol{\nabla }}_{x}\phi \right )\right )-fC, \end{align}

where $\tilde {D}$ is the effective diffusion coefficient and $f$ the effective adsorption coefficient defined as

(3.37) \begin{gather} \tilde {D}(\phi )=\frac {1}{{Pe}}\left (\unicode{x1D644}_d-\frac {1}{\phi }\int _{\omega _{\kern-1pt f}(\boldsymbol{x})}\unicode{x1D645}_\varGamma ^T\,\mathrm{d}\boldsymbol{y}\right )\!, \\[-28pt] \nonumber \end{gather}

(3.38) \begin{gather} f(\phi )=k\frac {|\partial \omega _s(\boldsymbol{x})|}{\phi }. \\[0pt] \nonumber \end{gather}

Figure 2.

Evaluation of the effective macroscopic transport coefficients as a function of porosity for a two-dimensional square lattice. (a) Diffusion coefficient. The numerical evaluation (dashed line) is obtained by solving the cell problem (3.25) for various porosity values within the range $\phi \in [1-\pi /4,1]$ using the FlexPDE software. The results are then plugged into (3.37) to obtain the plotted values. The closed-form approximation (solid line) represents evaluated expression (4.59). (b) Numerical solution of the cell problem (3.28) yielding the effective permeability coefficient $\unicode{x1D646}=K\unicode{x1D644}_d$ as a function of porosity. The numerical solution requires care and is described in Appendix B.

As an illustration, we consider spherical inclusions as a representation of filter geometry, that is, by a collection of $d$ -dimensional balls. In such a case, we can explicitly rewrite the function $f(\phi )$ as

(3.39) \begin{align} f(\phi )=k\frac {d(1-\phi )}{\phi }\left (\frac {V_d}{1-\phi }\right )^{\tfrac {1}{d}}. \end{align}

It is a direct consequence of the observations

(3.40) \begin{align} \phi (\boldsymbol{x})=|\omega _{\kern-1pt f}(\boldsymbol{x})|=1-|\omega _s(\boldsymbol{x})|=1-V_dR(\boldsymbol{x})^d \hspace {10pt}\text{and}\hspace {10pt} |\partial \omega _s(\boldsymbol{x})|= \textrm{d}V_dR(\boldsymbol{x})^{d-1}, \end{align}

where $V_d$ represents a volume of a unit ball in $d$ dimensions.

We illustrate the values of $\tilde {D}(\phi )$ for spherical inclusions (filter geometry) in two dimensions in figure 2. Note that analytical approximations for several two-dimensional (2-D) packing of circles is available due to Rayleigh, see § 4.3 for details.

3.2.2. Fluid flow

As the last step of the homogenisation process, we find macroscopic equations for the pressure $p$ and the velocity $\boldsymbol{U}$ . We obtain a relation between $p$ and $\boldsymbol{U}$ at the leading order by integrating (3.27a ):

(3.41) \begin{align} \boldsymbol{U}(\boldsymbol{x},t)&\sim \int _{\omega _{\kern-1pt f}(\boldsymbol{x})}\boldsymbol{u}_0(\boldsymbol{x},\boldsymbol{y},t)\,\mathrm{d}\boldsymbol{y}=-\int _{\omega _{\kern-1pt f}(\boldsymbol{x})}\unicode{x1D646}(\boldsymbol{x},\boldsymbol{y},t)\boldsymbol{\cdot }{\boldsymbol{\nabla }}_{x} p_0(\boldsymbol{x},t)\,\mathrm{d}\boldsymbol{y} \nonumber \\ &\quad =-\left (\int _{\omega _{\kern-1pt f}(\boldsymbol{x})}\unicode{x1D646}(\boldsymbol{x},\boldsymbol{y},t)\,\mathrm{d}\boldsymbol{y}\right )\boldsymbol{\cdot }{\boldsymbol{\nabla }}_{x} p_0(\boldsymbol{x},t)\sim -\unicode{x1D646}_{\textit{eff}}(\boldsymbol{x})\boldsymbol{\cdot }{\boldsymbol{\nabla }}_{x} p(\boldsymbol{x},t), \end{align}

where we define the effective permeability tensor as

(3.42) \begin{align} \unicode{x1D646}_{\textit{eff}}(\boldsymbol{x})=\int _{\omega _{\kern-1pt f}(\boldsymbol{x})}\unicode{x1D646}(\boldsymbol{x},\boldsymbol{y},t)\,\mathrm{d}\boldsymbol{y}. \end{align}

For a general cell, $\unicode{x1D646}_{\textit{eff}}$ is a symmetric positive-definite tensor, as follows from the reciprocal identity proved in Appendix B. Symmetry alone does not imply isotropy. For the centred circular/spherical obstacles in the square, cubic and hexagonal lattices considered in our computations, the reflection and rotational symmetries of the cell imply $\unicode{x1D646}_{\textit{eff}}(\boldsymbol{x})=K(\phi )\unicode{x1D644}_d$ . Only in those geometries do we use the scalar notation $K(\phi )$ in the following.

The macroscopic equation for the fluid flow is obtained in a similar way as the above-mentioned diffusion equation. We have

(3.43) \begin{align} \int _{\omega _{\kern-1pt f}(\boldsymbol{x})}{\boldsymbol{\nabla }}_{x}\boldsymbol{\cdot }\boldsymbol{u}_0\,\mathrm{d}\boldsymbol{y}=\frac {A}{{Pe}}\int _{\omega _{\kern-1pt f}(\boldsymbol{x})}{{\nabla }}_{x}^2c_0\,\mathrm{d}\boldsymbol{y}, \end{align}

by the virtue of ${\boldsymbol{\nabla }}_{y} c_0=0$ and ${\boldsymbol{\nabla }}_{y}\boldsymbol{\cdot }\boldsymbol{u}_0=0$ . Next, we use the Reynolds transport theorem (3.32) and the boundary condition (3.15) to rewrite the left-hand side of (3.43),

(3.44) \begin{align} \int _{\omega _{\kern-1pt f}(\boldsymbol{x})}{\boldsymbol{\nabla }}_{x}\boldsymbol{\cdot }\boldsymbol{u}_0\,\mathrm{d}\boldsymbol{y}={\boldsymbol{\nabla }}_{x}\boldsymbol{\cdot }\int _{\omega _{\kern-1pt f}(\boldsymbol{x})}\boldsymbol{u}_0\,\mathrm{d}\boldsymbol{y}+\int _{\partial \omega _s(\boldsymbol{x})}{\boldsymbol{\nabla }}_{x} R\boldsymbol{\cdot }\boldsymbol{u}_0\,{\mathrm ds}\sim {\boldsymbol{\nabla }}_{x}\boldsymbol{\cdot }\boldsymbol{U}(\boldsymbol{x},t), \end{align}

and $C(\boldsymbol{x},t)\sim c_0(\boldsymbol{x},t)\phi (\boldsymbol{x})$ to rewrite the right-hand side,

(3.45) \begin{align} \frac {A}{{Pe}}\int _{\omega _{\kern-1pt f}(\boldsymbol{x})}{{\nabla }}_{x}^2c_0\,\mathrm{d}\boldsymbol{y}=\frac {A}{{Pe}}\phi (\boldsymbol{x}){{\nabla }}_{x}^2c_0\sim \frac {A}{{Pe}}\phi (\boldsymbol{x}){{\nabla }}_{x}^2\left (\frac {C}{\phi (\boldsymbol{x})}\right )\!. \end{align}

Finally, by combining both (3.44) and (3.45), we obtain the last macroscopic governing equation in the form

(3.46) \begin{align} {\boldsymbol{\nabla }}_{x}\boldsymbol{\cdot }\boldsymbol{U}=\frac {A}{{Pe}}\phi {{\nabla }}_{x}^2\left (\frac {C}{\phi }\right )\!, \end{align}

describing the fluid flow through the filter.

4.1.1. Upstream equations

Having the inlet and outlet conditions (4.2) and (4.3), we can extend the system (4.1) to both reservoirs and find its solution. We first proceed for the upstream variables $C^{-},U^{-}$ . We apply limit $\phi \to 1$ to the system (4.1), which yields

(4.4a) \begin{align} \frac {{d}}{\mathrm{d}x}\left (\frac {1}{{Pe}}\frac {\mathrm{d}C^{-}}{\mathrm{d}x}-C^{-}U^{-}\right )=0&, \hspace {10pt} x\in (-\infty ,0), \\[-12pt] \nonumber \end{align}

(4.4b) \begin{align} \frac {\mathrm{d}U^{-}}{\mathrm{d}x}=\frac {A}{{Pe}}\frac {\mathrm{d}^2C^{-}}{\mathrm{d} x^2}&, \hspace {10pt} x\in (-\infty ,0). \\[10pt] \nonumber \end{align}

Both equations can be integrated directly, giving

(4.5a) \begin{align} \frac {1}{{Pe}}\frac {\mathrm{d}C^{-}}{\mathrm{d}x}-C^{-}U^{-}=\beta ^{-}&, \hspace {10pt} x\in (-\infty ,0), \\[-12pt] \nonumber \end{align}

(4.5b) \begin{align} U^{-}=\frac {A}{{Pe}}\frac {\mathrm{d}C^{-}}{\mathrm{d}x}+\alpha ^{-}&, \hspace {10pt} x\in (-\infty ,0), \\[10pt] \nonumber \end{align}

where $\alpha ^{-},\beta ^{-}$ are integration constants. Substituting (4.5b ) into (4.5a ), we obtain a first-order separable differential equation

(4.6) \begin{align} \frac {\mathrm{d}C^{-}}{\mathrm{d}x}={Pe}\,\frac {\beta ^{-}+\alpha ^{-}C^{-}}{1-AC^{-}} \end{align}

with the implicit solution

(4.7a) \begin{align} \alpha ^{-}\left ({Pe}\, x+\gamma ^{-}\right )=\left (1+\frac {\beta ^{-}}{\alpha ^{-}}A\right )\ln |\beta ^{-}+\alpha ^{-}C^{-}|-AC^{-}, \hspace {10pt} x\in (-\infty ,0), \end{align}

where $\gamma ^{-}$ is another integration constant. Now, we take the opposite direction and substitute (4.6) back into (4.5b ), which yields the expression for the fluid velocity

(4.7b) \begin{align} U^{-}=\frac {\alpha ^{-}+\beta ^{-} A}{1-AC^{-}}, \hspace {10pt} x\in (-\infty ,0). \end{align}

In deriving the solutions (4.7), we assumed that $\alpha ^{-}\neq 0$ , which can be justified as follows. If we take $\alpha ^{-}=0$ , then by integrating (4.6), we obtain a quadratic equation for $C^{-}$ :

(4.8) \begin{align} \frac {A}{2}\left (C^{-}\right )^2-C^{-}+{Pe}(\beta ^{-} x + \gamma ^{-})=0. \end{align}

The solution to (4.8) is

(4.9) \begin{align} C^{-}(x)=\frac {1\pm \sqrt {\gamma '-2A{Pe}\beta ^{-} x}}{A}, \hspace {10pt} x\in (-\infty ,0), \end{align}

where we have included all constants under $\gamma '$ . Now, considering the limit $x\to -\infty$ , the term $(2A{Pe}\beta ^{-} x)$ diverges, which contradicts the assumption (4.2a ). Therefore, the only possible correction is to set $\beta ^{-}=0$ , which reduces $C^{-}$ to the constant

(4.10) \begin{align} C^{-}=\frac {1\pm \sqrt {\gamma '}}{A}. \end{align}

However, looking back at (4.5b ), the flow velocity $U^{-}$ would be in that case zero for the whole upstream interval $x\in (-\infty ,0)$ , which is also a contradiction to (4.2b ). Therefore, the constant $\alpha ^{-}$ must be non-zero.

Moreover, we can determine the first two integration constants $\alpha ^{-},\beta ^{-}$ from the inlet conditions (4.2). For the limit $x\to -\infty$ , the left-hand side of (4.7a ) approaches $\pm \infty$ , depending on the sign of $\alpha ^{-}$ . However, on the right-hand side, the only term capable of going to infinity is $\ln |\beta ^{-}+\alpha ^{-}C^{-}|$ . Therefore, to satisfy equality for $x\to -\infty$ , the following limit must apply:

(4.11) \begin{align} C^{-}(x)\to -\frac {\beta ^{-}}{\alpha ^{-}}, \hspace {10pt} x\to -\infty . \end{align}

Then, from the inlet concentration value (4.2a ), we obtain

(4.12) \begin{align} \alpha ^{-}=-\beta ^{-}. \end{align}

Applying the limit $x\to -\infty$ also to the expression (4.7b ) yields

(4.13) \begin{align} 1=\frac {\alpha ^{-}+\beta ^{-} A}{1-A}, \end{align}

which together with the previously derived relation (4.12) gives

(4.14) \begin{align} \alpha ^{-}=1, \hspace {7pt} \beta ^{-}=-1. \end{align}

These values can be substituted back to the upstream solution (4.7), giving its final form

(4.15a) \begin{align} \textit{Pe} x+\gamma ^{-}=\left (1-A\right )\ln |C^{-}-1|-AC^{-}, \hspace {10pt} x\in (-\infty ,0), \\[-12pt] \nonumber \end{align}

(4.15b) \begin{align} U^{-}=\frac {1-A}{1-AC^{-}}, \hspace {10pt} x\in (-\infty ,0). \\[10pt] \nonumber \end{align}

Let us note that this result holds only under the assumption that $A\neq 1$ , which is always satisfied. In fact, $A\lt 1$ . To see this, consider the definition of $A$ , the condition $A\lt 1$ is equivalent to

(4.16) \begin{align} \left (\frac {1}{\rho _2^T}-\frac {1}{\rho _1^T}\right )c_\infty \lt 1. \end{align}

Provided that $\rho _1^T,\rho _2^T,c_\infty \gt 0$ , (4.16) can be rewritten as

(4.17) \begin{align} 1-\frac {\rho _2^T}{c_\infty }\lt \frac {\rho _2^T}{\rho _1^T}. \end{align}

We remind that the true density $\rho _2^T$ corresponds to the solute, as does the characteristic concentration $c_\infty$ . The true density is defined only through the volume related to the solute, whereas the concentration is defined throughout the entire volume. Since we chose the $c_\infty$ as the average concentration within the mixture, it must hold that $c_\infty \lt \rho _2^T$ . Therefore, the entire left-hand side of (4.17) is negative and the inequality is always satisfied, meaning the assumption $A\lt 1$ is justified. Finally note that in the limit $x\to -\infty$ , the left-hand side of (4.15a ) also approaches $-\infty$ . Therefore, the term $ (1-A )\ln |C^{-}-1|$ must also approach minus infinity, which is satisfied as $A\lt 1$ .

4.1.2. Downstream equations

So far, we have specified the constants $\alpha ^{-}$ and $\beta ^{-}$ of the solution (4.7). To specify the third yet undetermined constant $\gamma ^{-}$ , we would need yet another inlet condition. However, as shown in the following, we do not need the value of $\gamma ^{-}$ to determine the filter boundary conditions at $x\in \{0,1\}$ to have a well-defined problem.

In the same manner as in the upstream case, we arrive at the downstream equations:

(4.18a) \begin{align} \frac {{d}}{\mathrm{d}x}\left (\frac {1}{{Pe}}\frac {\mathrm{d}C^{+}}{\mathrm{d}x}-C^{+}U^{+}\right )=0&, \hspace {10pt} x\in (1,+\infty ), \\[-12pt] \nonumber \end{align}

(4.18b) \begin{align} \frac {\mathrm{d}U^{+}}{\mathrm{d}x}=\frac {A}{{Pe}}\frac {\mathrm{d}^2C^{+}}{\mathrm{d} x^2}&, \hspace {10pt} x\in (1,+\infty ). \\[10pt] \nonumber \end{align}

By integrating and combining them, we obtain

(4.19) \begin{align} \frac {1}{{Pe}}\frac {\mathrm{d}C^{+}}{\mathrm{d}x}-\frac {A}{{Pe}}C^{+}\frac {\mathrm{d}C^{+}}{\mathrm{d}x}-\alpha ^{+}C^{+}=\beta ^{+}, \end{align}

where $\alpha ^{+}$ and $\beta ^{+}$ are integration constants defined analogously to the upstream case. Assuming that $(\mathrm{d}C^{+}/\mathrm{d}x)\neq 0$ , solution to this equation is again in the implicit form (and same as the upstream solution):

(4.20) \begin{align} \alpha ^{+}\left ({Pe} x+\gamma ^{+}\right )=\left (1+\frac {\beta ^{+}}{\alpha ^{+}}A\right )\ln |\beta ^{+}+\alpha ^{+}C^{+}|-AC^{+}, \end{align}

with $\gamma ^{+}$ being the last integration constant. However, we show that this solution is not admissible for the upstream problem. To see this, we consider the limit $x\to +\infty$ in (4.19). The first two terms on the left-hand side approach zero due to $(\mathrm{d}C^{+}/\mathrm{d}x)\to 0$ , giving the limit value for $C^{+}$ :

(4.21) \begin{align} C^{+}(x)\to -\frac {\beta ^{+}}{\alpha ^{+}}, \hspace {10pt} x\to +\infty . \end{align}

Also applying the limit to the implicit solution (4.20), we see that the left-hand side approaches $+\infty$ . The positive sign is a result of the limit

(4.22) \begin{align} U^{+}(x)\to \alpha ^{+}, \hspace {10pt} x\to +\infty , \end{align}

which is a consequence of $U^{+}= [(A/{Pe}) (\mathrm{d}C^{+}/{\mathrm{d}x})]+\alpha ^{+}$ and $ {\mathrm{d}C^{+}}/{\mathrm{d}x}\to 0$ as $x\to +\infty$ . Additionally, because $U^{+}$ is the fluid velocity behind the filter, the value $\alpha ^{+}$ must be positive. However, on the right-hand side of (4.20), the term $\ln |\beta ^{+}+\alpha ^{+}C^{+}|$ approaches $-\infty$ as $x\to +\infty$ , due to (4.21). Therefore, we must consider

(4.23) \begin{align} 1+\frac {\beta +}{\alpha +}A\lt 0, \end{align}

which can be rewritten as

(4.24) \begin{align} \alpha ^{+}\lt -\beta ^{+} A. \end{align}

Going back to the limit value of $C^{+}$ , we have

(4.25) \begin{align} C^{+}(x)\to -\frac {\beta ^{+}}{\alpha ^{+}}\gt \frac {-\beta ^{+}}{-\beta ^{+} A}=\frac {1}{A}\gt 1, \hspace {10pt} x\to +\infty , \end{align}

where we used the previously derived condition $A\lt 1$ . This means that the outlet concentration $C^{+}$ is higher than the inlet concentration $C^{-}$ , which is physically incorrect, considering we have no sources of matter within the filter. Thus, the implicit solution (4.20) for the downstream concentration $C^{+}$ is indeed not admissible.

Note however that the inadmissible solution was derived based on the assumption $ {\mathrm{d}C^{+}}/{\mathrm{d}x}\neq 0$ . Returning now to the omitted special case, we can see that (4.5) can be solved supposing that $C^{+}$ and $U^{+}$ are constants. To match the previous notation, we denote them as

(4.26a) \begin{align} C^{+}&=-\frac {\beta ^{+}}{\alpha ^{+}}, \\[-12pt] \nonumber \end{align}

(4.26b) \begin{align} U^{+}&=\alpha ^{+}. \\[10pt] \nonumber \end{align}

However, since these constants are outgoing values of the concentration and fluid velocity from the filter, they remain unknown and must be determined as a part of the final solution to (4.1).

4.1.3. Final boundary conditions

We derived the solutions (4.7) and (4.26) for both reservoirs before and after the filter. Despite not having determined the constants $\alpha ^{+},\beta ^{+},\gamma ^{\pm }$ , we can use the knowledge of solutions to derive the boundary conditions for the filter itself and, hence, simplify the problem. For this purpose, we need a condition at the filter boundary $x\in \{0,1\}$ connecting the external variables $C^{-},C^{+},U^{-},U^{+}$ with the internal variables $C,U$ . For this, it is sufficient to take the following relations:

(4.27) \begin{align} C^{-}=\frac {C}{\phi }, \; U^{-}=U&, \hspace {10pt} x=0, \\[-12pt] \nonumber \end{align}

(4.28) \begin{align} C^{+}=\frac {C}{\phi }, \; U^{+}=U&, \hspace {10pt} x=1, \\[10pt] \nonumber \end{align}

which guarantee the continuity of solute concentration and fluid velocity at the filter boundaries. From a physical point of view, we also need to impose the continuity conditions of concentration flux

(4.29) \begin{align} \frac {1}{{Pe}}\frac {\mathrm{d}C^{-}}{\mathrm{d}x}-C^{-}U^{-}=\tilde {D}\frac {\mathrm{d}C}{\mathrm{d}x}-\frac {C}{\phi }\left (U+\tilde {D}\frac {\mathrm{d}\phi }{\mathrm{d}x}\right )&, \hspace {10pt} x=0, \\[-12pt] \nonumber \end{align}

(4.30) \begin{align} \frac {1}{{Pe}}\frac {\mathrm{d}C^{+}}{\mathrm{d}x}-C^{+}U^{+}=\tilde {D}\frac {\mathrm{d}C}{\mathrm{d}x}-\frac {C}{\phi }\left (U+\tilde {D}\frac {\mathrm{d}\phi }{\mathrm{d}x}\right )&, \hspace {10pt} x=1, \\[10pt] \nonumber \end{align}

which can be derived by comparing (4.1a ) with (4.4a ) or (4.18a ).

Using these relations, we can finally derive the filter boundary conditions for the set of equations (4.1). We take the concentration flux continuity condition (4.29), whose left-hand side is from (4.5a ) equal to $\beta ^{-}$ . Knowing that $\beta ^{-}=-1$ , we obtain the final condition on the left filter boundary:

(4.31a) \begin{align} \tilde {D}\frac {\mathrm{d}C}{\mathrm{d}x}-\frac {C}{\phi }\left (U+\tilde {D}\frac {\mathrm{d}\phi }{\mathrm{d}x}\right )=-1, \hspace {10pt} x=0. \end{align}

Similarly for the right boundary condition, we start from (4.30). Considering that $C^{+}$ is constant, we obtain

(4.31b) \begin{align} -C^{+}U^{+}=\tilde {D}\frac {\mathrm{d}C}{\mathrm{d}x}-\frac {C}{\phi }\left (U+\tilde {D}\frac {\mathrm{d}\phi }{\mathrm{d}x}\right )\!, \hspace {10pt} x=1. \end{align}

This can be further amended using (4.28), yielding the final boundary condition on the right side of the filter,

(4.31c) \begin{align} \frac {\mathrm{d}C}{\mathrm{d}x}-\frac {C}{\phi }\frac {\mathrm{d}\phi }{\mathrm{d}x}=0, \hspace {10pt} x=1. \end{align}

So far, we have derived two Robin boundary conditions for the solute concentration equation (4.1a ). To complete the system, we also need a boundary condition for the fluid velocity equation (4.1b ). It can be obtained directly from the expression (4.7b ), using the derived constants (4.14) together with the relations (4.27) as

(4.31d) \begin{align} U=\frac {1-A}{1-\frac {AC}{\phi }}, \hspace {10pt} x=0. \end{align}

That the above-mentioned boundary condition is well defined is clear from the definition of $C$ in (3.1), which, considering the already derived observation (3.22), leaves the denominator non-zero.

To summarise, we obtained the three desired boundary conditions at the filter boundary $x\in \{0,1\}$ , which completes the well-posed description of the unidirectionally graded filter (4.1). Moreover, this allows us to separate the filtering task from the entire region with reservoirs and limit ourselves to the solution inside the filter only.

4.2.1. Leading-order solution

We substitute the expansions (4.33) and (4.34) into the set of (4.1) with boundary conditions (4.31a ), (4.31c ) and (4.31d ), and solve order by order in $\epsilon$ . For the leading order, we obtain the equations

(4.39a) \begin{align} \frac {{\rm d}}{\mathrm{d}x}\left (\tilde {D_0} \frac {\mathrm{d}C_0}{\mathrm{d}x}-\frac {C_0}{\phi _0}U_0\right )=f_0 C_0, \hspace {10pt} x\in (0,1), \\[-12pt] \nonumber \end{align}

(4.39b) \begin{align} \frac {\mathrm{d}U_0}{\mathrm{d}x}=\frac {A}{{Pe}}\frac {\mathrm{d}^2 C_0}{\mathrm{d} x^2}, \hspace {10pt} x\in (0,1) \\[10pt] \nonumber \end{align}

with the respective boundary conditions

(4.40a) \begin{align} \tilde {D_0}\frac {\mathrm{d}C_0}{\mathrm{d}x}-\frac {C_0}{\phi _0}U_0=-1, \hspace {10pt} x=0, \\[-12pt] \nonumber \end{align}

(4.40b) \begin{align} U_0=\frac {1-A}{1-\frac {A C_0}{\phi _0}}, \hspace {10pt} x=0, \\[-12pt] \nonumber \end{align}

(4.40c) \begin{align} \frac {\mathrm{d}C_0}{\mathrm{d}x}=0, \hspace {10pt} x=1. \\[10pt] \nonumber \end{align}

After some algebra, the set of (4.39) leads to a second-order differential equation with polynomial coefficients. A solution to such an equation is in the form of hypergeometric functions and is therefore not suitable for our purpose to obtain a closed form solution, in which higher-order solutions are expressed in terms of the lower orders.

Therefore, to progress further, we resort to one further simplification, that is, the constant $A$ is of order $\textit {O}(\epsilon )$ , specifically

(4.41) \begin{align} A= \epsilon \bar {A}. \end{align}

Note that as $A= [({1}/{\rho _2^T})- ({1}/{\rho _1^T})]c_\infty$ , its smallness is justified, for example, when the solution is dilute, and/or the true concentrations of the solute and solvent are not different by orders of magnitude.

Equation (4.39b ) with boundary condition (4.40b ) then decouples into

(4.42a) \begin{align} \frac {\mathrm{d}U_0}{\mathrm{d}x}=0&, \hspace {10pt} x\in (0,1), \\[-12pt] \nonumber \end{align}

(4.42b) \begin{align} U_0=1&, \hspace {10pt} x=0 \\[10pt] \nonumber \end{align}

with the simple solution

(4.43) \begin{align} U_0=1, \hspace {10pt} x\in (0,1). \end{align}

The rest of the equations then form a second separate system,

(4.44a) \begin{align} \tilde {D_0}\frac {\mathrm{d}^2 C_0}{\mathrm{d} x^2}-\frac {1}{\phi _0}\frac {\mathrm{d}C_0}{\mathrm{d}x}-f_0 C_0=0&, \hspace {10pt} x\in (0,1), \\[-12pt] \nonumber \end{align}

(4.44b) \begin{align} \tilde {D_0}\frac {\mathrm{d}C_0}{\mathrm{d}x}-\frac {C_0}{\phi _0}=-1&, \hspace {10pt} x=0, \\[-12pt] \nonumber \end{align}

(4.44c) \begin{align} \frac {\mathrm{d}C_0}{\mathrm{d}x}=0&, \hspace {10pt} x=1, \\[10pt] \nonumber \end{align}

being a second-order differential equation with constant coefficients. Therefore, the leading-order solution reads

(4.45) \begin{align} C_0(x)=2\alpha \phi _0 e^{ax}\left [b\cosh (ab(1-x))+\sinh (ab(1-x))\right ], \end{align}

where

(4.46) \begin{align} a=\frac {1}{2\tilde {D_0}\phi _0}, \hspace {5pt} b=\sqrt {1+\frac {f_0}{a^2 \tilde {D_0}}} \end{align}

and

(4.47) \begin{align} \alpha =\frac {1}{(b^2+1)\sinh (ab)+2b\cosh (ab)}. \end{align}

4.2.2. First-order solution

The leading-order solution corresponds to the case when $A=0$ due to our simplifying assumption $A=\textit {O}(\epsilon )$ . Therefore, the first-order solution represents a correction due to the effect of mixture of two incompressible components.

Again, substituting the expansions (4.33) and (4.34) into the set (4.1) with boundary conditions (4.31a ), (4.31c ), (4.31d ) and equating terms of order $\textit {O}(\epsilon )$ , we obtain

(4.48a) \begin{align} \tilde {D_0}\frac {\mathrm{d}^2 C_1}{\mathrm{d} x^2}-\frac {1}{\phi _0}\frac {\mathrm{d}C_1}{\mathrm{d}x}-f_0 C_1 & =-\frac {{d}}{\mathrm{d}x}\left (\tilde {D_1}\frac {\mathrm{d}C_0}{\mathrm{d}x}-\frac {C_0}{\phi _0}\left (U_1+\tilde {D_0}\frac {\mathrm{d}\phi _1}{\mathrm{d}x}-\frac {\phi _1}{\phi _0}\right )\right ) \nonumber \\[5pt] & \quad +f_1 C_0, \hspace {10pt} x\in (0,1), \\[-12pt] \nonumber \end{align}

(4.48b) \begin{align} \frac {\mathrm{d}U_1}{\mathrm{d}x} & =\frac {\bar {A}}{{Pe}}\frac {\mathrm{d}^2 C_0}{\mathrm{d} x^2}, \hspace {10pt} x\in (0,1) . \\[10pt] \nonumber \end{align}

The boundary conditions are then

(4.49a) \begin{align} \tilde {D_0} \frac {\mathrm{d}C_1}{\mathrm{d}x}-\frac {C_1}{\phi _0}=-\tilde {D_1}\frac {\mathrm{d}C_0}{\mathrm{d}x}+\frac {C_0}{\phi _0}\left (U_1+\tilde {D_0}\frac {\mathrm{d}\phi _1}{\mathrm{d}x}-\frac {\phi _1}{\phi _0}\right )&, \hspace {10pt} x=0, \\[-12pt] \nonumber \end{align}

(4.49b) \begin{align} U_1=\bar {A}\left (\frac {C_0}{\phi _0}-1\right )&, \hspace {10pt} x=0, \\[-12pt] \nonumber \end{align}

(4.49c) \begin{align} \frac {\mathrm{d}C_1}{\mathrm{d}x}=\frac {C_0}{\phi _0}\frac {\mathrm{d}\phi _1}{\mathrm{d}x}&, \hspace {10pt} x=1. \\[10pt] \nonumber \end{align}

The flow equation (4.48b ) can be directly integrated, giving the solution

(4.50) \begin{align} U_1(x)=\frac {\bar {A}}{{Pe}}2\alpha \phi _0 a\big (1-b^2\big )e^{ax}\sinh (ab(1-x))+\theta , \end{align}

where $\theta$ follows from the boundary condition (4.49b ), yielding

(4.51) \begin{align} \theta =\bar {A}\left (\left (1-\frac {\phi _0}{{Pe}}a\left (1-b^2\right )\right )2\alpha \sinh (ab)+2\alpha b\cosh (ab)-1\right )\!. \end{align}

We note that with the velocity $U_1$ determined, all the terms on the right-hand side of (4.48a ) are explicit known functions of space $x$ . Therefore, we are left with an inhomogeneous second-order differential equation with constant coefficients. The solution can be readily found using the method of variation of parameters and reads

(4.52) \begin{align} C_1(x) & =\frac {e^{ax}}{2ab\tilde {D_0}}\left [\left (\int G(x)e^{-(a+ab)x}\,{\rm d}\textit{x}+\gamma _1\right )e^{abx} \right. \nonumber \\[5pt]& \quad \left. -\left (\int G(x)e^{-(a-ab)x}\,{\rm d}\textit{x}+\gamma _2\right )e^{-abx}\right ]\!, \end{align}

where $G(x)$ is defined as the right-hand side of (4.48a ). Further, the constants satisfy

(4.53) \begin{align} \gamma _1& =\frac {\alpha }{2}\left (-(b+1)^2 e^{ab}\left (\beta _1(1)-\beta _1(0)\right )-\frac {(b+1)^3}{b-1}e^{ab}\beta _2(0)+\left (b^2-1\right )e^{-ab}\beta _2(1) \right . \nonumber \\[5pt]&\quad \left . +2\alpha b\frac {b+1}{a}\frac {\mathrm{d}\phi _1}{\mathrm{d}x}\left (1\right )-2\frac {(b+1)^2}{b-1}\phi _0 e^{ab}N(0)\right )+\frac {b+1}{b-1}\beta _2(0)-\beta _1(0) \nonumber \\[5pt]& \quad +2\frac {\phi _0}{b-1}N(0), \\[-12pt] \nonumber \end{align}

(4.54) \begin{align} \gamma _2& =\frac {\alpha }{2}\,\Bigg ({-}\big (b^2-1\big ) e^{ab}\left (\beta _1(1)-\beta _1(0)\right )-(b+1)^2e^{ab}\beta _2(0)+(b-1)^2e^{-ab}\beta _2(1) \nonumber \\[5pt]&\quad +2\alpha b\frac {b-1}{a}\frac {\mathrm{d}\phi _1}{\mathrm{d}x}\left (1\right )-2(b+1)\phi _0 e^{ab}N(0)\Bigg ), \\[10pt] \nonumber \end{align}

with

(4.55) \begin{align} \beta _1(x)=\frac {1}{2ab\tilde {D_0}}\int G(x)e^{-(a+ab)x}\,\mathrm{d}x + \gamma _1, \\[-12pt] \nonumber \end{align}

(4.56) \begin{align} \beta _2(x)=-\frac {1}{2ab\tilde {D_0}}\int G(x)e^{-(a-ab)x}\,\mathrm{d}x + \gamma _2, \\[10pt] \nonumber \end{align}

and where we denoted the function $N(x)$ as the right-hand side of (4.49a ).

4.3.1. Square lattice

To evaluate the solution, we also need to specify the effective coefficients $\tilde {D}(\phi )$ and $f(\phi )$ which, as pointed out earlier, are dependent on the filter microscale geometry. The adsorption coefficient $f(\phi )$ is for spherical inclusions (as a filter geometry) in two dimensions from (3.39) equal to

(4.58) \begin{align} f(\phi )=k\frac {2\pi }{\phi }\sqrt {\frac {1-\phi }{\pi }}. \end{align}

For the diffusion coefficient $\tilde {D}(\phi )$ , we start from the closed-form approximation for a square lattice derived by Jóhannesson & Halle (Reference Jóhannesson and Halle1996):

(4.59) \begin{align} \tilde {D}(\phi )=\frac {1}{{Pe}\,\phi }\left (1-\frac {2(1-\phi )}{2-\phi -0.3058(1-\phi )^4}\right )\!, \hspace {10pt} \phi \gt 0.3, \end{align}

where the volume fraction value on the right indicates the values for which the approximation is valid. As shown in figure 2, the approximation provides very good results for porosity $\phi \gt 0.3$ .

Figure 3.

Numerically obtained (a) fluid-volume averaged solute concentration $C^{\kern-1pt f}=C/\phi$ and (b) fluid velocity $U=\phi U^{\kern-1pt f}$ for three gradients $m=-0.3,0,0.3$ . The evaluation is performed for a two-dimensional square lattice with the parameters $\phi _0=0.75,\, A=0.6,\, {Pe}=3$ and $k=1$ . Note that for the special case of divergence free flow field, $A=0$ , the solution is $U=1$ and, hence, it can be viewed as a measure of the effect of the novel term, being a result of treating the filtration problem as a true mixture.

Using the relations for the effective macroscale coefficients, we evaluate the solutions $C^{\kern-1pt f}$ and $U$ for various values of parameters $A,\,{Pe},\,k$ and $m$ in figure 3. We opt for the fluid-volume averaged quantity $C^{\kern-1pt f}$ , which needs to be monotonic and flatten out to zero derivative when leaving the filter. The fluid flow velocity, however, is chosen to be $U=\phi U^{\kern-1pt f}$ as it allows a simple comparison to the special case $A=0$ (fluid flow being divergence free), where the solution is known to be $U=1$ , see the following. As expected, the concentration value decreases with increasing $x$ , as the solute particles are being filtered out throughout the filter. We also observe that with increasing porosity gradient, the rate of concentration change decreases. This behaviour is also to be expected, because as the porosity gradient increases, more contaminant is filtered out towards the front of the filter, where the filter obstacle size is larger. Therefore, the concentration near the front should decrease faster. The same logic applies to the end of the filter, where thus, the overall concentration profile becomes more uniform. Note that such observations hold for $C^{\kern-1pt f}$ , not $C$ .

In addition, note that the solute concentration $C^{\kern-1pt f}$ does not enter the filter with the (normalised) free-flow upstream value, that is, $C^{\kern-1pt f}(x=0)\neq 1$ . The reason is that, as we have shown already in § 4.1, the experimentally controlled influx conditions at $x=-\infty$ translate into non-trivial Robin boundary conditions at the entry to the filter $x=0$ , (4.31a ). As can be seen from the boundary condition (4.31a ) and also from the numerical solutions in figure 3, the filter presence and geometry significantly affect the composition of the mixture even before (upstream) the filter and entails various entry values of the solute concentration.

Figure 4.

Numerically obtained concentration profiles for linear porosity (4.57) for $A=0$ and $A=0.6$ , hence, highlighting the effect of the new term. The evaluation is performed for a two-dimensional square lattice with the parameters $\phi _0=0.75,\, {Pe}=3,\, k=1$ and $m=-0.3,0.3$ . The lower two curves (solid and dashed) correspond to positive porosity gradient $m=0.3$ , while the upper two to a negative porosity gradient $m=-0.3$ . Note that the inlet solute concentration is affected (being larger for $A\neq 0$ ) as well as the amount of the particles being filtered out (again, being larger for $A\neq 0$ ).

To examine the effect of the non-zero mixture velocity divergence, that is, the right-hand side of (4.1b ), we already mentioned that we plot $U$ instead of $U^{\kern-1pt f}$ in figure 3(b). The velocity profile is affected significantly, as for $A=0$ , the solution from (4.1b ) to (4.31d ) is $U=1$ , regardless of the used porosity profile. In addition, we calculate and plot the effect on the concentration profile, see figure 4. We observe that increasing $A$ leads to a more pronounced concentration change during the filtration. Therefore, this represents the first major qualitatively different result of our filtration model.

4.3.2. Hexagonal lattice

Until now, we have only examined the filtration process in dependence on various porosity profiles $\phi (x)$ and parameters $A,{Pe},k$ . However, the effective coefficients $D(\phi ),f(\phi )$ from their definitions (3.37), (3.38) depend not only on the porosity $\phi$ , but also on the geometry of the filter. This can be seen from the cell problem (3.25) and constitutes one of the advantages of an asymptotic approach to deriving the macroscale equations, as the method provides a link between microscale and macroscale description. The solution to the cell problem explicitly reflects the geometry of the unit cell $\omega$ . To demonstrate this effect and link to the microscale, we introduce a new filter layout in the form of a hexagonal lattice in two dimensions (with the same spherical inclusions). We represent the hexagonal lattice by unit cells $\omega = [-3^{-1/4},3^{-1/4} ]\times [-3^{1/4},3^{1/4} ]$ with one solid disk placed in the centre and a quarter of a disk placed in each corner of the cell. By this definition, we achieve a uniform distribution of the disks, while keeping the cell size at $|\omega |=1$ .

To evaluate the solution for such layout, we first need to modify the effective coefficients $D,f$ for this specific geometry. Since every hexagonal cell contains two disks in total, we can rewrite the adsorption coefficient analogous to (3.39) as

(4.60) \begin{align} f_{\textit{H}}(\phi )=k\frac {4\pi }{\phi }\sqrt {\frac {1-\phi }{2\pi }}, \end{align}

where the index H denotes the hexagonal lattice. Next, we take the diffusion coefficient from Jóhannesson & Halle (Reference Jóhannesson and Halle1996) using the same steps as for the square lattice in (4.59), yielding

(4.61) \begin{align} \tilde {D}_{\textit{H}}(\phi )=\frac {1}{{Pe}\,\phi }\left (1-\frac {2(1-\phi )}{2-\phi -0.07542(1-\phi )^6}\right )\!, \hspace {10pt} \phi \gt 0.2. \end{align}

Having specified the above-mentioned coefficients, we evaluate concentration profiles for the linear porosity (4.57) in figure 5. As expected, we observe slightly different results for both lattices, despite using the same porosity values. Specifically, we see that the hexagonal layout causes higher solute adsorption at the front of the filter. This may be a consequence of a more uniform distribution of the balls (disks) within the filter, in comparison to the square lattice. However, to confirm this phenomenon, further analysis would be required.

Figure 5.

Numerically obtained concentration profiles for linear porosity (4.57) for square and hexagonal lattices in two dimensions. The parameters used are $\phi _0=0.75$ , $A=0.6,\, {Pe}=3,\, k=1$ and $m=-0.3,0.3$ . The effect of microscale geometry and the capability of the asymptotic method to reflect these effects in the derived macroscale model is apparent.

4.3.3. Asymptotic approximation

Taking the advantage of having both numerical and approximate solution, we check their validity by mutual comparison. In particular, we compare the asymptotic approximation of the solution for slowly varying porosity with the numerical results obtained here.

For this purpose, we once again use the liner porosity profile $\phi (x)=\phi _0+m(x-0.5)$ . Comparing it with the asymptotic expansion (4.33c ), we require

(4.62) \begin{align} m=\epsilon , \hspace {10pt} \phi _1(x)=x-0.5. \end{align}

We also need to specify the effective coefficients $\tilde {D}_0,\tilde {D}_1$ and $f_0,f_1$ . The adsorption coefficients $f_0$ and $f_1$ can be determined from prescription (4.38) using the previously derived value (4.58) as

(4.63) \begin{align} f_0&=k\frac {2\pi }{\phi _0}\sqrt {\frac {1-\phi _0}{\pi }}, \\[-12pt] \nonumber \end{align}

(4.64) \begin{align} f_1(x)&=k\phi _1(x)\frac {\phi _0-2}{\phi _0^2}\sqrt {\frac {\pi }{1-\phi _0}}. \\[10pt] \nonumber \end{align}

In the same manner, we obtain the diffusion coefficients from (4.37) using (4.59) as

(4.65) \begin{align} \tilde {D}_0&=\frac {1}{\phi _0}\left (1-\frac {2(1-\phi _0)}{2-\phi _0-0.3058(1-\phi _0)^4}\right )\!, \\[-12pt] \nonumber \end{align}

(4.66) \begin{align} \tilde {D}_1(x)&=\frac {\phi _1(x)}{{Pe}\,\phi _0^2}\left (-1+\frac {2(1-\phi _0)}{2-\phi _0-0.3058(1-\phi _0)^4}+\frac {\left (2+1.8348(1-\phi _0)^4\right )\phi _0}{\left (2-\phi _0-0.3058(1-\phi _0)^4\right )^2}\right )\!. \\[0pt] \nonumber \end{align}

With specified effective coefficients, we evaluate the asymptotic solution in figure 6. We see that the first-order approximation gives very good results in comparison with the numerical solution. We note that the deviation is slightly larger for the steeper gradient $m=-0.3$ in figure 6(a), as should be expected because $\epsilon =m$ .

Figure 6.

(a,b) Asymptotic and numerical concentration profiles for linear porosity given by (4.57). The evaluation is performed for a two-dimensional square lattice with the parameters $\phi _0=0.75$ , $A=0.6, \textit{Pe}=3,\, k=1,\, m=-0.3$ in panel (a) and $m=-0.1$ in panel (b).

This close form solution allows further qualitative understanding of the filtration’s dynamics. In particular, we may inspect the effect of key model parameters: $ \textit{Pe}$ , $A$ or the effect of porosity ( $\phi _{0}$ or $m=\epsilon$ ), see § 6.