## Abstract

Filters whose porosity decreases with depth are often more efficient at removing solute from a fluid than filters with a uniform porosity. We investigate this phenomenon via an extension of homogenization theory that accounts for a macroscale variation in microstructure. In the first stage of the paper, we homogenize the problems of flow through a filter with a near-periodic microstructure and of solute transport owing to advection, diffusion and filter adsorption. In the second stage, we use the computationally efficient homogenized equations to investigate and quantify why porosity gradients can improve filter efficiency. We find that a porosity gradient has a much larger effect on the uniformity of adsorption than it does on the total adsorption. This allows us to understand how a decreasing porosity can lead to a greater filter efficiency, by lowering the risk of localized blocking while maintaining the rate of total contaminant removal.

## 1. Introduction

Membrane separation is a vast industry with a wide range of applications, including water treatment [1,2], biopharmaceuticals [3,4] and food processing [5,6]. For example, filters are crucial to remove waste and excess water from the blood in kidney dialysis [7], and yeast and bacteria in beer production [8]. Despite the diverse industrial applications, an overarching goal of membrane design is to maximize the product yield. In applications where the solvent is the desired product, such as in water treatment, this is accomplished by maximizing both particle removal from the fluid suspension and filter lifespan. Mathematical modelling can offer key insight into the filtration process and operating conditions, and thus provide a cost-effective way to optimize filter design.

A common form of membrane separation is *depth filtration*, in which small particles (contaminants) are trapped within, and not just on the surface of, the porous filter material. Such filters often capture the majority of the contaminants in the initial portion of the filter while leaving the latter portions relatively unused, leading to premature clogging and reduced filtration efficiency [9]. In such cases, filters whose porosity decreases with depth, or *porosity-graded filters*, which are an example of a functionally graded material, can improve filtration efficiency, and so are often used experimentally [10–14]. Their increased efficiency can be qualitatively attributed to a decrease in porosity compensating for a reduction of contaminant concentration with depth, which occurs owing to prior filtering. However, the mechanism behind this observation is not fully understood. Experiments are costly, and, moreover, it is very difficult to observe the particle trapping within a filter during the filtration process directly. Instead, deductions are made only after dissecting the porous medium once filtration has ceased. For these reasons, optimizing a porosity distribution via a systematic series of experiments is impractical.

Mathematical and computational methods to model the filtration process are very useful for investigating this issue at a fraction of the full experimental cost. A summary table of computational fluid dynamics (CFD) models for filtration is given in [15]. In the same paper, the authors use scanning electron microscopes to obtain a full description of a given membrane microstructure, and implement a full CFD model of particles moving within the membrane. While a full CFD simulation gives us excellent insight into how an individual particle is trapped, there are large computational costs associated with keeping track of all the particles within a complicated pore structure. Additionally, as in the example given above, each membrane must be scanned to accurately represent its pore structure. Although filtration occurs on the scale of the particle or pore size, it is generally the overall macroscale behaviour, such as the total mass of particles removed, which is the main concern in filter design.

An alternative approach to CFD for complex heterogeneous materials, such as filters, is to consider upscaling methods. These reduce the complexity of an equation by averaging any microscale variation while retaining the important macroscale variation. Mathematical homogenization via the method of multiple scales is an asymptotic technique often used for this upscaling. Traditionally, for this technique to be appropriate to use, there must be a periodic microstructure and the ratio between the length of the periodic cell and the length of the important macroscale variation must be small. In fluid mechanics, where the underlying flow equations can be difficult to solve, this procedure is used to determine the macroscale behaviour of fluids in complicated geometries. For example, in saturated single-phase Stokes flow past a periodic array of inert obstacles, this homogenization procedure leads to Darcy's law [16], and, in non-saturated media, homogenization techniques can be used to derive Richards' equation [17].

While the majority of homogenization procedures are carried out in domains whose microstructure is fully periodic, and thus would only be applicable to filters with a constant macroscale porosity, extensions to non-periodic domains have been explored. In [18], the steady problem of nutrient uptake past randomly placed point sinks is considered in one spatial dimension. As the governing equations can be solved if the locations of the sinks are known, significant analytical progress is made into investigating the macroscale effect of different random distributions. A comprehensive extension to near-periodic domains is developed in [19], where the authors use a general curvilinear coordinate transform to homogenize the electric potential within a beating heart, mapping the near-periodic microscale to a periodic domain. Similar approaches are used to develop homogenized equations in [20], where the authors consider saturated flow in a poroelastic solid with surface growth of the solid phase, which is allowed to have a macroscopic variation, and in [21], where the authors consider blood and drug exchange across capillary walls in malignant tumours.

