1 Introduction
Filtration is a vital process in many industries, such as kidney dialysis (Lonsdale Reference Lonsdale1982), air purification (Barhate & Ramakrishna Reference Barhate and Ramakrishna2007), waste water treatment (Vandevivere, Bianchi & Verstraete Reference Vandevivere, Bianchi and Verstraete1998) and beer production (Fillaudeau & Carrère Reference Fillaudeau and Carrère2002). Although the industrial applications may vary widely, the main goal is often the same: to maximize the removal of contaminants or particulates entrained within the fluid that passes through the filter. Since a single experiment may take hours to complete, mathematical modelling provides a valuable tool for assisting with filter design.
Broadly speaking, filters are composed of either material containing a series of pores (e.g. track-etched membranes) or a fibrous mesh (e.g. so-called depth filters). In the former case, if the contaminants are larger than the pores, filtration occurs at the pore entrance via size exclusion. A cake layer of particulates builds up on the entrance, making it more difficult for fluid to pass through the filter. When the contaminants are smaller than the pores, they are able to penetrate into the pores, where they may become trapped in the tortuous pore network and adhere to the pore walls. This reduces the available space through which contaminated fluid can flow, and again makes it more difficult for fluid to pass through the filter (Griffiths, Kumar & Stewart Reference Griffiths, Kumar and Stewart2014). Filters composed of a fibrous mesh provide an internal structure of obstacles to which the contaminants can adsorb. As with filters composed of pores, the physical trapping of particulates within the filter will cause a reduction in the available space through which the fluid can flow. Since the contaminants are removed as the fluid flows into the filter depth, their concentration decreases with depth and so this constriction effect is usually more significant towards the filter entrance (Datta & Redner Reference Datta and Redner1998). Over time, this constriction may ultimately become a blockage, halting filtration. To prolong lifespan, filters are often designed such that their porosity decreases with depth, allowing areas of the filter away from the entrance to have trapped more contaminants by the time of blocking (Anderson Reference Anderson1951; Burggraaf & Keizer Reference Burggraaf and Keizer1991; Dickerson et al. Reference Dickerson, Monnin, Rickle, Borer, Stuart, Velu, Haberkamp and Graber2005; Barg, Koch & Grathwohl Reference Barg, Koch and Grathwohl2009; Vida Simiti et al. Reference Vida Simiti, Jumate, Moldovan, Thalmaier and Sechel2012). These are known as porosity-graded filters. Moreover, an initially homogenous filter will become graded over time as the local porosity changes due to local contaminant trapping.
While mathematical methods can provide a cost-effective way to predict solute transport without experimentally sweeping through parameter space, the explicit coupling of solute transport with filter evolution yields a complicated moving boundary problem. This significantly increases the mathematical and computational effort required to solve the system and, as such, many mathematical models of filtration consider only one of these mechanisms. That is, either the problem of solute transport, or the problem of filter evolution. For example, in Griffiths et al. (Reference Griffiths, Kumar and Stewart2014), the filter evolution is considered by simulating a microscale model (on the lengthscale of pore size) for the blocking of a network of pores from a fundamental mechanistic level. While the results are consistent with the standard set of constitutive macroscale laws (on the lengthscale of filter size) used to predict fluid throughput (often fitted with experimental data) (Bowen, Calvo & Hernandez Reference Bowen, Calvo and Hernandez1995), the flow problem is simplified by assuming Poiseuille flow through the constricting cylindrical pores. This idea is generalized to explore the effect of tapered pores and filters composed of multiple membrane layers in Griffiths, Kumar & Stewart (Reference Griffiths, Kumar and Stewart2016). A different approach to the problem of solute transport past sinks is used in Chernyavsky et al. (Reference Chernyavsky, Leach, Dryden and Jensen2011), where the authors considered blood flow and nutrient delivery in the placenta. Here, the transport of nutrient governed by diffusion and advection with unidirectional flow past randomly placed point sinks in the placenta is considered in one dimension. The effect of placental growth and the nature of the local flow are neglected to focus on the important features of the paper, such as quantifying the error in upscaling techniques. The authors do this by analytically determining the concentration field for different statistical distributions of the sinks using homogenization techniques and comparing these to exact solutions. In Krehel, Muntean & Knabner (Reference Krehel, Muntean and Knabner2015), the authors also use homogenization techniques to consider colloid particles diffusing around a fluid domain past a periodic array of solid obstacles. The particles are trapped on the surface of these obstacles, where they may aggregate.
Mathematical homogenization provides an upscaling method whereby the fundamental behaviour of concern, such as the contaminant removal rate, is captured without the computational expense of globally calculating the microscale behaviour. This is achieved by averaging the microscale variation while retaining the macroscale variation (see, for example, Bensoussan, Lions & Papanicolaou Reference Bensoussan, Lions and Papanicolaou1978, Sánchez-Palencia Reference Sánchez-Palencia1980, Mei & Vernescu Reference Mei and Vernescu2010, Hornung Reference Hornung2012). For this procedure to be valid, it is necessary for the ratio between microscale and macroscale lengthscales to be small. Although traditional homogenization techniques also require the microstructure to be strictly periodic, it can be extended to microstructures that are only near-periodic. Some relevant formal analysis results include Belyaev, Pyatnitskii & Chechkin (Reference Belyaev, Pyatnitskii and Chechkin1998) and Chechkin & Piatnitski (Reference Chechkin and Piatnitski1999), and more applied work using this homogenization method includes van Noorden (Reference van Noorden2009), Peter & Böhm (Reference Peter and Böhm2009), Fatima et al. (Reference Fatima, Arab, Zemskov and Muntean2011), van Noorden & Muntean (Reference van Noorden and Muntean2011), Richardson & Chapman (Reference Richardson and Chapman2011), Ray et al. (Reference Ray, van Noorden, Frank and Knabner2012) and Ray, Elbinger & Knabner (Reference Ray, Elbinger and Knabner2015). Of these, Ray et al. (Reference Ray, van Noorden, Frank and Knabner2012) is particularly relevant to the work in this paper. Here, the authors considered the dynamics of colloids suspended within a fluid moving past circular obstacles in a square array. The colloidal particles were transported by diffusion, advection and an interaction potential between particles, and were able to attach to and detach from the obstacles. The colloidal attachment to the obstacles caused the microstructure to vary from its initially strictly periodic state. The authors used near-periodic homogenization techniques to upscale the problem, deriving averaged equations for the flow problem, the colloidal transport and the microstructure evolution, which were all coupled. The solution scheme for these homogenized equations involved solving a cell problem at each point in time and space in the macroscale domain, and the results focused on the impact of the interaction potential between colloidal particles.
Whilst this extension of standard homogenization theory to near-periodic microstructures is very useful, the added generality of the extension comes with a drawback. A varying microstructure means that the cell problem must be solved at every point in the macroscale, increasing the computational expense required to solve the problem. This issue is bypassed in Bruna & Chapman (Reference Bruna and Chapman2015), by considering a microstructure consisting of an array of spherical obstacles whose radii vary spatially over the macroscale length. Imposing a specific one-parameter shape on the microstructure means that the cell problem is identified by a single parameter, the cell porosity, and the macroscale equation can be written explicitly as a function of the porosity. Moreover, as the cell problems are only dependent on the cell porosity, the coefficients for the macroscale problem can be calculated beforehand, yielding a macroscale problem that can be solved at a significantly reduced computational expense.
In a previous work (Dalwadi et al. Reference Dalwadi, Griffiths and Bruna2015), we considered solute transport through a porosity-graded filter. The filter was modelled as a series of fixed spherical obstacles with a near-periodic spatial variation on the macroscale, through which the contaminated fluid flowed. The techniques developed in Bruna & Chapman (Reference Bruna and Chapman2015) were used to derive a homogenized advection–diffusion–reaction macroscale equation, where the cell problems defining the macroscale coefficients were only dependent on the cell porosity. Although the particles that make up the solute would physically accumulate on the surface of each obstacle, causing the filter geometry to vary over time, we assumed that the particles were sufficiently small and dilute that this effect was negligible on the timescale of interest. Thus, we did not consider the filter evolution or the subsequent filter blocking. We showed that filters whose porosity decreases with depth have a much more uniform adsorption than those with increasing porosity. Moreover, we conjectured that such uniform adsorption would prolong filter lifespan, as contaminant removal would be distributed throughout the filter rather than near the filter entrance, and thus blocking may occur ‘all at once’ rather than in just one place. The aim of this paper is to confirm this conjecture, to quantify the time-dependent evolution of the filter and to determine how to construct a filter that will block all at once.
In this paper, we derive a generalized mathematical model for a porosity-graded filter which explicitly accounts for changes in the filter microstructure over time due to particle adsorption. The time evolution of the filter geometry means that the full fluid and transport problems are coupled. We use a homogenization method that accounts for a microstructure that varies over the macroscale to systematically determine an effective macroscale equation. The filter material is modelled as a lattice of obstacles, whose size varies both spatially over a long lengthscale and temporally, past which a fluid with suspended contaminant particles flows. The particles are transported via advection and diffusion, and can be trapped on the filter as they come into contact with the obstacles. The accumulation of contaminant particles via trapping modifies the filter geometry, and this process eventually leads to pore blockage within the filter. We assume that the contaminant particles are small and in dilute suspension within the fluid. Thus, particle–particle interactions are negligible, and the dominant interaction is the adsorption of contaminant on the obstacles, which we assume is proportional to the concentration of contaminant at the surface. Due to the slow accumulation of contaminant particles, the cell problems we obtain are only dependent on the cell porosity, allowing us to solve the macroscale problem we derive in a computationally efficient manner.
We present the full flow and contaminant transport problems in § 2, and homogenize these to obtain effective equations on the macroscale in § 3. In § 4 we investigate the effect of blocking by considering a filter whose porosity varies in one direction only, with a unidirectional flow in the same direction as the porosity gradient. We use asymptotic methods to significantly reduce the computational expense of solving our system, allowing us to perform an efficient parameter sweep, and we use our results to validate the conjectures made in Dalwadi et al. (Reference Dalwadi, Griffiths and Bruna2015). In addition, these asymptotic results allow us to solve the inverse problem of determining which initial porosity distributions result in filters that block everywhere at once. We conclude in § 5 with an overview and discussion of our results, and we discuss the challenge of determining the filter that globally optimizes contaminant removal.
2 Model description
We consider the transport of contaminant particles via advection and diffusion through a porosity-graded filter, and its resulting blocking. The filter is modelled as a collection of infinitely long solid cylindrical fibres with circular cross-section, whose axes are parallel and where the fluid flows normal to these axes. Contaminants can adsorb onto the solid fibre surface (obstacle surface) and, over time, the build-up of contaminant will result in a reduction in porosity. We describe the contaminant particle distribution in terms of the concentration
$\tilde{c}(\tilde{\boldsymbol{x}},\tilde{t})$
, where
$\tilde{c}$
is the number of moles of contaminant per unit volume,
$\tilde{\boldsymbol{x}}$
is the spatial vector coordinate and
$\tilde{t}$
is time.
Although the problem we have described is three-dimensional, the inherent symmetry allows us to consider a two-dimensional problem as follows. The concentration field is defined within the fluid phase of the domain
$\unicode[STIX]{x1D6FA}_{f}(\tilde{t})\subset \mathbb{R}^{2}$
. We denote the entire domain by
$\unicode[STIX]{x1D6FA}\subset \mathbb{R}^{2}$
, which we refer to as the porous medium. We define the solid phase of the domain as
$\unicode[STIX]{x1D6FA}_{s}\subset \mathbb{R}^{2}$
, noting that
$\unicode[STIX]{x1D6FA}_{s}(\tilde{t})=\unicode[STIX]{x1D6FA}\setminus \unicode[STIX]{x1D6FA}_{f}(\tilde{t})$
, and we define the boundary between fluid and solid phases,
$\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{f}(\tilde{t})$
, which we refer to as the obstacle surface.
The solid phase is modelled by a collection of fixed non-overlapping discs, whose centres are located on a hexagonal lattice at a distance
$\unicode[STIX]{x1D6FF}l$
apart, where
$\unicode[STIX]{x1D6FF}$
is a (small) dimensionless parameter and
$l$
is the characteristic filter depth. Thus,
$\unicode[STIX]{x1D6FF}\ll 1$
is the ratio of fibre centre separation and filter depth. We allow the radii of the circles to vary in both space and time, and a circle with centre at
$\tilde{\boldsymbol{x}}$
has radius
$\tilde{R}(\tilde{t};\tilde{\boldsymbol{x}})$
, where
$2\tilde{R}\leqslant \unicode[STIX]{x1D6FF}l$
and blocking occurs when equality is reached. Thus, blocking does not correspond to the fraction of solid phase reaching
$1$
. We use a hexagonal lattice in contrast to the cubic lattice used by Dalwadi et al. (Reference Dalwadi, Griffiths and Bruna2015) so that we can reach a higher solid fraction at the point of blocking: arranging the obstacles on a cubic lattice, we can only reach a blocking fraction of
$\unicode[STIX]{x03C0}/4\approx 0.785$
, whereas the hexagonal lattice at contact can reach a blocking fraction of
$\unicode[STIX]{x03C0}/(2\sqrt{3})=0.907$
. A schematic of this set-up is shown in figure 1(a).
Figure 1. Model schematic in two dimensions. (a) An example in which the macroscale porosity decreases in the direction of the flow. (b) A magnified view of a given cell
$\unicode[STIX]{x1D714}(\boldsymbol{x},t)$
, with microscale coordinate
$\boldsymbol{y}\in [-1/2,1/2]\times [-\sqrt{3}/2,\sqrt{3}/2]$
.
The pore space is assumed to be entirely saturated by an incompressible Newtonian fluid, which satisfies the Stokes equations,
$$\begin{eqnarray}\displaystyle & \displaystyle -\tilde{\unicode[STIX]{x1D735}}\tilde{p}+\unicode[STIX]{x1D707}\tilde{\unicode[STIX]{x1D6FB}}^{2}\tilde{\boldsymbol{u}}=\mathbf{0},\quad \tilde{\boldsymbol{x}}\in \unicode[STIX]{x1D6FA}_{f}(\tilde{t}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\unicode[STIX]{x1D735}}\boldsymbol{\cdot }\tilde{\boldsymbol{u}}=0,\quad \tilde{\boldsymbol{x}}\in \unicode[STIX]{x1D6FA}_{f}(\tilde{t}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\boldsymbol{u}}=-\frac{\unicode[STIX]{x2202}\tilde{R}}{\unicode[STIX]{x2202}\tilde{t}}\boldsymbol{n},\quad \tilde{\boldsymbol{x}}\in \unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{f}(\tilde{t}), & \displaystyle\end{eqnarray}$$
$\tilde{\unicode[STIX]{x1D735}}$
refers to the nabla operator with respect to
$\tilde{\boldsymbol{x}}$
,
$\unicode[STIX]{x1D707}$
is the constant fluid viscosity,
$\tilde{\boldsymbol{u}}(\tilde{\boldsymbol{x}},\tilde{t})$
is the fluid velocity,
$\tilde{p}(\tilde{\boldsymbol{x}},\tilde{t})$
is the fluid pressure and
$\boldsymbol{n}$
is the unit normal to the obstacle surface directed into the obstacle. Here, (2.1c
) arises from applying a standard no-slip boundary condition to a spherical obstacle, with radius
$\tilde{R}(\tilde{t};\tilde{\boldsymbol{x}})$
, growing radially.
We assume that the contaminant diffuses within the fluid with a constant diffusion coefficient,
$D$
, and is advected by the velocity field. The governing equation is thus
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\tilde{c}}{\unicode[STIX]{x2202}\tilde{t}}=\tilde{\unicode[STIX]{x1D735}}\boldsymbol{\cdot }(D\tilde{\unicode[STIX]{x1D735}}c-\tilde{\boldsymbol{u}}\tilde{c}),\quad \tilde{\boldsymbol{x}}\in \unicode[STIX]{x1D6FA}_{f}(\tilde{t}). & & \displaystyle\end{eqnarray}$$
To determine the correct boundary condition for the concentration on a moving boundary, it is helpful to consider the rate of change of the total number of moles of contaminant in an arbitrary time-dependent volume
$V(t)$
with boundary
$\unicode[STIX]{x2202}V(t)$
moving with the flow:
$$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}\tilde{t}}\int _{V(\tilde{t})}\tilde{c}\,\text{d}V=\int _{V(\tilde{t})}\frac{\unicode[STIX]{x2202}\tilde{c}}{\unicode[STIX]{x2202}\tilde{t}}\,\text{d}V+\int _{\unicode[STIX]{x2202}V(\tilde{t})}(\boldsymbol{n}\boldsymbol{\cdot }\tilde{\boldsymbol{v}})\tilde{c}\,\text{d}S=\int _{\unicode[STIX]{x2202}V(\tilde{t})}\boldsymbol{n}\boldsymbol{\cdot }(D\tilde{\unicode[STIX]{x1D735}}\tilde{c}+(\tilde{\boldsymbol{v}}-\tilde{\boldsymbol{u}})\tilde{c})\,\text{d}S, & & \displaystyle\end{eqnarray}$$
where
$\tilde{\boldsymbol{v}}$
is the velocity of the boundary and
$\boldsymbol{n}$
is the outward unit normal to
$\unicode[STIX]{x2202}V(t)$
. In (2.3), we have used the divergence theorem and (2.2). As
$\tilde{\boldsymbol{u}}=\tilde{\boldsymbol{v}}$
on
$\unicode[STIX]{x2202}V(t)$
, equation (2.3) reduces to
$$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}\tilde{t}}\int _{V(\tilde{t})}\tilde{c}\,\text{d}V=\int _{\unicode[STIX]{x2202}V(\tilde{t})}\boldsymbol{n}\boldsymbol{\cdot }D\tilde{\unicode[STIX]{x1D735}}\tilde{c}\,\text{d}S. & & \displaystyle\end{eqnarray}$$
We suppose that the rate of change of the concentration is linearly dependent on the concentration at the obstacle surface. Thus, using (2.4) in an infinitesimal fluid volume, partly bounded by some of the obstacle surface, we obtain a partially adsorbing Robin boundary condition:
$$\begin{eqnarray}\displaystyle -\unicode[STIX]{x1D6FE}\tilde{c}=\boldsymbol{n}\boldsymbol{\cdot }D\tilde{\unicode[STIX]{x1D735}}\tilde{c},\quad \tilde{\boldsymbol{x}}\in \unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{f}(\tilde{t}), & & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D6FE}\geqslant 0$
is the constant contaminant-adsorption coefficient. There is no adsorption when
$\unicode[STIX]{x1D6FE}=0$
, and the adsorption is instantaneous in the limit as
$\unicode[STIX]{x1D6FE}\rightarrow \infty$
.
Finally, we must couple the growth of the obstacles to the accumulation of contaminant at the obstacle boundary. This is in the form of a volume conservation law, where the volume of contaminant lost to adhesion is equal to the obstacle volume gain. If we multiply (2.3) by
$V_{m}$
, the molar volume of the contaminant in the filtrate (defined as the volume of contaminant per mole of contaminant), we now have an equation for the rate of change of the volume of contaminant in an arbitrary volume of fluid. Applying this to an infinitesimal volume, partially bounded by part of the obstacle surface, and equating the loss of contaminant at the surface to the gain of obstacle volume, we obtain
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\tilde{R}}{\unicode[STIX]{x2202}\tilde{t}}=-V_{m}\boldsymbol{n}\boldsymbol{\cdot }D\tilde{\unicode[STIX]{x1D735}}\tilde{c},\quad \tilde{\boldsymbol{x}}\in \unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{f}(\tilde{t}). & & \displaystyle\end{eqnarray}$$
We have included the effect of filter porosity changing due to contaminant particle adsorption but neglected the effect of particle–particle interactions because the contaminant particles are small compared with the obstacles and in dilute suspension within the fluid, hence particle–particle interactions occur far less frequently than particle–obstacle interactions. The interfacial condition (2.6) is similar to a Stefan condition for phase-change problems in heat transport (Gupta Reference Gupta2003).
2.1 Dimensionless equations
We scale the variables via
$\tilde{\boldsymbol{x}}=l\boldsymbol{x}$
,
$\tilde{\boldsymbol{u}}={\mathcal{U}}\boldsymbol{u}$
,
$\tilde{t}=(\unicode[STIX]{x1D6FF}l/(c_{\infty }V_{m}\unicode[STIX]{x1D6FE}))t$
,
$\tilde{R}=\unicode[STIX]{x1D6FF}lR$
,
$\tilde{c}=c_{\infty }c$
and
$\tilde{p}=(\unicode[STIX]{x1D707}{\mathcal{U}}/(\unicode[STIX]{x1D6FF}^{2}l))p$
, where
${\mathcal{U}}$
and
$c_{\infty }$
are characteristic velocity and concentration scales respectively. The time scale is chosen to balance obstacle growth with contaminant concentration in (2.6), and the pressure scale is chosen to balance the pressure gradient over the macroscale with viscous forces over the obstacle microscale. The flow problem (2.1) then transforms to
$$\begin{eqnarray}\displaystyle & \displaystyle -\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D6FF}^{2}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}=\mathbf{0},\quad \boldsymbol{x}\in \unicode[STIX]{x1D6FA}_{f}(t), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0,\quad \boldsymbol{x}\in \unicode[STIX]{x1D6FA}_{f}(t), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}=-\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6FC}\frac{\unicode[STIX]{x2202}R}{\unicode[STIX]{x2202}t}\boldsymbol{n},\quad \boldsymbol{x}\in \unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{f}(t), & \displaystyle\end{eqnarray}$$
$\unicode[STIX]{x1D6FC}=c_{\infty }V_{m}\unicode[STIX]{x1D6FE}/(\unicode[STIX]{x1D6FF}{\mathcal{U}})$
.
For the concentration problem given by (2.2), (2.5) and (2.6), we obtain the dimensionless solute-transport equation
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FC}\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}t}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }(Pe^{-1}\unicode[STIX]{x1D735}c-\boldsymbol{u}c),\quad \boldsymbol{x}\in \unicode[STIX]{x1D6FA}_{f}(t), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle -\unicode[STIX]{x1D6FF}kc=\boldsymbol{n}\boldsymbol{\cdot }Pe^{-1}\unicode[STIX]{x1D735}c,\quad \boldsymbol{x}\in \unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{f}(t), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}R}{\unicode[STIX]{x2202}t}=c,\quad \boldsymbol{x}\in \unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{f}(t), & \displaystyle\end{eqnarray}$$
$Pe={\mathcal{U}}l/D$
is the Péclet number and
$k=\unicode[STIX]{x1D6FE}/(\unicode[STIX]{x1D6FF}{\mathcal{U}})$
. The boundary condition (2.8b
) expresses the flux across the moving boundary relative to the boundary, and as such does not contain any velocity component. This system reduces to that considered in Dalwadi et al. (Reference Dalwadi, Griffiths and Bruna2015) when
$R$
is independent of
$t$
, i.e. when
$V_{m}=0$
, which can be seen by scaling
$t\sim \unicode[STIX]{x1D6FC}$
and then taking
$\unicode[STIX]{x1D6FC}\rightarrow 0$
. We assume that
$\unicode[STIX]{x1D6FC}$
,
$Pe$
and
$k$
are all
$O(1)$
parameters, which corresponds to the richest asymptotic limit. That is, all mechanisms contribute at leading order, from which all asymptotic sublimits may be distilled. In practice, the contaminant accumulates on the obstacles over a much longer timescale than the fluid takes to travel through the porous medium. This manifests in the smallness of the parameter
$\unicode[STIX]{x1D6FC}$
, and we exploit this feature in § 4.2 to consider the physically relevant effect of slow contaminant accumulation.
In dimensionless units, the obstacles now form a two-dimensional hexagonal lattice of circles whose centres are a distance of
$\unicode[STIX]{x1D6FF}$
apart, and a circle with centre at
$\boldsymbol{x}$
has radius
$\unicode[STIX]{x1D6FF}R(\boldsymbol{x},t)$
(figure 1).
3 Homogenization
The complexity of the problem geometry is reduced by homogenizing the governing equations (2.7)–(2.8) using the method of multiple scales. This provides effective equations on a simpler macroscale domain, which formally capture the relevant information about the microscale geometry (see, for example, Sánchez-Palencia (Reference Sánchez-Palencia1980) for a standard derivation of Darcy’s law, Conca (Reference Conca1985) for a homogenization of the Navier–Stokes equations for various boundary conditions and Auriault & Adler (Reference Auriault and Adler1995) for a derivation of an advection–diffusion solute transport model in porous media). Following standard homogenization theory, we introduce a microscale variable
$\boldsymbol{y}=\boldsymbol{x}/\unicode[STIX]{x1D6FF}$
and treat
$\boldsymbol{x}$
and
$\boldsymbol{y}$
as independent. This adds an extra degree of freedom which is later removed by imposing that the solution is periodic in
$\boldsymbol{y}$
. Hence, any small variation between unit cells is captured through the macroscale variable
$\boldsymbol{x}$
. As shown in figure 1(b), the microscale variable
$\boldsymbol{y}$
is defined in the unit cell
$\unicode[STIX]{x1D714}(\boldsymbol{x},t)$
, whereas the macroscale variable
$\boldsymbol{x}$
spans across the whole filter. The solid portion of the cell, occupied by obstacles, is denoted by
$\unicode[STIX]{x1D714}_{s}(\boldsymbol{x},t)$
. The fluid portion of the cell is
$\unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t)=\unicode[STIX]{x1D714}(\boldsymbol{x},t)\backslash \unicode[STIX]{x1D714}_{s}(\boldsymbol{x},t)$
. The boundary between solid and fluid portions within the cell is denoted by
$\unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t)$
, while the outer boundary of the unit cell is
$\unicode[STIX]{x2202}\unicode[STIX]{x1D714}(\boldsymbol{x},t)$
(see figure 1). Further, we treat each dependent variable as a function of both
$\boldsymbol{x}$
and
$\boldsymbol{y}$
. Thus, spatial derivatives transform in the following manner
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}\mapsto \unicode[STIX]{x1D735}_{\boldsymbol{x}}+{\displaystyle \frac{1}{\unicode[STIX]{x1D6FF}}}\unicode[STIX]{x1D735}_{\boldsymbol{y}}, & & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D735}_{\boldsymbol{x}}$
and
$\unicode[STIX]{x1D735}_{\boldsymbol{y}}$
refer to the nabla operator in the
$\boldsymbol{x}$
- and
$\boldsymbol{y}$
-coordinate systems, respectively. The introduction of this new spatial variable also changes the normal vector
$\boldsymbol{n}$
, used to evaluate the boundary conditions (2.7c
) and (2.8b
), to
$$\begin{eqnarray}\displaystyle \boldsymbol{n}={\displaystyle \frac{\boldsymbol{n}_{\boldsymbol{y}}+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D735}_{\boldsymbol{x}}R}{\Vert \boldsymbol{n}_{\boldsymbol{y}}+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D735}_{\boldsymbol{x}}R\Vert }}, & & \displaystyle\end{eqnarray}$$
where
$\boldsymbol{n}_{\boldsymbol{y}}=-\unicode[STIX]{x1D735}_{\boldsymbol{y}}\Vert \boldsymbol{y}\Vert =-\boldsymbol{y}/\Vert \boldsymbol{y}\Vert$
is the geometric outward unit normal on the obstacle boundary
$\unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t)$
, and
$\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D735}_{\boldsymbol{x}}R$
accounts for the macroscale effect of varying obstacle size. The details of this transformation are described in Dalwadi et al. (Reference Dalwadi, Griffiths and Bruna2015).
As our goal is to derive effective governing equations that are valid in the macroscale domain, we consider variables averaged over an entire cell
$\unicode[STIX]{x1D714}(\boldsymbol{x},t)$
. To this end, we define the macroscale porosity
$\unicode[STIX]{x1D719}(\boldsymbol{x},t)$
to be
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}(\boldsymbol{x},t)={\displaystyle \frac{|\unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t)|}{|\unicode[STIX]{x1D714}(\boldsymbol{x},t)|}}, & & \displaystyle\end{eqnarray}$$
and the volumetric average concentration
$C$
and volumetric average fluid velocity
$\boldsymbol{U}$
(known as the Darcy velocity in porous-media formulations (Kaviany Reference Kaviany2012)) as follows
$$\begin{eqnarray}\displaystyle & \displaystyle C(\boldsymbol{x},t)={\displaystyle \frac{1}{|\unicode[STIX]{x1D714}(\boldsymbol{x},t)|}}\int _{\unicode[STIX]{x1D714}(\boldsymbol{x},t)}c(\boldsymbol{x},\boldsymbol{y},t)\,\text{d}\boldsymbol{y}={\displaystyle \frac{1}{|\unicode[STIX]{x1D714}(\boldsymbol{x},t)|}}\int _{\unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t)}c(\boldsymbol{x},\boldsymbol{y},t)\,\text{d}\boldsymbol{y}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{U}(\boldsymbol{x},t)={\displaystyle \frac{1}{|\unicode[STIX]{x1D714}(\boldsymbol{x},t)|}}\int _{\unicode[STIX]{x1D714}(\boldsymbol{x},t)}\boldsymbol{u}(\boldsymbol{x},\boldsymbol{y},t)\,\text{d}\boldsymbol{y}={\displaystyle \frac{1}{|\unicode[STIX]{x1D714}(\boldsymbol{x},t)|}}\int _{\unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t)}\boldsymbol{u}(\boldsymbol{x},\boldsymbol{y},t)\,\text{d}\boldsymbol{y}, & \displaystyle\end{eqnarray}$$
$c=0$
and
$\boldsymbol{u}\equiv \mathbf{0}$
in
$\unicode[STIX]{x1D714}_{s}(\boldsymbol{x},t)$
.
The homogenization in this section proceeds in a similar manner to other homogenization work, such as van Noorden (Reference van Noorden2009), Ray et al. (Reference Ray, van Noorden, Frank and Knabner2012, Reference Ray, Elbinger and Knabner2015), Bruna & Chapman (Reference Bruna and Chapman2015) and Dalwadi et al. (Reference Dalwadi, Griffiths and Bruna2015). We first consider the flow problem (2.7), and use the results in the solute-transport problem (2.8).
3.1 Flow problem
Under the spatial transforms (3.1) and (3.2), the flow equations (2.7) become
$$\begin{eqnarray}\displaystyle & \displaystyle -(\unicode[STIX]{x1D6FF}^{-1}\unicode[STIX]{x1D735}_{\boldsymbol{y}}+\unicode[STIX]{x1D735}_{\boldsymbol{x}})p+(\unicode[STIX]{x1D735}_{\boldsymbol{y}}+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D735}_{\boldsymbol{x}})^{2}\boldsymbol{u}=\mathbf{0},\quad \boldsymbol{y}\in \unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle (\unicode[STIX]{x1D735}_{\boldsymbol{y}}+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D735}_{\boldsymbol{x}})\boldsymbol{\cdot }\boldsymbol{u}=0,\quad \boldsymbol{y}\in \unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}=-\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6FC}\frac{\unicode[STIX]{x2202}R}{\unicode[STIX]{x2202}t}\boldsymbol{n}_{\boldsymbol{y}}+O(\unicode[STIX]{x1D6FF}^{2}),\quad \boldsymbol{y}\in \unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u},p\text{ periodic},\quad \boldsymbol{y}\in \unicode[STIX]{x2202}\unicode[STIX]{x1D714}(\boldsymbol{x},t). & \displaystyle\end{eqnarray}$$
$\unicode[STIX]{x1D6FF}$
as usual in the multiple-scales method (see Dalwadi et al.
Reference Dalwadi, Griffiths and Bruna2015 for details), we find that the leading-order pressure,
$p_{0}$
, is independent of the microscale, i.e.
$p_{0}=p_{0}(\boldsymbol{x},t)$
. Considering the next order in (3.5) allows us to write 
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}_{0}=-\unicode[STIX]{x1D646}(\boldsymbol{x},\boldsymbol{y},t)\unicode[STIX]{x1D735}_{\boldsymbol{x}}p_{0}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle p_{1}=-\unicode[STIX]{x1D72B}(\boldsymbol{x},\boldsymbol{y},t)\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\boldsymbol{x}}p_{0}+\overline{p}(\boldsymbol{x},t), & \displaystyle\end{eqnarray}$$
$\overline{p}$
is an arbitrary function of the macroscale and time only, and the matrix function
$\unicode[STIX]{x1D646}$
and the vector function
$\unicode[STIX]{x1D72B}$
satisfy the so-called cell problem: 
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D644}-\unicode[STIX]{x1D735}_{\boldsymbol{y}}\unicode[STIX]{x1D72B}+\unicode[STIX]{x1D6FB}_{\boldsymbol{y}}^{2}\unicode[STIX]{x1D646}=\mathsf{0},\quad \boldsymbol{y}\in \unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}_{\boldsymbol{y}}\boldsymbol{\cdot }\unicode[STIX]{x1D646}=\mathbf{0},\quad \boldsymbol{y}\in \unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D646}=\mathsf{0},\quad \boldsymbol{y}\in \unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D646},\unicode[STIX]{x1D72B}\text{ periodic},\quad \boldsymbol{y}\in \unicode[STIX]{x2202}\unicode[STIX]{x1D714}(\boldsymbol{x},t). & \displaystyle\end{eqnarray}$$
$\unicode[STIX]{x1D644}$
is the two-dimensional identity matrix and the dependence of
$\unicode[STIX]{x1D646}$
and
$\unicode[STIX]{x1D72B}$
on
$\boldsymbol{x}$
and
$t$
is due to the dependence of the microscale boundary
$\unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t)$
on the macroscale variable and time.
Integrating (3.6a
) over
$\unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t)$
and using (3.4b
), we obtain the homogenized Darcy relation
$$\begin{eqnarray}\displaystyle \boldsymbol{U}(\boldsymbol{x},t)=-{\mathcal{K}}(\unicode[STIX]{x1D719})\unicode[STIX]{x1D735}_{\boldsymbol{x}}p & & \displaystyle\end{eqnarray}$$
at leading order, where
${\mathcal{K}}(\unicode[STIX]{x1D719})$
is a scalar function that contains the necessary microscale structure information, and is defined by
$$\begin{eqnarray}\displaystyle {\mathcal{K}}(\unicode[STIX]{x1D719})\unicode[STIX]{x1D644}={\displaystyle \frac{1}{|\unicode[STIX]{x1D714}(\boldsymbol{x},t)|}}\int _{\unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t)}\unicode[STIX]{x1D646}\,\text{d}\boldsymbol{y}. & & \displaystyle\end{eqnarray}$$
We note that the integral of
$\unicode[STIX]{x1D646}$
in (3.8b
) is a multiple of the identity matrix due to the symmetry of the cell problem described by (3.7). This would not be the case if, instead of circles, we had considered obstacles with fewer than two orthogonal planes of symmetry. Finally, to obtain a closed homogenized flow problem, we consider the
$O(\unicode[STIX]{x1D6FF})$
terms in (3.5b–d
), given by
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}_{\boldsymbol{y}}\boldsymbol{\cdot }\boldsymbol{u}_{1}+\unicode[STIX]{x1D735}_{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{u}_{0}=0,\quad \boldsymbol{y}\in \unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}_{1}=-\unicode[STIX]{x1D6FC}\frac{\unicode[STIX]{x2202}R}{\unicode[STIX]{x2202}t}\boldsymbol{n}_{\boldsymbol{y}},\quad \boldsymbol{y}\in \unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}_{1}~\text{periodic},\quad \boldsymbol{y}\in \unicode[STIX]{x2202}\unicode[STIX]{x1D714}(\boldsymbol{x},t). & \displaystyle\end{eqnarray}$$
$\unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t)$
and using Reynolds transport theorem in conjunction with the periodic and no-slip boundary conditions (3.9b,c
) yields a continuity equation for the macroscale fluid velocity 
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}_{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{U}=\unicode[STIX]{x1D6FC}{\displaystyle \frac{|\unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t)|}{|\unicode[STIX]{x1D714}(\boldsymbol{x},t)|}}\frac{\unicode[STIX]{x2202}R}{\unicode[STIX]{x2202}t}. & & \displaystyle\end{eqnarray}$$
The system (3.8) and (3.10) determines the flow problem given the obstacle structure, which manifests itself through the source term on the right-hand side of (3.10). This reflects the fact that reducing the available pore space within a cell for an incompressible fluid will push the fluid out of that cell. If
$\unicode[STIX]{x2202}R/\unicode[STIX]{x2202}t\equiv 0$
, we recover the system derived in Dalwadi et al. (Reference Dalwadi, Griffiths and Bruna2015).
3.2 Solute-transport problem
Under the spatial transforms (3.1) and (3.2), the solute-transport problem (2.8) becomes
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FF}^{2}\unicode[STIX]{x1D6FC}\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}t}=(\unicode[STIX]{x1D735}_{\boldsymbol{y}}+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D735}_{\boldsymbol{x}})\boldsymbol{\cdot }(Pe^{-1}(\unicode[STIX]{x1D735}_{\boldsymbol{ y}}+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D735}_{\boldsymbol{x}})c-\unicode[STIX]{x1D6FF}\boldsymbol{u}c),\quad \boldsymbol{y}\in \unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle -\unicode[STIX]{x1D6FF}^{2}kc=(\boldsymbol{n}_{\boldsymbol{y}}+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D735}_{\boldsymbol{x}}R)\boldsymbol{\cdot }(Pe^{-1}(\unicode[STIX]{x1D735}_{\boldsymbol{y}}+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D735}_{\boldsymbol{x}})c)+O(\unicode[STIX]{x1D6FF}^{3}),\quad \boldsymbol{y}\in \unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle c\text{ periodic},\quad \boldsymbol{y}\in \unicode[STIX]{x2202}\unicode[STIX]{x1D714}(\boldsymbol{x},t). & \displaystyle\end{eqnarray}$$
$c(\boldsymbol{x},\boldsymbol{y},t)$
in powers of
$\unicode[STIX]{x1D6FF}$
, we find that the leading-order concentration,
$c_{0}$
, is independent of the microscale, i.e.
$c_{0}=c_{0}(\boldsymbol{x},t)$
. The
$O(\unicode[STIX]{x1D6FF})$
terms in (3.11) allow us to write
$c_{1}$
as 
$$\begin{eqnarray}\displaystyle c_{1}(\boldsymbol{x},\boldsymbol{y},t)=-\unicode[STIX]{x1D71E}(\boldsymbol{x},\boldsymbol{y},t)\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\boldsymbol{x}}c_{0}(\boldsymbol{x},t)+\bar{c}(\boldsymbol{x},t), & & \displaystyle\end{eqnarray}$$
where
$\bar{c}$
is an arbitrary function of the macroscale and time only, and the components of the function
$\unicode[STIX]{x1D71E}$
satisfy the cell problem
$$\begin{eqnarray}\displaystyle & \displaystyle 0=\unicode[STIX]{x1D6FB}_{\boldsymbol{y}}^{2}\unicode[STIX]{x1D6E4}_{i},\quad \boldsymbol{y}\in \unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle n_{y,i}=\boldsymbol{n}_{\boldsymbol{y}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\boldsymbol{y}}\unicode[STIX]{x1D6E4}_{i},\quad \boldsymbol{y}\in \unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6E4}_{i}~\text{periodic},\quad \boldsymbol{y}\in \unicode[STIX]{x2202}\unicode[STIX]{x1D714}(\boldsymbol{x},t), & \displaystyle\end{eqnarray}$$
$n_{y,i}$
is the
$i$
th component of the unit vector
$\boldsymbol{n}_{\boldsymbol{y}}$
.
The effect of the obstacle growth only appears at the next order. Namely, the
$O(\unicode[STIX]{x1D6FF}^{2})$
terms in (3.11) are
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FC}\frac{\unicode[STIX]{x2202}c_{0}}{\unicode[STIX]{x2202}t} & = & \displaystyle \unicode[STIX]{x1D735}_{\boldsymbol{y}}\boldsymbol{\cdot }(Pe^{-1}(\unicode[STIX]{x1D735}_{\boldsymbol{ y}}c_{2}+\unicode[STIX]{x1D735}_{\boldsymbol{x}}c_{1})-\boldsymbol{u}_{1}c_{0}-\boldsymbol{u}_{0}c_{1})\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D735}_{\boldsymbol{x}}\boldsymbol{\cdot }(Pe^{-1}(\unicode[STIX]{x1D735}_{\boldsymbol{y}}c_{1}+\unicode[STIX]{x1D735}_{\boldsymbol{x}}c_{0})-\boldsymbol{u}_{0}c_{0}),\quad \boldsymbol{y}\in \unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle -kc_{0} & = & \displaystyle \boldsymbol{n}_{\boldsymbol{y}}\boldsymbol{\cdot }(Pe^{-1}(\unicode[STIX]{x1D735}_{\boldsymbol{y}}c_{2}+\unicode[STIX]{x1D735}_{\boldsymbol{x}}c_{1}))\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D735}_{\boldsymbol{x}}R\boldsymbol{\cdot }(Pe^{-1}(\unicode[STIX]{x1D735}_{\boldsymbol{y}}c_{1}+\unicode[STIX]{x1D735}_{\boldsymbol{x}}c_{0})),\quad \boldsymbol{y}\in \unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle c_{2}\text{ periodic},\quad \boldsymbol{y}\in \unicode[STIX]{x2202}\unicode[STIX]{x1D714}(\boldsymbol{x},t). & & \displaystyle\end{eqnarray}$$
Integrating (3.14a
) over
$\unicode[STIX]{x1D714}_{f}$
and applying the boundary conditions (3.5c
), (3.9b
), and (3.14b–c
) gives
$$\begin{eqnarray}\displaystyle \int _{\unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t)}\unicode[STIX]{x1D6FC}\frac{\unicode[STIX]{x2202}c_{0}}{\unicode[STIX]{x2202}t}\,\text{d}\boldsymbol{y} & = & \displaystyle \int _{\unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t)}\unicode[STIX]{x1D735}_{\boldsymbol{x}}\boldsymbol{\cdot }(Pe^{-1}(\unicode[STIX]{x1D735}_{\boldsymbol{y}}c_{1}+\unicode[STIX]{x1D735}_{\boldsymbol{x}}c_{0})-\boldsymbol{u}_{0}c_{0})\,\text{d}\boldsymbol{y}\nonumber\\ \displaystyle & & \displaystyle -\,\int _{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t)}\unicode[STIX]{x1D735}_{\boldsymbol{x}}R\boldsymbol{\cdot }(Pe^{-1}(\unicode[STIX]{x1D735}_{\boldsymbol{y}}c_{1}+\unicode[STIX]{x1D735}_{\boldsymbol{x}}c_{0})-\boldsymbol{u}_{0}c_{0})\,\text{d}s\nonumber\\ \displaystyle & & \displaystyle -\,\int _{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t)}kc_{0}\,\text{d}s+\int _{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t)}\unicode[STIX]{x1D6FC}\frac{\unicode[STIX]{x2202}R}{\unicode[STIX]{x2202}t}c_{0}\,\text{d}s,\end{eqnarray}$$
where
$\text{d}s$
denotes the differential element of the obstacle surface
$\unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t)$
. Using Reynolds transport theorem, equation (3.15) reduces to
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FC}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(|\unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t)|c_{0})=\unicode[STIX]{x1D735}_{\boldsymbol{x}}\boldsymbol{\cdot }\int _{\unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t)}(Pe^{-1}(\unicode[STIX]{x1D735}_{\boldsymbol{y}}c_{1}+\unicode[STIX]{x1D735}_{\boldsymbol{x}}c_{0})-\boldsymbol{u}_{0}c_{0})\,\text{d}\boldsymbol{y}-\int _{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t)}kc_{0}\,\text{d}s,\quad & & \displaystyle\end{eqnarray}$$
noting that
$c_{0}$
is independent of
$\boldsymbol{y}$
. From now on, we simply write
$\unicode[STIX]{x1D719}$
but it should be understood that the porosity will be, in general, a function of
$\boldsymbol{x}$
and
$t$
. In a similar manner, henceforth we use
$\unicode[STIX]{x1D714}_{f}(\unicode[STIX]{x1D719})$
instead of
$\unicode[STIX]{x1D714}_{f}(\boldsymbol{x},t)$
and likewise for
$\unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{f}$
. Finally, using (3.12), (3.16) reduces to
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FC}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(|\unicode[STIX]{x1D714}_{f}(\unicode[STIX]{x1D719})|c_{0}) & = & \displaystyle \unicode[STIX]{x1D735}_{\boldsymbol{x}}\boldsymbol{\cdot }\left(Pe^{-1}\left(\int _{\unicode[STIX]{x1D714}_{f}(\unicode[STIX]{x1D719})}(\unicode[STIX]{x1D644}-\unicode[STIX]{x1D645}_{\unicode[STIX]{x1D71E}}^{T})\,\text{d}\boldsymbol{y}\right)\unicode[STIX]{x1D735}_{\boldsymbol{ x}}c_{0}-\left(\int _{\unicode[STIX]{x1D714}_{f}(\unicode[STIX]{x1D719})}\boldsymbol{u}_{0}\,\text{d}\boldsymbol{y}\right)c_{0}\right)\nonumber\\ \displaystyle & & \displaystyle -\,|\unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{f}(\unicode[STIX]{x1D719})|kc_{0},\end{eqnarray}$$
where
$(\unicode[STIX]{x1D645}_{\unicode[STIX]{x1D71E}}^{T})_{ij}=\unicode[STIX]{x2202}\unicode[STIX]{x1D6E4}_{j}/\unicode[STIX]{x2202}y_{i}$
is the transpose of the Jacobian matrix of
$\unicode[STIX]{x1D71E}$
, the solution to the cell problem (3.13).
To express (3.17) in terms of the volumetric average concentration
$C(\boldsymbol{x},t)$
defined in (3.4a
), we note that
$C(\boldsymbol{x},t)\sim |\unicode[STIX]{x1D714}|^{-1}\int _{\unicode[STIX]{x1D714}_{f}(\unicode[STIX]{x1D719})}c_{0}\,\text{d}\boldsymbol{y}=\unicode[STIX]{x1D719}c_{0}(\boldsymbol{x},t)$
at leading order in
$\unicode[STIX]{x1D6FF}$
. Using this relation and (3.8a
), we find that
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FC}\frac{\unicode[STIX]{x2202}C}{\unicode[STIX]{x2202}t}=\unicode[STIX]{x1D735}_{\boldsymbol{x}}\boldsymbol{\cdot }\left({\mathcal{D}}(\unicode[STIX]{x1D719})\unicode[STIX]{x1D735}_{\boldsymbol{x}}C-{\displaystyle \frac{C}{\unicode[STIX]{x1D719}}}(\boldsymbol{U}(\unicode[STIX]{x1D719})+{\mathcal{D}}(\unicode[STIX]{x1D719})\unicode[STIX]{x1D735}_{\boldsymbol{x}}\unicode[STIX]{x1D719})\right)-f(\unicode[STIX]{x1D719})C, & & \displaystyle\end{eqnarray}$$
at leading order in
$\unicode[STIX]{x1D6FF}$
, where the effective diffusion coefficient is
$$\begin{eqnarray}\displaystyle {\mathcal{D}}(\unicode[STIX]{x1D719})\unicode[STIX]{x1D644}=Pe^{-1}\left(\unicode[STIX]{x1D644}-{\displaystyle \frac{1}{|\unicode[STIX]{x1D714}_{f}(\unicode[STIX]{x1D719})|}}\int _{\unicode[STIX]{x1D714}_{f}(\unicode[STIX]{x1D719})}\unicode[STIX]{x1D645}_{\unicode[STIX]{x1D71E}}^{T}\,\text{d}\boldsymbol{y}\right), & & \displaystyle\end{eqnarray}$$
and the effective adsorption coefficient is
$$\begin{eqnarray}\displaystyle f(\unicode[STIX]{x1D719})=k{\displaystyle \frac{|\unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{f}(\unicode[STIX]{x1D719})|}{|\unicode[STIX]{x1D714}_{f}(\unicode[STIX]{x1D719})|}}={\displaystyle \frac{4k\unicode[STIX]{x03C0}R}{\unicode[STIX]{x1D719}\sqrt{3}}}, & & \displaystyle\end{eqnarray}$$
using the fact that
$\unicode[STIX]{x1D719}=1-2\unicode[STIX]{x03C0}R^{2}/\sqrt{3}$
, and
$|\unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{f}(\unicode[STIX]{x1D719})|=4\unicode[STIX]{x03C0}R$
. As is the case for the permeability
${\mathcal{K}}$
, we find that the diffusion tensor is a multiple of the identity due to the symmetry of the cell problem (3.13). Finally, the condition to couple contaminant concentration with obstacle growth (2.8c
) becomes
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}\frac{\unicode[STIX]{x2202}R}{\unicode[STIX]{x2202}t}=C. & & \displaystyle\end{eqnarray}$$
The effective permeability
${\mathcal{K}}(\unicode[STIX]{x1D719})$
(from (3.8b
)) and the effective diffusion
${\mathcal{D}}(\unicode[STIX]{x1D719})$
(from (3.18b
)) that encapsulate the microscale physics can be computed by solving the respective cell problems (3.7) and (3.13) for a given cell porosity
$\unicode[STIX]{x1D719}$
, determined by the size of the obstacles. We calculate this numerically using the finite-element software Comsol Multiphysics and find that
${\mathcal{K}}$
and
${\mathcal{D}}$
are both monotonically increasing with the porosity, as expected, and that there is a sharp decrease in both functions as the porosity tends down towards the critical porosity where blocking occurs (see figure 2).
The homogenization procedure has significantly reduced the mathematical complexity of the growing multiply-connected domain in the full problem (2.8), with only a marginal increase in the complexity of the coefficients in the resulting governing equations (3.18), which capture the microscale structure of the problem in a systematic manner.
4 A unidirectionally graded filter
4.1 Model set-up
We are now in a position to use the homogenized equations (3.8), (3.10) and (3.18) to quantify the effect of blocking and porosity gradients on filter efficiency. We consider a canonical industrial process of interest whereby a three-dimensional filter separates two reservoirs, and a unidirectional flow is induced through the filter. As is common in industrial set-ups, we assume that the filter is graded in the same direction as the fluid flow. This reduces the problem to a one-dimensional description.
We define the direction of porosity variation as
$x$
, where
$x\in (0,1)$
within the filter (since the dimensional characteristic length
$l$
is chosen to be the depth of the filter). Thus, the set-up is similar to that illustrated in figure 1. The upstream is defined for
$x\in (-\infty ,0)$
and the downstream for
$x\in (1,\infty )$
.
Provided the boundary conditions allow for unidirectionality, equations (3.8) and (3.10) yield a unidirectional macroscale flux
$\boldsymbol{U}=U(x,t)\boldsymbol{e}_{x}$
(where
$\boldsymbol{e}_{x}$
is the unit vector in the
$x$
-direction). When filtering at constant (dimensionless) pressure drop
$\unicode[STIX]{x0394}P$
, we find that
$$\begin{eqnarray}\displaystyle U(x,t) & = & \displaystyle \overline{K}(t)\left(\unicode[STIX]{x0394}P-\unicode[STIX]{x1D6FC}\int _{0}^{1}\frac{1}{{\mathcal{K}}(\unicode[STIX]{x1D719}(v,t))}\left(\int _{0}^{v}|\unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{f}(u,t)|\frac{\unicode[STIX]{x2202}R}{\unicode[STIX]{x2202}t}(u,t)\,\text{d}u\right)\,\text{d}v\right)\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D6FC}\int _{0}^{x}|\unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{f}(\unicode[STIX]{x1D709},t)|\frac{\unicode[STIX]{x2202}R}{\unicode[STIX]{x2202}t}(\unicode[STIX]{x1D709},t)\,\text{d}\unicode[STIX]{x1D709},\end{eqnarray}$$
where
$$\begin{eqnarray}\displaystyle \overline{K}(t)=\left(\int _{0}^{1}{\displaystyle \frac{\text{d}x}{{\mathcal{K}}(x,t)}}\right)^{-1}. & & \displaystyle\end{eqnarray}$$
Thus, the governing equation (3.18a ) becomes
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FC}\frac{\unicode[STIX]{x2202}C}{\unicode[STIX]{x2202}t}=\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left[{\mathcal{D}}(\unicode[STIX]{x1D719})\frac{\unicode[STIX]{x2202}C}{\unicode[STIX]{x2202}x}-{\displaystyle \frac{C}{\unicode[STIX]{x1D719}}}\left(U(x,t)+{\mathcal{D}}(\unicode[STIX]{x1D719})\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}x}\right)\right]-f(\unicode[STIX]{x1D719})C,\quad x\in (0,1), & & \displaystyle\end{eqnarray}$$
and the relationship between obstacle growth and concentration (3.18d ) is
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}\frac{\unicode[STIX]{x2202}R}{\unicode[STIX]{x2202}t}=C,\quad x\in (0,1). & & \displaystyle\end{eqnarray}$$
Considering the limit of (3.18b
) and (4.3) as
$\unicode[STIX]{x1D719}\rightarrow 1$
, and using continuity of fluid flux, we obtain the following governing equations for the upstream
$C^{-}$
and downstream
$C^{+}$
concentrations
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FC}\frac{\unicode[STIX]{x2202}C^{-}}{\unicode[STIX]{x2202}t}=\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left(Pe^{-1}\frac{\unicode[STIX]{x2202}C^{-}}{\unicode[STIX]{x2202}x}-U(0,t)C^{-}\right),\quad x\in (-\infty ,0), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FC}\frac{\unicode[STIX]{x2202}C^{+}}{\unicode[STIX]{x2202}t}=\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left(Pe^{-1}\frac{\unicode[STIX]{x2202}C^{+}}{\unicode[STIX]{x2202}x}-U(1,t)C^{+}\right),\quad x\in (1,\infty ). & \displaystyle\end{eqnarray}$$
We impose continuity of concentration and concentration flux at the boundaries between the filter and the reservoirs. In the far-field of the reservoirs, the concentration tends to a constant value. We may take the upstream concentration
$C^{-}\rightarrow 1$
(by choice of our non-dimensionalization), while the downstream concentration
$C^{+}$
tends to a constant value that must be determined as part of the solution. Mathematically, this corresponds to
$$\begin{eqnarray}\displaystyle & \displaystyle C^{-}\rightarrow 1,\quad x\rightarrow -\infty , & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle C^{-}={\displaystyle \frac{C}{\unicode[STIX]{x1D719}}},\quad x=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle Pe^{-1}\frac{\unicode[STIX]{x2202}C^{-}}{\unicode[STIX]{x2202}x}-UC^{-}={\mathcal{D}}(\unicode[STIX]{x1D719})\frac{\unicode[STIX]{x2202}C}{\unicode[STIX]{x2202}x}-{\displaystyle \frac{C}{\unicode[STIX]{x1D719}}}\left(U+{\mathcal{D}}(\unicode[STIX]{x1D719})\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}x}\right),\quad x=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle C^{+}={\displaystyle \frac{C}{\unicode[STIX]{x1D719}}},\quad x=1, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle Pe^{-1}\frac{\unicode[STIX]{x2202}C^{+}}{\unicode[STIX]{x2202}x}-UC^{+}={\mathcal{D}}(\unicode[STIX]{x1D719})\frac{\unicode[STIX]{x2202}C}{\unicode[STIX]{x2202}x}-{\displaystyle \frac{C}{\unicode[STIX]{x1D719}}}\left(U+{\mathcal{D}}(\unicode[STIX]{x1D719})\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}x}\right),\quad x=1, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}C^{+}}{\unicode[STIX]{x2202}x}\rightarrow 0,\quad x\rightarrow \infty . & \displaystyle\end{eqnarray}$$
A suitable measure of filter efficiency is the cumulative contaminant removal, defined by
$$\begin{eqnarray}\displaystyle {\mathcal{T}}(t):=\int _{0}^{1}(\unicode[STIX]{x1D719}(x,0)-\unicode[STIX]{x1D719}(x,t))\,\text{d}x={\displaystyle \frac{1}{k}}\int _{0}^{t}\int _{0}^{1}f(\unicode[STIX]{x1D719}(x,s))C(x,s)\,\text{d}x\,\text{d}s, & & \displaystyle\end{eqnarray}$$
where the second equality can be deduced from using the relationship
$\unicode[STIX]{x1D719}=1-2\unicode[STIX]{x03C0}R^{2}/\sqrt{3}$
to obtain
$\unicode[STIX]{x2202}\unicode[STIX]{x1D719}/\unicode[STIX]{x2202}t=-(|\unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{f}|/|\unicode[STIX]{x1D714}|)\unicode[STIX]{x2202}R/\unicode[STIX]{x2202}t$
, in conjunction with (3.18c
) and (4.4). A measure of the corresponding structural changes within the filter is captured by the quantity
$$\begin{eqnarray}\displaystyle P(x,t)=1-2R(x,t) & & \displaystyle\end{eqnarray}$$
which expresses the ratio of the distance between obstacles and the obstacle centres, or an effective pore size in the filter.
Figure 3. Filtration through an initially uniform porosity filter,
$\unicode[STIX]{x1D719}(x,0)=0.8$
, with
$Pe=3$
,
$k=1$
. (a) The pore size,
$P$
, defined by (4.8), throughout the filter at four snapshots in time (shown by asterisks in b–d). We start with
$t=0$
and end with the first time at which a pore size in the filter reaches
$0.01$
, defined as
$t_{0.01}$
. The two intermediate snapshots are at time
$t_{0.01}/3$
and
$2t_{0.01}/3$
, where we use
$t_{0.01}$
calculated from the solution to the full homogenized equations for all cases. (b) The cumulative contaminant adsorption
${\mathcal{T}}$
, defined by (4.7), as a function of time. (c) The time taken until the pore size reaches a minimum value. (d) The cumulative contaminant adsorption
${\mathcal{T}}$
, defined by (4.7), as a function of minimum pore size. The solid black curves denote the solutions to the full homogenized equations (4.1)–(4.6), the dashed red curves denote the asymptotic approximation when obstacle growth is small (4.12) and the dotted yellow curves denote the asymptotic solution when we further assume a weak porosity gradient (4.13). We see excellent agreement throughout, with only small deviations that occur when the pore size becomes small.
We begin by considering the case of a filter whose initial porosity is uniform. Solving (4.6) we find that as filtration proceeds the pore size
$P(x,t)$
decreases in time throughout the filter (figure 3
a) and the rate of contaminant removal
${\mathcal{T}}^{\prime }(t)$
reduces (figure 3
b). The reduction in pore size is more rapid close to the filter entrance than towards the filter exit (figure 3
a), which causes the filter to block at the entrance while the rest of the filter is still fit for purpose. This has significant practical disadvantages, and so we now use our model to explore strategies that minimize this wastefulness by varying the initial filter porosity. To do so, we first note that the homogenized system (4.6) can be further simplified by exploiting two relevant asymptotic regimes. We discuss these in the following section.
4.2 Asymptotic reductions of the homogenized equations
In the majority of filtration applications the growth of the obstacles, and thus the blocking of the filter, takes place over a timescale that is much longer than the time taken for the fluid to flow through the filter. Mathematically, this corresponds to the quasi-steady limit in which
$\unicode[STIX]{x1D6FC}\rightarrow 0$
. At leading order in this limit, returning to the three-dimensional equations for generality, the flow equations (3.8) and (3.10) become
$$\begin{eqnarray}\displaystyle \boldsymbol{U}(\boldsymbol{x},t)=-{\mathcal{K}}(\unicode[STIX]{x1D719})\unicode[STIX]{x1D735}_{\boldsymbol{x}}p,\quad \unicode[STIX]{x1D735}_{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{U}=0, & & \displaystyle\end{eqnarray}$$
where
${\mathcal{K}}(\unicode[STIX]{x1D719})$
is defined in (3.8b
); and the transport equations (3.18) become
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}_{\boldsymbol{x}}\boldsymbol{\cdot }\left({\mathcal{D}}(\unicode[STIX]{x1D719})\unicode[STIX]{x1D735}_{\boldsymbol{x}}C-{\displaystyle \frac{C}{\unicode[STIX]{x1D719}}}(\boldsymbol{U}(\unicode[STIX]{x1D719})+{\mathcal{D}}(\unicode[STIX]{x1D719})\unicode[STIX]{x1D735}_{\boldsymbol{x}}\unicode[STIX]{x1D719})\right)=f(\unicode[STIX]{x1D719})C,\quad \unicode[STIX]{x1D719}\frac{\unicode[STIX]{x2202}R}{\unicode[STIX]{x2202}t}=C, & & \displaystyle\end{eqnarray}$$
where
${\mathcal{D}}(\unicode[STIX]{x1D719})$
and
$f(\unicode[STIX]{x1D719})$
are defined in (3.18b,c
), respectively. In this regime, the growth in
$R$
decouples from the flow and transport problems, which are now quasi-steady.
For a unidirectional porosity variation, the flow velocity becomes independent of
$x$
and so (4.9a
) simplifies to
$$\begin{eqnarray}\displaystyle U(t)=\overline{K}(t)\unicode[STIX]{x0394}P, & & \displaystyle\end{eqnarray}$$
where
$\overline{K}$
is given in (4.2). Moreover, the upstream and downstream concentrations can be solved to obtain
$$\begin{eqnarray}\displaystyle C^{-}=1-A(t)\exp (PeU(t)x),\quad C^{+}=B(t), & & \displaystyle\end{eqnarray}$$
where
$A(t)$
and
$B(t)$
are arbitrary functions of time. Substituting (4.11) into (4.6), the system for
$C$
becomes
$$\begin{eqnarray}\displaystyle & \displaystyle f(\unicode[STIX]{x1D719})C=\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left[{\mathcal{D}}(\unicode[STIX]{x1D719})\frac{\unicode[STIX]{x2202}C}{\unicode[STIX]{x2202}x}-{\displaystyle \frac{C}{\unicode[STIX]{x1D719}}}\left(U(t)+{\mathcal{D}}(\unicode[STIX]{x1D719})\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}x}\right)\right],\quad x\in (0,1), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle -U(t)={\mathcal{D}}(\unicode[STIX]{x1D719})\frac{\unicode[STIX]{x2202}C}{\unicode[STIX]{x2202}x}-{\displaystyle \frac{C}{\unicode[STIX]{x1D719}}}\left(U(t)+{\mathcal{D}}(\unicode[STIX]{x1D719})\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}x}\right),\quad x=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle 0=\frac{\unicode[STIX]{x2202}C}{\unicode[STIX]{x2202}x}-{\displaystyle \frac{C}{\unicode[STIX]{x1D719}}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}x},\quad x=1. & \displaystyle\end{eqnarray}$$
The governing equation for the porosity (4.4) is unchanged:
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}\frac{\unicode[STIX]{x2202}R}{\unicode[STIX]{x2202}t}=C,\quad x\in (0,1). & & \displaystyle\end{eqnarray}$$
This problem is similar to that considered in Dalwadi et al. (Reference Dalwadi, Griffiths and Bruna2015), where the absence of microstructural changes due to filter clogging meant that
$U(t)\equiv 1$
in (4.12), and
$\unicode[STIX]{x2202}R/\unicode[STIX]{x2202}t\equiv 0$
instead of (4.12d
). We can solve the system (4.12) numerically using the method of lines, discretizing in space with a second-order finite difference scheme and using the MATLAB program ode15s to solve in time.
We can make further analytic progress by considering the case where, in addition to the slow obstacle growth, the deviation of the porosity
$\unicode[STIX]{x1D719}(x,t)$
from its average value in space
$\bar{\unicode[STIX]{x1D719}}(t)=\int _{0}^{1}\unicode[STIX]{x1D719}(x,t)\,\text{d}x$
is small. In this case, we express the porosity as its average in space plus the deviation, as follows
$\unicode[STIX]{x1D719}(x,t)=\bar{\unicode[STIX]{x1D719}}(t)+\unicode[STIX]{x1D6F7}(x,t)$
, and assume that
$\max _{x\in (0,1)}|\unicode[STIX]{x1D6F7}|\ll \bar{\unicode[STIX]{x1D719}}$
for any given time. Thus, the leading-order coefficients in the linear governing equation (4.12a
) are constant in
$x$
, and the solution can be written in terms of hyperbolic functions. The
$O(\unicode[STIX]{x1D6F7})$
correction terms in (4.12a
) can then be solved using the method of variation of parameters (a similar calculation is carried out in Dalwadi et al.
Reference Dalwadi, Griffiths and Bruna2015). From this, we may express the concentration
$C$
as a sum of the concentration distribution in a filter with constant porosity,
$C_{0}(x,t)$
, and the adjustment due to the porosity gradient,
$C_{1}(x,t)$
, i.e.
$$\begin{eqnarray}\displaystyle C(x,t)\sim C_{0}(x,t)+C_{1}(x,t), & & \displaystyle\end{eqnarray}$$
where
$|C_{1}|\ll C_{0}$
;
$C_{0}$
and
$C_{1}$
are functions of
$\unicode[STIX]{x1D719}$
that can be solved for explicitly and are given in appendix A. Using this analytic solution, we are able to reduce the entire problem to solving (4.12d
). We solve this system numerically using the method of lines, discretizing in space with a second-order finite difference scheme and using the MATLAB program ode15s to solve in time.
In figure 3, we compare the results from these two asymptotic results, (4.12) and (4.13), with the full numerical solution, (4.1)–(4.6) (solid black lines), as dashed red and dotted yellow lines, respectively. Both asymptotic solutions offer very good agreement in the pore size
$P$
with the full numerical solution (figure 3
a), and in
${\mathcal{T}}$
until around
$t=0.45$
(figure 3
b).
The discrepancy between numeric and asymptotic solutions appears when considering the time taken to reach a minimum pore size rather than in the cumulative contaminant removal metric
${\mathcal{T}}$
. We are able to contextualize this metric by introducing relative timings in a filter lifespan. That is, we define the minimum pore size
$$\begin{eqnarray}\displaystyle P^{\ast }(t)=\min _{x\in (0,1),t\geqslant 0}P(x,t), & & \displaystyle\end{eqnarray}$$
and the time taken until the minimum pore size reaches a given value,
$$\begin{eqnarray}\displaystyle t_{P^{\ast }}=\min \{t\geqslant 0:P(x,t)\leqslant P^{\ast },\text{ for some }x\in (0,1)\}. & & \displaystyle\end{eqnarray}$$
Blocking occurs at
$t_{0}$
, but we must consider small finite values of
$P^{\ast }$
due to the singularity in (4.12) when
$\unicode[STIX]{x1D719}\rightarrow 0$
, and hence when
$P^{\ast }\rightarrow 0$
, at the point of blocking. For example, in figure 3 the simulations are stopped at
$t_{0.01}$
.
From figure 3(b), we see that although the maximum time to removal is different between the asymptotic and numerical solutions, the total cumulative removals are very close. We show the former in more detail in figure 3(c), where we see that the numerical and asymptotic solutions agree well in
$t_{P^{\ast }}$
until
$P^{\ast }$
falls below around
$P^{\ast }=0.05$
. We show the latter in more detail in figure 3(d), where we use
$P^{\ast }$
as the abscissa instead of time, against cumulative contaminant removal. We see that there is excellent agreement between numeric and asymptotic solutions throughout, and henceforth use the full asymptotic approximation (4.13), which offers a very accurate representation of the solution to the full homogenized system (4.1)–(4.6) while enabling further analytic study of the system. Moreover, we note that we halt simulations when
$P^{\ast }$
reaches
$0.01$
, and refer to this as ‘blocking’, and we see from figure 3(d) that this is unlikely to yield major differences in values of total contaminant removal.
4.3 Linearly graded filters
As shown in figure 3, in a depth filter with a uniform porosity the early portion of the media will block before the later portion has been used to its full effect. This can be counteracted by using a filter whose porosity decreases with depth.
In this section we consider linearly graded filters at
$t=0$
, characterized by
$\unicode[STIX]{x1D719}(0,x)=\unicode[STIX]{x1D719}_{0}+m(x-0.5)$
, where we vary the initial average porosity,
$\unicode[STIX]{x1D719}_{0}$
, and the initial porosity gradient,
$m$
. In a similar manner to Dalwadi et al. (Reference Dalwadi, Griffiths and Bruna2015), we assume that the pressure drop is kept constant across each filter considered and, noting that the system is invariant to the scaling
$$\begin{eqnarray}\displaystyle (U(0),Pe,k)\mapsto (\unicode[STIX]{x1D714}U(0),Pe/\unicode[STIX]{x1D714},\unicode[STIX]{x1D714}k), & & \displaystyle\end{eqnarray}$$
for any constant
$\unicode[STIX]{x1D714}$
, we impose
$U(0)=1$
and vary
$Pe$
and
$k$
accordingly. For short times, a filter with spatially uniform porosity (i.e.
$m=0$
) provides the largest cumulative removal for a given average porosity, and there is an optimum value for this average porosity (figure 4
a). This optimum exists because a highly porous filter has less surface area with which to remove contaminant, but a less porous filter will admit a lower flow rate through the filter, relying more on diffusion than advection for contaminant transport and thus lowering the contaminant removal. The result that there is no apparent difference in the total contaminant removal of a filter with a positive or negative porosity gradient for small time replicates the results for the initial instantaneous adsorption obtained in Dalwadi et al. (Reference Dalwadi, Griffiths and Bruna2015). However, in contrast to Dalwadi et al. (Reference Dalwadi, Griffiths and Bruna2015) we are now able to capture the asymmetry in removal that arises over longer times as a result of blockages within the filter. In particular, negative initial porosity gradients do provide an inherently larger cumulative removal at a given time than positive initial porosity gradients (figure 4
b,c), and their final cumulative removal is larger (figure 4
d). We note that whilst negative initial porosity gradients do last for a longer time before blocking, the filters that take the longest time to block do not automatically remove the most contaminant (figure 5).
In Dalwadi et al. (Reference Dalwadi, Griffiths and Bruna2015), a further metric was introduced to measure the uniformity of uptake, thus allowing us to account for the superior performance of filters with a negative porosity gradient without explicitly accounting for filter evolution. As we are able to account for filter evolution in this work,
${\mathcal{T}}$
is now the better measure of filtration. We show in appendix B that this alternative metric can only provide qualitative insights into filtration, and should not be used when the filter structure evolves in time.
Figure 4. Cumulative contaminant removal,
${\mathcal{T}}$
, defined by (4.7), for varying
$\unicode[STIX]{x1D719}_{0}$
and
$m$
in initial porosity distributions of the form
$\unicode[STIX]{x1D719}(x)=\unicode[STIX]{x1D719}_{0}+m(x-0.5)$
, given by the solution to (4.1)–(4.6): (a)
$t=10^{-4}$
; (b)
$t=0.1$
; (c)
$t=0.7$
; (d)
$t=2$
. Each curve represents a different value of
$\unicode[STIX]{x1D719}_{0}$
, which increments in steps of
$0.1$
from
$\unicode[STIX]{x1D719}_{0}=0.4$
to
$\unicode[STIX]{x1D719}_{0}=0.8$
and the arrow in (d) denotes curves of increasing
$\unicode[STIX]{x1D719}_{0}$
. We vary
$m$
for a given
$\unicode[STIX]{x1D719}_{0}$
such that
$\unicode[STIX]{x1D719}\in [0.2,0.95]$
. Therefore, the available range of
$m$
varies with
$\unicode[STIX]{x1D719}_{0}$
. We use the reference values
$\unicode[STIX]{x1D719}=0.6$
,
$Pe=5$
,
$k=0.1$
from which to modify appropriate parameters. The dashed curves denote that the minimum pore size has reached
$0.01$
and the filtration has stopped.
Figure 5. (a) The time until the minimum pore size reaches 0.01 for varying
$\unicode[STIX]{x1D719}_{0}$
and
$m$
in initial porosity distributions of the form
$\unicode[STIX]{x1D719}(x)=\unicode[STIX]{x1D719}_{0}+m(x-0.5)$
, determined by solving (4.1)–(4.6). The dashed black lines denote the times at which snapshots are taken in figures 4 and 7. (b) The total cumulative adsorption over the lifespan of a filter versus the time taken to block. Each line represents a different initial average porosity, and the unlabelled arrows denote the direction of increasing
$m$
for each line. In both figures, each line represents a different value of
$\unicode[STIX]{x1D719}_{0}$
, which increments in steps of
$0.1$
from
$\unicode[STIX]{x1D719}_{0}=0.4$
to
$\unicode[STIX]{x1D719}_{0}=0.8$
and the labelled arrow denotes lines of increasing
$\unicode[STIX]{x1D719}_{0}$
. We vary
$m$
for a given
$\unicode[STIX]{x1D719}_{0}$
such that
$\unicode[STIX]{x1D719}\in [0.2,0.95]$
. Therefore, the available range of
$m$
varies with
$\unicode[STIX]{x1D719}_{0}$
. We use the reference values
$\unicode[STIX]{x1D719}=0.6$
,
$Pe=5$
,
$k=0.1$
from which to modify appropriate parameters.
For a given initial average porosity, the largest cumulative removal does not necessarily correspond to a filter with the most negative initial porosity gradient (figure 4
d). If the porosity gradient is too negative, the removal is skewed towards the filter exit, and thus blocking occurs at the filter exit too quickly. This effect causes the sharp drop-off in filter performance seen in figure 4(c,d) for large negative values of
$m$
. Thus, an optimum initial porosity gradient exists for a given initial average porosity. Whilst this ‘optimal’ initial porosity gradient maximizes contaminant removal over the space of initially linear graded filters with a given initial average porosity, it will not necessarily maximize contaminant removal over the space of all filters with a given initial average porosity. We explore this idea in more detail in the next section.
4.4 Filters that block everywhere at once
Using the framework we have set up in this paper, we can address the question of which initial porosity distributions maximize the cumulative removal for given operating conditions. We are able to bypass much of the difficulty involved in the full optimal control problem by noting that, for a given initial average filter porosity, the cumulative adsorption (4.7) is maximized when a filter blocks all at once, rather than in just one place. Most filters will only block in one place, and this is the case for all of the filters represented in figure 4, where the filters block at either the filter entrance or the filter exit.
By starting with a filter that is blocked everywhere and running our simulations backwards in time, we can obtain initial porosity distributions whose lifespans end with blocking occurring everywhere in the filter at once. To avoid the computational issues arising when blocking occurs, we approximate the concept of a filter that blocks everywhere at once by a filter that reaches a pore size of
$0.01$
everywhere at once. An issue with this procedure is that the full homogenized system for the concentration (4.1)–(4.6) is parabolic in time, and thus running the simulations backwards in time results in an ill-posed problem. However, the evolution equation for the filter porosity (4.12d
) is well-posed if run backwards in time, and only relies on knowing the concentration for a given porosity. Thus, if we use the asymptotic solutions for the concentration derived in § 4.2, we do not have to run the concentration evolution equation backwards in time and we are able to side-step the issue with the ill-posedness of this problem.
A related issue, due to solving the system of equations backwards in time, is that only the final flow velocity can be imposed, and the initial flow velocity is now an output to the calculation. Thus, a general final flow velocity
$U(t_{0.01})$
will not necessarily result in
$U(0)=1$
, the condition we prescribe in the forward problem. Moreover, we vary
$Pe$
and
$k$
according to the initial porosity in the forward problem, using the invariant scaling (4.16), and thus it is not immediately clear what values of
$Pe$
and
$k$
we should impose for the inverse problem. We overcome these problems by exploiting the fact that we impose a constant pressure drop across the filter, and thus
$\overline{K}(0)=U(0)\overline{K}(t_{0.01})/U(t_{0.01})=\overline{K}(t_{0.01})/U(t_{0.01})$
, where the first equality holds from (4.10), and the second equality holds by using
$U(0)=1$
(the filter for which we are aiming). This allows us to move between the initial and final average permeabilities as a function of the final flow velocity, as well as imposing values of
$Pe$
and
$k$
using the final porosity. Hence, we are able to use a shooting method for just one parameter, the final flow velocity
$U(t_{0.01})$
, shooting to achieve
$U(0)=1$
.
To summarize, we impose a final constant pore size with
$P(x,t_{0.01})\equiv 0.01$
and a final flow velocity and run our simulation in reverse, halting the process when the porosity at any point reaches a set maximum. We iterate this process, varying the final flow velocity, until the initial flow velocity
$U(0)=1$
to within an accuracy of
$2\times 10^{-3}$
. We repeat this process for different maximum porosity constraints to obtain a family of initial porosity distributions that block everywhere at once when contaminant flows through, each with a different maximum (and average) porosity. The maximum porosities we choose are given by
$1-\unicode[STIX]{x03C0}R^{2}$
, where
$R$
takes the values
$0.1174$
,
$0.15$
,
$0.2$
,
$0.25$
and
$0.3$
. We produce this range of initial porosity distributions as there may be physical constraints on how large/small the porosity/fibres can be.
Each resulting filter has a decreasing porosity with depth and, compared with the region near the filter entrance, the porosity gradient in each filter is more negative towards the centre of the filter and smaller towards the end of the filter (figure 6
a). This effect is more pronounced for larger maximum porosities. We compare how these filters remove contaminant against the linear filters we considered in § 4.3 and see that the ideal filters do remove around
$20\,\%$
more contaminant for a given initial average porosity (figure 6
b). However, these filters cannot have arbitrarily large initial average porosity, as a filter whose initial porosity is everywhere close to
$1$
will not block everywhere at once, and thus a linearly graded filter with a large enough initial average porosity can remove more contaminant than a filter with a lower initial average porosity that blocks everywhere at once. As a result, we may term these pseudo-optimal filters. As a filter may require a maximum porosity to maintain structural integrity, the choice of initial porosity distribution will depend on the operating conditions.
Figure 6. (a) The initial porosity distributions for filters that block everywhere at once with different initial average porosities,
$\unicode[STIX]{x1D719}_{0}$
. The initial average porosities are: 0.320, 0.390, 0.465, 0.548 and 0.607 (all to three significant figures). (b) The cumulative contaminant removal
${\mathcal{T}}$
for the pseudo-optimal filters (red asterisks) compared with the linear filters considered in § 4.3 (black crosses). The data for the black crosses is from figure 4. We use dimensionless parameters corresponding to the reference values
$\unicode[STIX]{x1D719}(x,0)\equiv 0.6$
,
$Pe=5$
,
$k=0.1$
in the manner described in the text. In (b) we see that the pseudo-optimal filters do provide increased contaminant removal for a given initial average porosity, but the average porosity of the pseudo-optimal filters cannot be increased to an arbitrarily large value.
5 Discussion
We have systematically derived a macroscopic model for the dynamic blocking of a porosity-graded filter from microscale information by a generalization of standard homogenization theory for near-periodic systems. The result is an advection–diffusion–reaction equation for the solute concentration within the filter, a modified Darcy law for the fluid transport and an evolution equation for the underlying porosity of the filter. This system depends on the initial porosity distribution of the filter and the operating conditions. In particular, our theory has allowed us to determine solute trapping across the lifetime of a filter whose pores are constricting due to contaminant removal. We have also solved the inverse problem of calculating the initial porosity distributions that lead to a filter that blocks uniformly for given operating conditions.
To track the evolution of filter porosity, we have accounted for a macroscale variation in our filter, using the generalized homogenization technique of Bruna & Chapman (Reference Bruna and Chapman2015). The resulting macroscale equations allow us to efficiently quantify filter performance. By performing a parameter sweep over the functional space of initially linear porosity functions, we have quantified the experimental observation that filters with an initially negative porosity gradient are more effective at removing contaminant over time than filters with an initially zero (or, worse, positive) porosity gradient.
In addition, our asymptotic reductions in the limit of slow filter blocking and weak spatial porosity variation have allowed us to solve an inverse problem to calculate the initial filter porosity distributions that lead to a filter that blocks uniformly for given operating conditions. We show that these filters provide a greater contaminant removal than a linearly graded filter with the same initial average porosity, but there is a maximum initial average porosity that these filters can reach. As a result, we term these pseudo-optimal filters.
In Dalwadi et al. (Reference Dalwadi, Griffiths and Bruna2015) we showed that, for contaminants that did not alter the geometry of the filter as they adhered, a negative porosity gradient could result in a near-uniform contaminant removal in space, compared with a significantly asymmetric removal for a positive porosity gradient. We conjectured that, were this contaminant removal to result in pore constriction over time, an initially negative porosity gradient would outlast an initially positive porosity gradient. In this paper, we have confirmed and quantified this conjecture. Furthermore, we were able to deduce that a filter with an initially negative porosity gradient yields a larger cumulative contaminant adsorption than an initially positive porosity gradient at a given time, even before blocking occurs.
In this paper we considered a microscale structure consisting of cylindrical fibres, resulting in a method of varying the local porosity through a single parameter (the fibre radius), and an explicit macroscale equation. However, the homogenization procedure we use can be readily extended to a general microstructure, and one could combine this work with that of Richardson & Chapman (Reference Richardson and Chapman2011), who considered a general curvilinear coordinate transform to map a near-periodic microscale to a periodic domain, thus allowing a homogenization method to be applied. In general, the coefficients within the resulting governing equations would have to be determined at each point in space, and thus the resulting macroscale equations are more complicated to solve than in the case presented here. However, there is a computationally efficient middle ground between cylinders and a general microstructure. For instance, one could choose a one-parameter family for the potential microstructure, which allows the macroscale coefficients to be written as functions of the porosity. However, we note that obstacle shape is largely unimportant for high porosities (Bruna & Chapman Reference Bruna and Chapman2015) and is unlikely to be preserved as particles are adsorbed. Indeed, one may intuitively expect initially non-circular obstacles to tend to a more circular shape as particles are adsorbed, and thus the system would tend towards the one considered in this paper.
One simplifying assumption that we have made in this work is that the contaminant trapping rate is linearly dependent on contaminant concentration, and that the trapping will continue until the microscale cell is fully blocked. In many adsorption models, different trapping conditions, for example, the Langmuir adsorption model, are used (see, for example, Morel & Hering Reference Morel and Hering1993). This model assumes that there are a finite number of adsorption sites on the obstacle surface (thus, the number of sites is proportional to the obstacle surface area), and that the adsorption layer is only one contaminant molecule thick. Thus, while the adsorption rate will be approximately linearly dependent on the contaminant concentration when few adsorption sites are in use, the adsorption rate will decrease to zero as more adsorption sites are used. This adsorption model can easily be included in this work, and the restriction on the number of adsorption sites means that blocking would not occur in most cases. In this case, the filter that maximizes adsorption would be the one with the largest number of adsorption sites, that is, the filter with the maximum surface area.
While we have shown how to maximize contaminant removal across the functional space of filters with the same initial average porosity, our results only apply up to some maximum initial average porosity. This is because filters that block everywhere at once will have an initially decreasing porosity, and the porosity can never exceed one. Thus, the full control problem for any given initial average porosity is still an open problem. This is an important question that would be useful to tackle in future work.
Finally, we note that this work has the potential not only to guide filter manufacture and operating conditions, but also to provide assistance to many other industries. We have introduced a general framework in this paper, which applies to any problem where the underlying transport satisfies an advection–diffusion equation with a general adsorption condition on the microstructure surface. For example, our model can predict drug transport and delivery to tumours, and a simple change of sign to the evolution equation for the porosity will allow prediction of the resultant tumour shrinkage. Using appropriate parameter values, one can predict the effect of various drugs, significantly aiding the task of testing new drug therapies.
Acknowledgements
This work was carried out whilst M.P.D. was funded by an EPSRC Impact Acceleration Account Award (grant no. EP/K503769/1) at the University of Oxford. I.M.G. is supported by a Royal Society Fellowship. M.B. is funded by a Junior Research Fellowship of St John’s College, Oxford. This work was also supported by the 2020 Science programme, which was funded through the EPSRC Cross-Disciplinary Interface Programme (grant no. EP/I017909/1).
Appendix A. Asymptotic results for the concentration field
The functions
$C_{0}$
and
$C_{1}$
are defined as
$$\begin{eqnarray}\displaystyle C_{0}(x,t) & = & \displaystyle 2\unicode[STIX]{x1D6FD}(t)\bar{\unicode[STIX]{x1D719}}(t)\exp (a(t)x)[b(t)\cosh a(t)b(t)(x-1)-\sinh a(t)b(t)(x-1)],\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle C_{1}(x,t) & = & \displaystyle {\displaystyle \frac{\exp (a(t)x)}{\sinh a(t)b(t)}}\left[(A_{1}(x,t)+\unicode[STIX]{x1D6FD}(t)B_{1}(t))\sinh (a(t)b(t)x)\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\,(A_{2}(x,t)+\unicode[STIX]{x1D6FD}(t)B_{2}(t))\sinh (a(t)b(t)(x-1))\right],\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle a(t)={\displaystyle \frac{U(t)}{2\bar{\unicode[STIX]{x1D719}}(t){\mathcal{D}}(\bar{\unicode[STIX]{x1D719}})}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle b(t)=\sqrt{1+f(\bar{\unicode[STIX]{x1D719}})/\left(a^{2}(t){\mathcal{D}}(\bar{\unicode[STIX]{x1D719}})\right)}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FD}(t)={\displaystyle \frac{1}{(1+b^{2}(t))\sinh a(t)b(t)+2b(t)\cosh a(t)b(t)}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle A_{1}(x,t)={\displaystyle \frac{1}{a(t)b(t){\mathcal{D}}(\bar{\unicode[STIX]{x1D719}})}}\int _{x}^{1}g(\unicode[STIX]{x1D709},t,C_{0}(\unicode[STIX]{x1D709},t))\exp (-a\unicode[STIX]{x1D709})\sinh ab(\unicode[STIX]{x1D709}-1)\,\text{d}\unicode[STIX]{x1D709}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle A_{2}(x,t)={\displaystyle \frac{1}{a(t)b(t){\mathcal{D}}(\bar{\unicode[STIX]{x1D719}})}}\int _{0}^{x}g(\unicode[STIX]{x1D709},t,C_{0}(\unicode[STIX]{x1D709},t))\exp (-a\unicode[STIX]{x1D709})\sinh ab\unicode[STIX]{x1D709}\,\text{d}\unicode[STIX]{x1D709}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \left(\begin{array}{@{}c@{}}\displaystyle B_{1}(t)\\ \displaystyle B_{2}(t)\end{array}\right)=\left(\begin{array}{@{}cc@{}}\displaystyle -q(t) & \displaystyle b(t)\\ b(t) & \displaystyle -q(t)\end{array}\right)\left(\begin{array}{@{}c@{}}\displaystyle b(t)\left(A_{2}(1,t)-2\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F7}(1,t)}{\unicode[STIX]{x2202}x}\right)\\ \displaystyle b(t)A_{1}(0,t)+N(0,C_{0}(0,t))/(a(t){\mathcal{D}}(\bar{\unicode[STIX]{x1D719}}))\end{array}\right),\quad & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle g(x,t,C)=-\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}N(x,t,C)+\unicode[STIX]{x1D6F7}f(\bar{\unicode[STIX]{x1D719}})C, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle N(x,t,C)=\unicode[STIX]{x1D6F7}{\mathcal{D}}(\bar{\unicode[STIX]{x1D719}})\frac{\unicode[STIX]{x2202}C}{\unicode[STIX]{x2202}x}-{\displaystyle \frac{C}{\bar{\unicode[STIX]{x1D719}}}}\left({\mathcal{D}}(\bar{\unicode[STIX]{x1D719}})\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F7}}{\unicode[STIX]{x2202}x}-{\displaystyle \frac{\unicode[STIX]{x1D6F7}}{\bar{\unicode[STIX]{x1D719}}}}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle q(t)=\sinh (a(t)b(t))+b(t)\cosh (a(t)b(t)). & \displaystyle\end{eqnarray}$$
Appendix B. Uniformity of removal
In Dalwadi et al. (Reference Dalwadi, Griffiths and Bruna2015), blocking was not explicitly accounted for, and thus the superior performance of filters with a negative porosity gradient was accounted for by introducing a metric that measured the uniformity of uptake. The equivalent metric in our time-dependent case is
$$\begin{eqnarray}\displaystyle {\mathcal{S}}(t)={\displaystyle \frac{1}{k}}\int _{0}^{1}\left|f(\unicode[STIX]{x1D719}(x,t))C(x,t)-\int _{0}^{1}f(\unicode[STIX]{x1D719}(s,t))C(s,t)\,\text{d}s\right|\,\text{d}x\geqslant 0, & & \displaystyle\end{eqnarray}$$
where a lower value of
${\mathcal{S}}$
represents a more uniform removal and
${\mathcal{S}}=0$
represents uniform uptake across the filter. As with the cumulative removal, for short times we observe similar results for
${\mathcal{S}}$
as obtained in Dalwadi et al. (Reference Dalwadi, Griffiths and Bruna2015) (figure 7
a), with initial negative porosity gradients generally yielding more uniform removal than positive porosity gradients. However, as time increases, this measure becomes much less useful (figure 7
b–d). Whilst lower values of
${\mathcal{S}}(t)$
are useful indicators of larger removal, they do not necessarily correspond to the largest values of removal (comparing figure 7
c with figure 4
c). The instantaneous nature of
${\mathcal{S}}$
does not convey enough global information about the cumulative removal, and factors such as a lower average porosity are also important.
Figure 7. The instantaneous uniformity of adsorption,
${\mathcal{S}}$
, defined by (B 1), for varying
$\unicode[STIX]{x1D719}_{0}$
and
$m$
in initial porosity distributions of the form
$\unicode[STIX]{x1D719}(x)=\unicode[STIX]{x1D719}_{0}+m(x-0.5)$
, determined by solving (4.1)–(4.6): (a)
$t=10^{-4}$
; (b)
$t=0.1$
; (c)
$t=0.7$
; (d)
$t=2$
. Each curve represents a different value of
$\unicode[STIX]{x1D719}_{0}$
, which increments in steps of
$0.1$
from
$\unicode[STIX]{x1D719}_{0}=0.4$
to
$\unicode[STIX]{x1D719}_{0}=0.8$
and the arrow in (a) denotes curves of increasing
$\unicode[STIX]{x1D719}_{0}$
. We vary
$m$
for a given
$\unicode[STIX]{x1D719}_{0}$
such that
$\unicode[STIX]{x1D719}\in [0.2,0.95]$
. Therefore, the available range of
$m$
varies with
$\unicode[STIX]{x1D719}_{0}$
. We use the reference values
$\unicode[STIX]{x1D719}=0.6$
,
$Pe=5$
,
$k=0.1$
from which to modify appropriate parameters. The dashed curves denote that the minimum pore size has reached
$0.01$
and the filtration has stopped.