Although these extensions to non-periodic domains are certainly significant advancements to the field, their generality means that, usually, the cell problem must be solved at every point in the macroscale, reducing the computational efficiency. This issue is bypassed in [22], where the multiple scales homogenization method is extended to quasi-periodic structures where the microstructure is allowed to vary slowly. In particular, the authors consider the problem of diffusion in a porous medium whose solid phase is modelled as an array of inert spherical inclusions (obstacles). The radii of these obstacles are allowed to vary spatially over the macroscale length. By imposing a specific one-parameter form on the obstacles, an explicit formulation of the macroscale equation is obtained in terms of the cell-averaged porosity, which varies in the macroscale. A notable aspect of this analysis is that a porosity variation induces a macroscale particle advection in the direction of decreasing porosity.

In this paper, we use the homogenization method developed in [19,22] to investigate porosity-graded filters. We model the filter material as a collection of spherical obstructions, whose radii vary spatially over a long lengthscale, past which a fluid with suspended contaminant particles flows owing to an applied transmembrane pressure. The particles are transported via advection and diffusion, and can be trapped on the filter microstructure. In general, the accumulation of contaminant particles via trapping modifies the filter geometry, and this process eventually leads to pore blockage within the filter. Because this effect occurs within and not just on the inlet surface of the filter, the blockage is difficult to remove. Dealing with a blockage causes a significant reduction in efficiency and an increase in maintenance cost, as the filtration procedure must be halted to perform a back-flow procedure to dislodge the blockage or, more drastically, a filter replacement.

The underlying trapping or adsorption mechanism of a particle to the filter structure will, in principle, depend on the solute solubility and the filter material adsorption. Many existing models for trapping in filters use constitutive macroscale laws to determine how fluid volume throughput varies in time, which are fitted with experimental data [23]. In recent work by Griffiths *et al*. [24], the authors perform simulations of a microscale model for the blocking of a network of pores formulated from a more fundamental mechanistic level, and their results are consistent with the set of constitutive laws. However, the model in [24] does not account for the flow profile and is computationally expensive as it keeps track of each pore. Incorporating pore blockage and the subsequent filter failure is mathematically challenging to include in homogenized models, because pore growth leads to difficult free-boundary problems.

The main objective of this work is to understand and quantify how porosity-graded filters improve filter efficiency. To this end, starting from the flow and particle transport problems in the complicated domain described by a given microstructure, we employ the homogenization method developed in [19,22] to systematically determine an effective macroscale equation. Our model takes into account the effects of the microstructure on the fluid flow, as well as the transport by diffusion and advection of the contaminant particles. The trapping of particles by the filter microstructure is also included as an adsorption process. As in [22], we consider a solid structure composed of spherical inclusions whose radii vary slowly in the macroscale. However, now the inclusions are partially adsorbing sinks for the particle transport and inert obstacles to the flow. In order to focus on the effect of variations in porosity within a filter, in this work, we assume that the contaminant particles are small and in dilute suspension within the fluid, and thus their trapping has a negligible effect on the obstacle size. This means that our model will not be able to explicitly capture pore blocking. As a result of this simplification, the flow problem decouples from the particle transport problem (but not vice versa), and we do not have a free-boundary problem. As we consider a dilute suspension, we follow the lead of Polyakov [25] and impose that particle adsorption is linearly dependent on the number of particles.

This paper is divided into two stages. The first stage consists of the presentation of the full problem in §2, and the homogenization of the flow and particle transport problems to systematically obtain effective equations on the macroscale in §3. This allows us to investigate the effect of a varying porosity on filtration in the second stage of the paper, where the steady-state effective equations are applied to a filter whose porosity varies in one direction only and the flow is unidirectional, in the same direction as the porosity variation. We consider a general operating regime where the filter is placed between two reservoirs. We solve the derived system numerically and analyse our results to investigate how a variation in filter porosity affects particle trapping in steady-state operation. We introduce two metrics to quantify the effectiveness of a filter and, through these, we are able to determine how a porosity gradient can improve filter efficiency. Additionally, we perform an asymptotic analysis on the derived system, exploiting the small porosity gradient (a common feature in many industrial operating regimes) to determine very accurate analytical expressions at a much reduced computational cost. These tasks are all carried out in §4. We conclude in §5 with an overview and discussion of the results from this paper, ending with relevant extensions based on this work.

## 2. Model description

We consider the transport of particles via advection and diffusion through a porosity-graded filter. The filter is modelled as a collection of solid obstacles, to which particles can adsorb. We describe the particle locations in terms of the concentration

The concentration field is defined within the fluid phase of the domain, where *d*=2 or 3 the number of spatial dimensions. We define the solid phase of the domain as *Ω*_{s}=*Ω*∖*Ω*_{f}. The solid phase is modelled by a collection of fixed non-overlapping *d*-dimensional balls, whose centres are located on a square or cubic lattice at a distance *δl* apart, where *δ* is a (small) dimensionless parameter and *l* is the characteristic filter length. We allow the radii of the balls to vary in space, and a ball with centre at

A consequence of the assumption that the contaminants are small particles which are in dilute suspension is that particle–particle interactions are negligible, and the dominant interaction is the adsorption of particles on the obstacles. Thus, the particle concentration is governed by the standard advection–diffusion equation with a reactive (or partially adsorbing) Robin boundary condition, which represents a linear adsorption rate:
*D*>0 is the constant diffusion coefficient, ∂*Ω*_{f} is the interface between the solid obstacles and the fluid region (hereafter known as the solid–fluid interface), ** n** is the outward unit normal to ∂

*Ω*

_{f},

*γ*≥0 is the constant particle-adsorption coefficient. There is no adsorption when

*γ*=0, and the adsorption is instantaneous in the limit as

Additionally, we assume that we have an incompressible Newtonian fluid, which satisfies Stokes equations,
*μ* is the fluid viscosity and

### (a) Dimensionless equations

We scale the variables via *O*(1) parameters. As we see in §3, this choice yields a distinguished limit in which all mechanisms are present at leading order in the macroscale. We note that *k* may be expressed in terms of the Péclet number and the Damköhler number *Da* (which relates adsorption to diffusion) via *k*=*Da*/(*δPe*). The parameter *δ* that premultiplies the adsorption term in (2.3b) is required to ensure a dominant balance between particle transport and adsorption, because an *O*(*δ*) adsorption rate over an *O*(1/*δ*) number of obstacles will lead to an overall *O*(1) adsorption rate. Chernyavsky *et al*. [18] show that an asymptotic analysis of a similar advection–reaction–diffusion equation with *Pe*=*O*(1) and *k*=*O*(1/*δ*) leads to a breakdown in the asymptotic series considered. In the limit of large *k*, the solute is adsorbed very quickly and, for this work, could be tracked via an initial boundary layer in time of *O*(1/*k*) and, if necessary, diffusive boundary layers in space of

The dimensionless version of (2.2) is
*d*-cubic lattice of *d*-balls a distance of *δ* apart, and a *d*-ball with centre at ** x** has radius

*δR*(

**). A schematic of the dimensionless geometry is shown in figure 1**

*x**a*.

## 3. Homogenization

In order to reduce the complexity of the problem geometry, we homogenize the governing equations (2.3)–(2.4) using the method of multiple scales. This will allow us to obtain effective equations on a simpler macroscale domain while formally capturing the relevant information about the microscale geometry. As is standard in homogenization theory, we introduce a microscale variable ** y**=

**/**

*x**δ*and treat

**and**

*x***as independent. The extra degree of freedom this allows is removed by imposing that the solution is exactly periodic in**

*y***, and hence any small variation between unit cells is captured through the macroscale variable**

*y***. As shown in figure 1, the microscale variable**

*x***is defined in the unit cell**

*y**ω*(

**), whereas the macroscale variable**

*x***spans the entire filter. The solid portion of the cell, occupied by the obstacle, is denoted by**

*x**ω*

_{s}(

**), and the fluid portion is**

*x**ω*

_{f}(

**)=**

*x**ω*(

**)∖**

*x**ω*

_{s}(

**). The boundary of the unit cell owing to the obstacle is denoted by ∂**

*x**ω*

_{s}(

**), and the outer boundary of the unit cell is ∂**

*x**ω*(

**). Further, we consider each dependent variable as a function of both**

*x***and**

*x***. Using the new variable**

*y***, we transform spatial derivatives in the following manner:**

*y*_{x}and ∇

_{y}refer to the nabla operator in the

**- and**

*x***-coordinate systems, respectively.**

*y*Our goal is to systematically reduce the geometrical complexity of the problem by deriving effective governing equations valid in the macroscale domain. For this purpose, the variables we use in our macroscale equations are quantities averaged over an entire cell *ω*(** x**). To this end, we define the porosity

*ϕ*(

**) to be**

*x**ω*(

**)|=1), and the volumetric average concentration**

*x**C*and volumetric average fluid velocity

**(also known as the Darcy velocity in porous-media formulations) as follows**

*U**c*≡0 and

**≡**

*u***in**

*0**ω*

_{s}(

**). We note that, alternatively, we could have averaged over the fluid portion of the cell (using |**

*x**ω*

_{f}(

**)| in the denominator of (3.3)) rather than over the whole cell; this average, known as the intrinsic average, is most commonly used in the method of volume averaging [26]. This would give the particle concentration in terms of the available fluid space of the membrane, rather than the overall particle concentration within a membrane. As we will see later, the volume average concentration is most convenient in our case because it allows us to directly deduce the effect of macroscopic changes in the porosity.**

*x*We first consider the flow problem (2.4), because it has no dependence on the concentration field, and use the ensuing results in the solute-transport problem (2.3).

### (a) Flow problem

Under the spatial transform (3.1), the flow equations (2.4) become
*ω*_{s} is the obstacle boundary, and ∂*ω* is the outer cell boundary as shown in figure 1. Expanding the flow velocity and pressure in powers of *δ* as follows:
*G*∈{** u**,

*p*}, the flow equations (3.4) at

*O*(

*δ*

^{−1}) are

*p*

_{0}=

*p*

_{0}(

**,**

*x**t*), a direct consequence of our pressure scaling choice in §2a.

The *O*(1) terms in (3.4) are given by
**K** and the vector function ** Π** satisfy the problem

**I**

_{d}is the

*d*-dimensional identity matrix, and the dependence of

**K**and

**on**

*Π***is due to the dependence of the microscale boundary ∂**

*x**ω*

_{s}(

**) on the macroscale variable.**

*x*To obtain the effective flow equation for the Darcy velocity (3.3b), we integrate (3.8a) over *ω*_{f}(** x**) to deduce

*K*(

*ϕ*) is a scalar function defined by

**K**in (3.10b) is a multiple of the identity matrix owing to the symmetry of the cell problem described by (3.9). This would no longer be the case if, instead of a spherical obstacle, we had considered an anisotropic obstacle. Integrating (3.4b) over

*ω*

_{f}, and using the transport theorem (outlined in appendix A) in conjunction with the periodic and no-slip boundary conditions (3.4c)–(3.4d), we obtain the following incompressibility condition for the macroscale flow velocity:

**as expected [27].**

*U*### (b) Solute-transport problem

We now perform a similar homogenization for the solute-transport problem. While the Dirichlet boundary conditions for the flow problem (2.4c) are invariant to the multiple-scales transformation, this is not true for the transformation of the Robin boundary condition (2.3b) in the solute-transport problem. Following [22,19], we address this issue by introducing the relation
*χ*(** x**,

**)=0 defines the solid–fluid interface.**

*y*^{1}Using the gradient transform (3.1), the normal vector

**in (2.3b) becomes**

*n*

*n*_{y}=−

**/∥**

*y***∥ is the outward unit normal on the obstacle boundary ∂**

*y**ω*

_{s}(

**), and**

*x**δ*∇

_{x}

*R*accounts for the macroscale effect of varying obstacle size (see [19] for details).

Under the spatial transforms (3.1) and (3.12), the solute-transport problem (2.3) becomes
*c*(** x**,

**,**

*y**t*) in powers of

*δ*as in (3.5) and equating powers of

*δ*, the

*O*(1) terms are

*c*

_{0}is independent of the microscale, i.e.

*c*

_{0}=

*c*

_{0}(

**,**

*x**t*). The

*O*(

*δ*) terms in (3.13) are

*c*

_{1}, which allows us to write

*c*

_{1}as

**satisfy**

*Γ**n*

_{y,i}is the

*i*th component of the unit vector

*n*_{y}.

The *O*(*δ*^{2}) terms in (3.13) are
*ω*_{f} and applying the boundary conditions (3.18*b*,*c*) gives
*s* denotes the differential element of the obstacle surface ∂*ω*_{s}(** x**). The first two terms on the right-hand side of (3.19) can be simplified using the transport theorem as stated in (A 1), giving

*c*

_{0}=

*c*

_{0}(

**,**

*x**t*), |

*ω*

_{f}(

**)|=**

*x**ϕ*(

**) and (3.16), (3.20) reduces to**

*x***J**

^{T}

_{Γ})

_{ij}=∂

*Γ*

_{j}/∂

*y*

_{i}is the transpose of the Jacobian matrix of

**, the solution to the cell problem (3.17). From now on, we simply write**

*Γ**ϕ*but it should be understood that the porosity will be, in general, a function of

**. In a similar manner, we use**

*x**ω*

_{f}(

*ϕ*) instead of

*ω*

_{f}(

**) henceforth.**

*x*To express (3.21) in terms of the volumetric average concentration *C*(** x**,

*t*) defined in (3.3a), we note that

*δ*. Using this relation and (3.10a), we find that

*δ*, where the effective diffusion coefficient is

*ϕ*(

**)=1−**

*x**V*

_{d}

*R*(

**)**

*x*^{d}for spherical obstacles, where

*V*

_{d}is the volume of the unit ball with dimension

*d*, so

*V*

_{2}=

*π*and

*V*

_{3}=4

*π*/3, and |∂

*ω*

_{s}(

*ϕ*)|=

*dV*

_{d}

*R*(

**)**

*x*^{d−1}. We recall that the Darcy velocity that appears in (3.22a) is given by

**(**

*U***,**

*x**t*)=−

*K*(

*ϕ*)∇

_{x}

*p*, with the effective permeability

*K*defined in (3.10b). As for the permeability

*K*, we note the right-hand side of (3.22b) is a multiple of the identity matrix owing to the symmetry of the cell problem (3.17).

The effective permeability *K*(*ϕ*) and the effective diffusion *ϕ*, determined by the size of the spherical obstacle. We do this numerically using the finite-element software Comsol Multiphysics and find that both *K* and *b* is equivalent to that calculated in [22], which considered the limit of ** U**→

**0**,

*k*→0.

The homogenization procedure has significantly decreased the geometrical complexity of the multiply-connected domain in the full problem (2.3), while only slightly increasing the complexity of the governing equations for the homogenized problem (3.22). This complexity manifests in the variable coefficients of the governing equation; these coefficients reflect the microscale structure of the problem in a systematic manner.

In addition to the increase in computational efficiency, our homogenized model gives insight into the physical effect of a non-uniform porosity. Namely, the homogenized governing equation (3.22a) indicates that a gradient in the porosity enters the equation as a term similar to the advection term owing to the Darcy velocity. In other words, a non-uniform porosity in the membrane gives us an extra degree of freedom with which to either enhance or reduce the effect of a fluid flow on the solute concentration distribution. In §4, we determine how exactly to vary a macroscale porosity gradient to improve the effectiveness of a membrane filter.

## 4. A unidirectionally graded filter

### (a) Model set-up

Here, we use the homogenized equations to understand and quantify the effect of porosity gradients on filter efficiency. We model a filtration process using (3.22) in the steady state, i.e. setting ∂*C*/∂*t*=0, coupled with appropriate boundary conditions at the inlet and outlet that represent two reservoirs. Because we are interested in the interaction between the fluid flow and the porosity gradient in the advection term, we consider a filter that is graded in the same direction as the fluid flow. While it can be difficult to fine-tune filter porosity owing to manufacturing constraints, it is acknowledged that a negative porosity gradient can improve filter efficiency. As discussed in the Introduction, a quantitative understanding of this behaviour is an open question.

Specifically, we consider an industrial process whereby a filter with a *unidirectional* constant porosity gradient separates two reservoirs, and a flow is induced through the filter in the same direction as the gradient in porosity. We are interested in the steady-state operation of this process, and how the porosity gradient affects particle adsorption. We define the direction of porosity gradient as *x*, where *x*∈(0,1) within the filter (because the dimensional characteristic length *l* is chosen to be the length of the filter). Thus, the set-up is similar to that illustrated in figure 1. The upstream is defined for *three-dimensional* domain (and take *d*=3 henceforth) and it is only the variation in porosity that is one-dimensional.

Provided the boundary conditions allow for unidirectionality, (3.10) yields a unidirectional constant macroscale flux ** U**=

*U*

*e*_{x}(where

*e*_{x}is the unit vector in the

*x*-direction), at all points in the filter, regardless of the porosity. We are able to take

*U*=1 without loss of generality by choosing the characteristic velocity scale

*x*-direction, a quantity that is easily measured experimentally. The homogenized governing equation for the concentration

*C*within the filter (3.22a) is then

*ϕ*→1 provides governing equations for the upstream

*C*

^{−}and downstream

*C*

^{+}concentrations,

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*^{−}→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

Equations (4.2), (4.3a) and (4.3f) give
*A*,*B*∈[0,1] are two constants to be determined. These expressions may be substituted into the remaining boundary conditions (4.3b)–(4.3e) to yield two boundary conditions for *C*, thus closing the elliptic problem (4.1),
^{3} grid points. To ensure that our results are experimentally relevant, we make comparisons between simulations that use the same transmembrane pressure difference. This corresponds to modifying the dimensionless parameters *Pe* and *k* between simulations, and is discussed in more detail in appendix B.

We consider a linear porosity model of the form
*ϕ*_{0} and gradients *m*. This corresponds to a uniform filter when *m*=0. As we expect, the concentration *C* decreases with *x* as the contaminants are trapped by the filter medium (figure 3*a*). We also observe that the rate at which the concentration changes in *x* (which is related to the rate of particle removal) decreases as the porosity gradient increases (figure 3*a*). We probe this result further in §4b to explore which filter set-up is most efficient at contaminant removal.

### (b) Filter efficiency

Because the goal of filtration is to maximize particle adsorption, it is helpful to consider the particle adsorption distribution, described by the right-hand side of (4.1), within the filter. To this end, we introduce
*x* in the filter. For the concentration profiles depicted in figure 3*a*, we show the corresponding *Λ* in figure 3*b* for comparison. A natural measure to evaluate filter efficiency is the total contaminant removal by the filter per unit filter area per unit time, given by

We investigate how a porosity variation affects the total particle removal within a filter by calculating *T* as we vary the mean porosity *ϕ*_{0} and gradient *m* in the linear porosity function (4.6). Although *T* is observed to increase with decreasing *ϕ*_{0}, as we would expect, there are diminishing returns as we continue to decrease *ϕ*_{0} (figure 4*a*). However, we observe only a weak dependence on *m* (figure 4*a*) whereby, upon including a porosity gradient, there is an increase in the total particle adsorption *T*. Interestingly, we find that the sign of *m* does not affect *T*, exhibited through the symmetry about *m*=0. While it may appear surprising that the orientation of the filter has no effect on total particle adsorption, this phenomenon arises as a result of the commutativity of the obstacle layout: the concentration after each subsequent obstacle is a proportion of the incoming concentration. Hence, the contaminant removal by the entire filter is a product of all these proportions and, moreover, any permutation of the sinks will result in the same total contaminant removal. If our model included pore blocking as a result of particle adsorption, this symmetry feature would no longer hold.

Thus, while the metric *T* is useful, it cannot account for the superior performance often observed for filters with a negative porosity gradient. To explain this experimental observation, we must investigate *Λ* further. Although long-term blockage is inevitable, a particular aim of an efficient filter is to avoid local blockages while the rest of the filter is still functional, i.e. to have the filter fail ‘all at once’ rather than in one place. Thus, while the total contaminant removal is important, we must also consider the uniformity in which it is removed by the filter. To quantify this, we introduce the metric
*S* is minimized (and zero) for uniform contaminant removal and so we posit that a smaller value of *S* corresponds to a more efficient filter.

In contrast to our analysis of the total particle adsorption, *T*, we find that *S* strongly depends on the orientation of the filter (figure 4*b*). Additionally, a negative porosity gradient *m* yields a concentration profile that is closer to providing spatially uniform contaminant removal than a positive *m*. We note however that the minimum value of *S* for a given mean porosity *ϕ*_{0} (shown as asterisks in figure 4*b*) does not necessarily correspond to choosing the porosity gradient as negative as possible. This feature occurs because a porosity gradient that is too negative will have less particle adsorption near the entrance to the filter than towards the exit, thus also creating an uneven adsorption distribution.

Our analysis therefore indicates that, for a given mean porosity, the porosity gradient has a much larger effect on the distribution of particle adsorption within the filter than it does on the total particle adsorption. This provides a quantitative interpretation for the experimental observation that efficiency in contaminant removal is improved when a filter with a negative porosity gradient is used. As the shape of the concentration profile relates to the problem of localized blocking, we hypothesize that a negative porosity gradient allows for a more uniform contaminant removal. This increases the time until localized blocking becomes a problem, and thus increases filter efficiency.

### (c) Slowly varying porosity

The results we have presented so far provide a very useful quantitative analysis on the role of porosity gradients in filtration. However, because they were obtained numerically, it is still rather cumbersome to perform major parameter sweeps through the space of porosity functions and the dimensionless parameters *Pe* and *k* (although we emphasize that results we presented in §4b are still orders of magnitude cheaper than solving the full equations in a complex geometry).

Here, we take a further simplifying step by exploiting the fact that manufacturing constraints limit the possible variations in porosity. We determine approximate analytical solutions to the equations presented, significantly increasing the computational efficiency and allowing us to gain insight into the concentration profiles we obtain. To implement this, we explore the limit of a weakly varying porosity such that d*ϕ*/d*x*=*O*(*ϵ*), where *ϵ*≪1.

We expand *C* and *ϕ*(*x*) in powers of *ϵ* as follows:
*ϕ*_{0} is the mean porosity, thus constant in *x*. We expand *f*(*ϕ*(*x*)) in the same manner as *ϕ*(*x*), i.e. where the leading-order term is constant. The coefficients are given by
*ϕ*(*x*), *f*(*ϕ*(*x*)) are known. For the linear porosity functions (4.6) considered in §4b, we have *ϵ*=*m*, *ϕ*_{1}(*x*)=*x*−0.5, and the asymptotic limit we consider is equivalent to *m*→0, i.e. close to uniform porosity. Note that the theoretical maximum value of *ϵ* for a linear porosity function is *π*/6≈0.52 in three dimensions. However, the maximum value of *ϵ* will be smaller in practice.

Substituting the asymptotic series (4.11) into the system (4.1), (4.5) and equating powers of *ϵ*, the *O*(1) system is given by

The *O*(*ϵ*) terms in the system (4.1) and (4.5) are given by

In addition to the increase in computational efficiency, this analytical result also allows us to determine closed-form approximations up to *O*(*ϵ*) for the metrics *T* and *S*,

The asymptotic concentration profile (up to two terms) agrees very well with the profile obtained numerically, even when the porosity gradient takes physically extreme values (figure 5). Further, the closed-form nature of (4.18) allows us to produce approximate results for *T* and *S* (shown as dashed black lines in figure 4) at a fraction of the computational cost. As with all the analytic results in this section, we observe excellent agreement with the numerical results.

## 5. Discussion

We have systematically derived a macroscopic model for a porosity-graded filter from microscale information by generalizing standard homogenization theory for near-periodic systems. The result is an advection–diffusion–reaction equation for the solute concentration within the filter as a function of the porosity distribution and operating conditions. In particular, this equation has allowed us to investigate how porosity gradients affect solute trapping and filter efficiency. To this end, we have defined two suitable metrics to quantify the implicit effect of blocking. The first corresponds to the rate of total contaminant removal. The second measures how close to uniform the adsorption of the contaminant is within the filter and is an indication of the propensity of a filter to localized blocking. By performing a computationally efficient parameter sweep, we have been able to determine and quantify that, for a given mean porosity, a porosity gradient has a much larger effect on the second metric than on the first metric. In general, we have found that a filter with a decreasing porosity will be less prone to localized blocking and we have quantified the optimal gradient as a function of the operating conditions. This allows us to understand the experimental observation that a decreasing porosity can lead to a greater filter efficiency; the porosity gradient lowers the risk of localized blocking, while maintaining the rate of total contaminant removal.

The homogenization procedure we used allowed a near-periodic, rather than a strictly periodic, microstructure to be considered. The macroscale variation that near-periodicity allows is a vital feature of a porosity-graded filter. The resulting macroscale equations are computationally cheap to solve, allowing us to efficiently explore the experimental parameter space to quantify filter efficiency.

With regards to an ‘ideal’ filter, we note that it is impossible to globally optimize both of the two metrics (*T* and *S*) that we created. Indeed, we found that porosity functions that increase total adsorption were further from uniform adsorption. Thus, there is a trade-off to be made between maximizing adsorption and minimizing potential blockage issues, an idea supported in the network model for trapping [24]. While the question of determining an ‘optimum’ porosity function will depend on the specifications of the end user, our work provides an inexpensive method to increase filter efficiency given these requirements. However, our model does suggest that, in the absence of additional constraints, it is wiser to use a filter with a negative rather than a positive porosity gradient. While this will not affect the total adsorption, it will provide a more uniform adsorption rate throughout the filter, which should reduce blocking issues.

This work can also be used to investigate filters comprising a series of material layers, each with a different mean porosity. One way to achieve this would be to derive suitable boundary conditions to couple each of the regions, using the method in §4a. Another (simpler) idea would be to use a differentiable approximation of the porosity (e.g. cubic splines defined between the centres of neighbouring layers) in the macroscale model defined in the first part of this paper.

Although we considered a filter whose microscale structure consisted of balls, it is possible to consider an arbitrary shape for the microscale obstacle. However, we chose our particular microstructure to allow an explicit macroscale equation, and thus significant analytical progress to be made, and the same will not hold for the general microstructure case. It is a simple task to combine this work with that of Richardson & Chapman [19], 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. However, as the coefficients within the resulting governing equations would have to be determined at each point in space, the ‘reduced’ equations are still computationally expensive.

One drawback of the method presented here is the restriction to a near-periodic microscale, meaning that materials that vary wildly on the microscale cannot be treated similarly. A viable extension to this work would be to test how robust the assumptions on a near-periodic microscale structure are, and how far away one could get before the macroscopic description we derived broke down. Although the near-periodic assumption is required to derive the macroscale equations, it may be the case that it is only the average properties of the microstructure that must possess this property. This could be tested by performing numerical simulations on the full problem for randomly placed spheres. Similar simulations have been carried out in [22] for diffusing particles only, and it was found that the concentration distribution for a random and square array of circles agree in the high-porosity limit. On physical grounds, we expect the same result for our work in the high-porosity limit. Furthermore, in the same limit, we expect the results for a general microscale obstacle to coincide with the high-porosity limit of the governing equations derived in this paper, using appropriate results for pore volumes and obstacle surface areas.

One feature that is lacking in our model is the dynamic effect of blockage (membrane fouling). As blocking would temporally change the geometry of the problem, the flow and particle transport problems would be fully coupled, and hence difficult to investigate. Despite the lack of explicit blocking in our model, we were able to implicitly infer the long-term effect this would have by quantifying the particle adsorption distribution within the filter. If blocking was also considered, we would expect the results in this paper to provide a quasi-static description of that model. Additionally, the linear trapping rate that we imposed is a generalized approximation of a feature that will depend on the solute solubility and the filter material adsorption. Intuitively, one would expect a Michaelis–Menten-type dependence for the adsorption rate on the number of particles, i.e. approximately linear for a small number of particles and bounded above for a large number of particles. Our approach is valid because of the dilute suspension assumption but, if one is in the same asymptotic limit, it is a trivial task to extend the analysis in this paper to a general nonlinear adsorption rate if required.

Finally, we note that this work has the potential not only to guide filter manufacture and operating conditions, but also to provide assistance to other industries that use functionally graded materials, such as heterogeneous artificial body tissue in tissue engineering [28,29], or graded electrodes in lithium-ion batteries [30]. Furthermore, as the technology to produce functionally graded materials grows [31], the potential experimental parameter space will also increase. Forming appropriate mathematical models and maximizing the analytical progress that can be made will significantly expedite the exploration of this parameter space and result in faster technological growth.

## Data accessibility

The corresponding author may be contacted for the numerical implementations described in this work.

## Author' contributions

M.P.D., I.M.G. and M.B. derived the mathematical model. M.P.D. wrote the paper and performed numerical simulations. I.M.G. and M.B. designed the study and co-drafted the paper. All authors gave final approval for publication.

## Competing interests

We declare we have no competing interests.

## Funding

M.P.D. is supported by an EPSRC Impact Acceleration Account Award (grant no. EP/K503769/1). I.M.G. is supported by a Royal Society Fellowship. M.B. is funded by a Junior Research Fellowship of St Johns College, Oxford. This work was also supported by the 2020 Science programme, which is funded through the EPSRC Cross-Disciplinary Interface Programme (grant no. EP/I017909/1).

## Appendix A. Transport theorem for a varying domain

During the homogenization procedure, we use the transport theorem to evaluate integrals over microscale cells which vary in the macroscale variable ** x**. To do this, we must determine the rate of change of the fluid domain

*ω*

_{f}(

**) with respect to**

*x***.**

*x*We note that *ω*_{f}(** x**) has a fixed outer boundary, ∂

*ω*, and an interior boundary on the surface of the sphere of radius

*R*(

**), ∂**

*x**ω*

_{s}(

**). The rate of change of ∂**

*x**ω*

_{s}(

**) with respect to**

*x***is ∇**

*x*_{x}

*R*. To see this, consider the difference between integrals over the domains

*ω*

_{f}(

**+**

*x**ξ*

*e*_{i}) and

*ω*

_{f}(

**) as**

*x**ξ*→0. The resulting domain of integration is a shell whose thickness is approximately

*ξ*∂

*R*/∂

*x*

_{i}as

*ξ*→0. As we are considering an integral over

*ω*

_{f}(

**)=**

*x**ω*(

**)∖**

*x**ω*

_{s}(

**), the relevant velocity of the interior boundary is −∇**

*x*_{x}

*R*.

Therefore, the transport theorem for the geometry considered in this paper is as follows. Let ** v**(

**,**

*x***,**

*y**t*) be a continuous vector field and periodic in

**on the outer boundary ∂**

*y**ω*(

**). Then**

*x*## Appendix B. Parameter modification

While the two dimensionless parameters *Pe* and *k* are useful for a mathematical analysis of the governing equations, within an experimental set-up it is not these parameters that will, in general, be kept constant as we vary the membrane porosity. When analysing results from our model, it is useful to consider mathematical variations that correspond to experimental variations in order to maximize the experimental relevance of our work.

In practice, the control parameter within a filtration system is typically the transmembrane pressure difference (note that the same calculation applies if the volumetric flux is the control parameter that is used instead). This corresponds to a linear variation in the average velocity in the *x*-direction, *P*, we can solve (3.10a) and (3.10c) in one dimension to deduce that (in dimensional variables)
*ϕ*(*x*).

Therefore, we are able to determine how *ϕ*(*x*). That is, for a given reference function *ϕ*_{ref}(*x*), and parameter values *Pe*[*ϕ*_{ref}] and *k*[*ϕ*_{ref}], we can determine the effect of changing *ϕ*(*x*) on *Pe*[*ϕ*] and *k*[*ϕ*] in the dimensionless model that correspond to the new function *ϕ* while keeping the transmembrane pressure difference constant. Accordingly, when we discuss *S* and *T* in §4b, we take *Pe*[*ϕ*] and *k*[*ϕ*] to vary appropriately, relative to a reference point, defined by *ϕ*_{0}=0.75, *m*=0, *Pe*=3, *k*=1 (illustrated in figure 6).

## Footnotes

↵1 Writing

*R*() in (3.11) assumes that*x**R*is a continuous function of the macroscale, rather than a piecewise constant function evaluated at the centre of the relevant unit cell. This simplifies the subsequent analysis while affecting only the boundary condition (3.13b) at higher orders than we need to consider. As a result, our final leading-order macroscale equation (3.22a) is unchanged by employing this simplification.*x*

- Received July 7, 2015.
- Accepted August 20, 2015.

- © 2015 The Author(s)

Published by the Royal Society. All rights reserved.