1 Introduction
Membrane filters are used in a wide variety of applications to remove particles and undesired impurities of a certain size range from a fluid. Membrane filtration is used in applications as diverse as water purification (Iritani Reference Iritani2013), treatment of radioactive sludge (Daniel et al. Reference Daniel, Billing, Russell, Shimskey, Smith and Peterson2011), various purification processes in the biotech industry (Ho & Zydney Reference Ho and Zydney1999, Reference Ho and Zydney2000; Bolton, LaCasse & Kuriyel Reference Bolton, LaCasse and Kuriyel2006a ; Bolton et al. Reference Bolton, Boesch and Lazzara2006b ), the cleaning of air or other gases (Brown et al. Reference Brown, Levison, Titchener-Hooker and Lye2009), and beer clarification (Van der Sman et al. Reference Van der Sman, Vollebregt, Mepschen and Noordman2012). The type of membrane used depends on the specific application, but an overarching requirement is to have fine control over particle removal from the feed solution, while keeping energy requirements to a minimum. Membrane filters used in microfiltration can have rather varied structure (see, for example, figure 1), but may generally be understood to be porous media, with characteristic pore size, shape, and void fraction determined by the manufacturer. Separation of particles from the feed solution may occur in two basic ways: (i) particles larger than pores cannot pass through pores and hence are sieved out; and (ii) particles smaller than pores may be adsorbed within pores and retained within the membrane. With this in mind, the issue of energy requirements for filtration may be understood: a membrane with tiny pores guarantees removal of all particles suspended in the feed solution (by sieving), but provides extremely high resistance to flow, so that a very large pressure drop is required to filter the fluid within a reasonable time frame. In practice, therefore, it is desirable that adsorption be responsible for a significant proportion of the filtration, so that membranes with larger pores operating at lower pressures can be used.
Figure 1. Magnified membranes with various pore distributions and sizes ((a) is from Apel (Reference Apel2001), (b,c) are from Ho & Zydney (Reference Ho and Zydney1999)). Photographs (b,c) have width
$10~\unicode[STIX]{x03BC}\text{m}$
.
In addition, the system resistance changes significantly during the course of filtration, as the pores of the membrane become fouled with impurities, which are carried by the flow. Filter performance thus ultimately deteriorates, via a combination of mechanisms (alluded to above). (i) Particles larger than the pores cannot pass through the membrane. Assuming that such particles follow streamlines (advection-dominated flow with large-particle Péclet number), they will be deposited on top of pores, blocking them. (ii) Particles smaller than the membrane-pore size are deposited (or adsorbed) within the pores, shrinking the pore diameter and increasing membrane resistance. (iii) Once pores are blocked, other particles can form a cake on top of the membrane, adding additional resistance via another porous layer on top. Mathematical models for all three fouling mechanisms have been proposed, based mostly on empirical laws of how membrane resistance relates to total volume of filtrate processed, or net flow rate through the membrane, in the different fouling regimes (see, for example, Bolton et al. Reference Bolton, LaCasse and Kuriyel2006a ; Meng et al. Reference Meng, Chae, Drews, Kraume, Shin and Yang2009; Daniel et al. Reference Daniel, Billing, Russell, Shimskey, Smith and Peterson2011; Giglia & Straeffer Reference Giglia and Straeffer2012; Van der Sman et al. Reference Van der Sman, Vollebregt, Mepschen and Noordman2012, among many others). In this paper we take a different approach, which accounts for the fluid dynamics through idealized pores of specified geometries, and models from first principles the fouling due to adsorption and sieving.
Various models for filtration and fouling, which attempt to address aspects of the effect of the pores’ size, geometry and distribution within the membrane, have been formulated and examined by researchers to date (e.g. Mochizuki & Zydney Reference Mochizuki and Zydney1993; Polyakov, Maksimov & Polyakov Reference Polyakov, Maksimov and Polyakov1995; Zeman & Zydney Reference Zeman and Zydney1996; Pujar & Zydney Reference Pujar and Zydney1997; Mehta & Zydney Reference Mehta and Zydney2005, Reference Mehta and Zydney2006; Hwang, Liao & Tung Reference Hwang, Liao and Tung2007; Polyakov Reference Polyakov2008; Kanani et al.
Reference Kanani, Fissell, Roy, Dubnisheva, Fleischman and Zydney2010; Zydney Reference Zydney2011; Polyakov & Zydney Reference Polyakov and Zydney2013; Griffiths, Kumar & Stewart Reference Griffiths, Kumar and Stewart2014; Jackson et al.
Reference Jackson, Bakhshayeshi, Zydney, Mehta, van Reis and Kuriyel2014; Dalwadi, Griffiths & Bruna Reference Dalwadi, Griffiths and Bruna2015; Griffiths, Kumar & Stewart Reference Griffiths, Kumar and Stewart2016). Several models have been proposed for describing the internal stenosis of membrane pores by deposition of small particles (so-called ‘standard blocking’). Most such models are based on simplifying assumptions such as uniform deposition of particles on pore walls, and round cylindrical pores that traverse the membrane depth. However, as particles deposit on the pore walls, their concentration decreases along the pore depth, and therefore the deposition rate, which is necessarily proportional to local particle concentration, decreases as the feed passes towards the pore outlet. In other words, particle deposition is greater at the upstream side of the filter (pore inlet) than at the downstream side (pore outlet). Experimentally, rather steep particle deposition profiles across the depth of the membrane have been observed; see Jackson et al. (Reference Jackson, Bakhshayeshi, Zydney, Mehta, van Reis and Kuriyel2014) for recent results. To account for such effects, Polyakov et al. (Reference Polyakov, Maksimov and Polyakov1995) and Polyakov (Reference Polyakov1998) proposed what they call an
$m$
-model, based on the assumption that particles can deposit (uniformly) only over the inlet portion of the pore walls characterized by the parameter
$m$
, the ratio of the length of this portion to the whole pore length. This model was further refined by Polyakov (Reference Polyakov2008), who modified it to account for non-uniform deposition of particles within pores. However, this work still assumes an initially uniform pore profile, and takes no account of additional blockage due to sieving of particles larger than pores. Depth-dependent filtration was also considered by Dalwadi et al. (Reference Dalwadi, Griffiths and Bruna2015), using rather different methods from those we use here. These authors use homogenization theory to model a membrane filter as a layered series of spherical obstacles around which the filtered liquid must flow, and which expand as fouling occurs. Griffiths et al. (Reference Griffiths, Kumar and Stewart2016) also made further contributions to understanding the depth variation, formulating a discrete ‘network’ model that treats a membrane as a series of layers, each of which contains cylindrical channels that may shrink under the action of adsorption (or be blocked from above by deposition of a large particle).
The goal of the present paper is to extend the scope of the work outlined above, deriving a continuum model that accounts for membrane internal geometry, and that allows fouling by both particle sieving and particle adsorption to operate simultaneously. We use first-principles modelling to make general predictions about how pore geometry affects filtration performance of a membrane filter. The paper is laid out as follows: in § 2 we introduce a mathematical model for flow through a single pore of specified geometry. Scenarios where flow is driven by specified pressure drop (§ 2.2) and constant flux (§ B.1) are considered. Some sample solutions, which demonstrate the features of fouling and separation, are presented in § 3. We also discuss which initial pore profile, in the restricted class of linear pore profiles, gives the best filtration performance (in a sense that we will make precise). Finally we conclude in § 4 with a discussion of our model and results in the context of real membrane filters.
We acknowledge, of course, that membrane filtration and fouling is a much more complex process than the model assumptions (both of our work and others’) allow for. It is clear from figure 1 that many membrane filters are porous media of very complex microstructure, containing many interconnected pores, possibly winding and tortuous, with varying cross-section. Such complexity makes detailed modelling very challenging, We nonetheless believe that reduced models, of the type considered here, can play a valuable role in guiding filter design. We return to this issue in § 3.2 in our discussion of how pore profiles within the membrane might be optimized, and in § 4.
2 Darcy flow model of filtration
The modelling throughout this section assumes that the membrane is flat and lies in the
$(Y,Z)$
-plane, with unidirectional Darcy flow through the membrane in the positive
$X$
-direction (so-called ‘dead-end’ filtration). The membrane properties and flow are assumed homogeneous in the
$(Y,Z)$
-plane, but membrane structure may vary internally in the
$X$
-direction (depth-dependent permeability), thus we seek a solution in which properties vary only in
$X$
and in time
$T$
. Our model may be considered as a representation of the average state (averaged across the
$(Y,Z)$
-cross-section) of a real membrane of thickness
$D$
, in which spatial fluctuations in the plane of the membrane are present. Throughout this section we use uppercase fonts to denote dimensional quantities; lowercase fonts, introduced subsequently in §§ 2.2 and B.1, will be dimensionless.
The superficial Darcy velocity
$\boldsymbol{U}=(U(X,T),0,0)$
within the membrane is given in terms of the pressure
$P$
by
$$\begin{eqnarray}U=-\frac{K(X,T)}{\unicode[STIX]{x1D707}}\frac{\unicode[STIX]{x2202}P}{\unicode[STIX]{x2202}X},\quad \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}X}\left(K(X,T)\frac{\unicode[STIX]{x2202}P}{\unicode[STIX]{x2202}X}\right)=0,\quad 0\leqslant X\leqslant D,\end{eqnarray}$$
where
$K(X,T)$
is the membrane permeability at depth
$X$
. We consider two driving mechanisms: (i) constant pressure drop across the membrane specified; and (ii) constant flux through the membrane specified. In the former case the flux will decrease in time as the membrane becomes fouled; in the latter, the pressure drop required to sustain the constant flux will rise as fouling occurs. We will focus primarily on case (i) in this paper, and so assume this in the following model description; our simulations for the constant flux scenario shown later require minor modifications to the theory (relegated to appendix B). With constant pressure drop, the conditions applied at the upstream and downstream membrane surfaces are
$$\begin{eqnarray}P(0,T)=P_{0},\quad P(D,T)=0.\end{eqnarray}$$
The key modelling challenge lies in linking the permeability
$K(X,T)$
to measurable membrane characteristics that evolve in time, in order to obtain a predictive model. In this paper we consider a simple model in which the membrane consists of a series of identical axisymmetric pores of variable radius
$A(X,T)$
, which traverse the entire membrane. While this may seem a poor approximation to some types of membrane, it is in fact a rather good description of a track-etched membrane filter of the type shown in figure 1(a). We further suppose the pores to be arranged in a square repeating lattice, with period
$2W$
. The basic set-up is schematized in figure 2 (described in more detail below): we consider a feed solution laden with particles, some of which are large and, if larger than the pores, will block them (sieving); and some of which are small, and are transported down the pore and may be deposited on its walls (adsorption, also referred to as ‘standard blocking’ in the literature).
Mass conservation shows that the pore velocity,
$U_{p}$
(the cross-sectionally averaged axial velocity within each pore), satisfies
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}(\unicode[STIX]{x03C0}A^{2}U_{p})}{\unicode[STIX]{x2202}X}=0,\end{eqnarray}$$
while Darcy’s law for the averaged superficial velocity
$U$
within the pore plus its period box gives
$$\begin{eqnarray}U=-\frac{\unicode[STIX]{x1D719}_{m}K_{p}(X,T)}{\unicode[STIX]{x1D707}}\frac{\unicode[STIX]{x2202}P}{\unicode[STIX]{x2202}X},\end{eqnarray}$$
where
$\unicode[STIX]{x1D719}_{m}=\unicode[STIX]{x03C0}A(X,T)^{2}/(2W)^{2}$
is the local membrane porosity at depth
$X$
, and
$K_{p}(X,T)=A(X,T)^{2}/8$
is the local permeability of an isolated pore. The pore and superficial velocities are related by
$$\begin{eqnarray}U=\unicode[STIX]{x1D719}_{m}U_{p},\end{eqnarray}$$
by a simple flux-balance argument. This flow model is completed by assumptions on how the membrane permeability changes in time due to fouling by particles, discussed below. The key nomenclature used both here and below is summarized for easy reference in table 1.
Figure 2. Schematic showing the single unit of membrane, assumed repeated in a square lattice. Small particles, at concentration
$C(X,T)$
, which enter pores and deposit within, are indicated, as are large particles, which block the pore inlet.
Table 1. Key nomenclature used in the model.
2.1 Fouling model: particle adsorption and sieving
We consider the effects of fouling by the two primary mechanisms discussed in the Introduction: (i) fouling by pore blocking (sieving of particles too large to pass through membrane pores, which thus deposit on the membrane’s upstream surface); and (ii) adsorption of small particles on pore walls. As noted earlier, although cake formation may be an important fouling mechanism in the late stages of the filtration, it is not considered explicitly in this paper. The ‘blocking’ mechanism that we do model is a necessary preliminary to cake formation, and may be considered as the first step in the caking process. The fouling modelling is similar in spirit to that used in our earlier work on pleated membrane filters (Sanaei et al. Reference Sanaei, Richardson, Witelski and Cummings2016); however, that work focused on the effects of the pleating, and took no account of depth-dependent structure within the membrane (which is the specific focus of the present paper). To model the two distinct fouling mechanisms, we treat the large blocking and small adsorbing particles as separate populations: a bimodal particle size distribution. The actual dimensions (of both particles and pores) that might be relevant will change from one situation to another, but it is reasonable to think of our ‘small’ particles as having diameter of no more than 0.1 times the pore diameter, while the ‘large’ particles are around the same diameter as the pore (or larger). We first discuss the pore-blocking mechanism (i).
2.1.1 Fouling by pore blocking
As noted above, we assume that pores are slender tubes spanning the membrane, of length
$D$
and variable radius
$A(X,T)$
, arranged in a
$2W$
-periodic square lattice (see figure 2). Initially all pores have specified radius,
$A(X,0)=A_{0}(X)<W$
. Where an individual pore (at time
$T$
) is unblocked the total flux through it,
$Q_{u,pore}(T)$
, is given (approximately) by
$$\begin{eqnarray}Q_{u,pore}=-\frac{1}{R_{u}}\frac{\unicode[STIX]{x2202}P}{\unicode[STIX]{x2202}X}\quad \text{where }R_{u}=\frac{8\unicode[STIX]{x1D707}}{\unicode[STIX]{x03C0}A^{4}},\end{eqnarray}$$
and
$R_{u}$
is the pore resistance per unit of the membrane depth. Blocking occurs when a large particle becomes trapped at the entrance to a pore, obstructing the flow. Instead of treating such pores as completely closed, permitting no further flow (as many authors do), we instead assume that the blocking increases the pore’s resistance to flow, so that the flux through it is decreased. We model this effect by adding an extra resistance, characterized by the dimensionless parameter
$\unicode[STIX]{x1D70C}_{b}$
, in series with the resistance
$R_{u}$
of the unblocked pore. The flux through a blocked pore,
$Q_{b,pore}(X,T)$
, is then given by
$$\begin{eqnarray}Q_{b,pore}=-\frac{1}{R_{b}}\frac{\unicode[STIX]{x2202}P}{\unicode[STIX]{x2202}X}\quad \text{where }R_{b}=\frac{8\unicode[STIX]{x1D707}}{\unicode[STIX]{x03C0}W^{4}}\left(\left(\frac{W}{A}\right)^{4}+\unicode[STIX]{x1D70C}_{b}\right).\end{eqnarray}$$
Here,
$R_{b}$
is the resistance per unit length of the blocked pore. The dimensionless parameter
$\unicode[STIX]{x1D70C}_{b}$
characterizes the tightness of the seal formed when a large particle sits over a pore: for large values of
$\unicode[STIX]{x1D70C}_{b}$
pore resistance increases dramatically after blocking (a tight seal, permitting only a small fraction of the original flux through the pore), while for small values resistance is almost unchanged (a poor seal, permitting nearly the same flux as the unblocked pore). We can now relate the number densities of unblocked and blocked pores per unit area,
$N(T)$
and
$N_{0}-N(T)$
respectively (where
$N_{0}=N(0)$
and
$N_{0}(2W)^{2}=1$
), to the superficial Darcy velocity by noting that the flux of fluid per unit area is
$$\begin{eqnarray}N_{0}(2W)^{2}U=N(T)Q_{u,pore}+(N_{0}-N(T))Q_{b,pore},\end{eqnarray}$$
so that, on substituting for
$Q_{u,pore}$
from (2.6) and for
$Q_{b,pore}$
from (2.7) in the above, we obtain
$$\begin{eqnarray}U=-\frac{\unicode[STIX]{x03C0}W^{4}}{8\unicode[STIX]{x1D707}}\frac{\unicode[STIX]{x2202}P}{\unicode[STIX]{x2202}X}\left(\frac{N}{(W/A)^{4}}+\frac{N_{0}-N}{(W/A)^{4}+\unicode[STIX]{x1D70C}_{b}}\right),\end{eqnarray}$$
an expression in which, comparing with (2.1), the depth-dependent membrane permeability
$K(X,T)$
is implicit.
To close the blocking model we require an equation describing the evolution of
$N(T)$
, the instantaneous number of unblocked pores. (We use a one-dimensional model in which quantities vary only in the depth of the membrane,
$X$
, and time,
$T$
, intended to represent a spatial average of a real system in which variation in the plane
$(Y,Z)$
of the membrane may be present. A different (much more computationally intensive) approach would be to model individual particles landing on pores at specific
$(Y,Z)$
locations, which requires stochastic considerations. Such an approach was used by Griffiths et al. (Reference Griffiths, Kumar and Stewart2014).) We assume a pore is blocked whenever a particle from the large-particle population, with radius
$S>A(0,T)$
is advected to the pore entrance (we assume that large particles follow the streamlines and do not interact with each other). If we assume a cumulative large-particle size distribution function
$G(S)$
, giving the number of large particles per unit volume of fluid with radius smaller than
$S$
, then the concentration of particles of size
$S>A(0,T)$
is
$G_{\infty }-G(A)$
(where
$G_{\infty }=\lim _{S\rightarrow \infty }G(S)$
is the total large-particle concentration). The probability that a particular pore become blocked (per unit time) is thus
$(G_{\infty }-G(A))$
multiplied by the flux through the pore,
$Q_{u,pore}$
:
$$\begin{eqnarray}\left[\begin{array}{@{}c@{}}\text{Probability per unit time that}\\ \text{pore of radius }A\text{ is blocked}\end{array}\right]=-\left.\frac{\unicode[STIX]{x03C0}A^{4}}{8\unicode[STIX]{x1D707}}\frac{\unicode[STIX]{x2202}P}{\unicode[STIX]{x2202}X}(G_{\infty }-G(A))\right|_{X=0}.\end{eqnarray}$$
It follows that
$N(T)$
, the number density of unblocked particles per unit area, evolves according to the equation
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}N}{\unicode[STIX]{x2202}T}=\left.N\frac{\unicode[STIX]{x03C0}A^{4}}{8\unicode[STIX]{x1D707}}\frac{\unicode[STIX]{x2202}P}{\unicode[STIX]{x2202}X}(G_{\infty }-G(A))\right|_{X=0}.\end{eqnarray}$$
Note that this model predicts that
$N\rightarrow 0$
as
$T\rightarrow \infty$
, so that eventually all pores will block. This is not unexpected: since blocked pores acquire significant additional resistance, flow will preferentially be diverted to unblocked pores (which admit higher flux), advecting the large blocking particles to those yet-unblocked sites. The model assumes blocking only by large particles, and only at the pore inlet, since terms are evaluated at
$X=0$
. Strictly speaking, a pore will be blocked in this way by any particle larger than its narrowest point, so one could argue that the right-hand side should be evaluated at the value
$X=X^{\ast }(T)$
where
$A(X,T)$
achieves its minimum at each instant. We will see, however, that, due to the adsorption occurring preferentially at the pore inlet, for all parameter sets we consider,
$X^{\ast }(T)\rightarrow 0$
quite quickly. Therefore, we do not anticipate that results would change significantly if we took careful account of this effect.
Another potential deficiency of this blocking model is that any particles from the ‘large-particle’ population that are smaller than the pore inlet will not be captured by the membrane, but simply pass through it. This seems possible for particles with significant inertia passing through a simple, track-etched membrane of the type shown in figure 1(a,b), but for membranes of more complex structure a scenario in which a significant proportion of such ‘large’ particles pass through the membrane is likely both unrealistic and undesirable. For this reason, almost all simulations presented in this paper are for the case where all large particles are bigger than the initial pore inlet radius,
$A(0,0)=A_{0}(0)$
. In this case the cumulative large-particle distribution may be taken as
$$\begin{eqnarray}G(S)=\left\{\begin{array}{@{}ll@{}}G_{\infty } & \text{if }S>A_{0}(0),\\ 0 & \text{if }S\leqslant A_{0}(0);\\ \end{array}\right.\end{eqnarray}$$
all particles from this distribution will be sieved by the membrane. For comparison, we will also show some results with an exponential cumulative large-particle distribution of the form
$$\begin{eqnarray}G(S)=G_{\infty }(1-\text{e}^{-BS}),\end{eqnarray}$$
where
$B^{-1}$
is a characteristic particle size in the feed solution. In this case some of these ‘large’ particles are smaller than the pore inlet radius; as noted above, such particles will escape capture altogether and simply be advected straight through the pore.
2.1.2 Fouling by adsorption
To account also for the effects of membrane fouling by particle adsorption, we must specify how the population of small particles is deposited within pores. As indicated in figure 2 we consider these small particles independently of the large blocking particles discussed above, and track the concentration,
$C(X,T)$
, of small particles, averaged over the pore cross-section, as the feed passes down the pore. In general the small particles are advected and diffuse within the flow (Dechadilok & Deen Reference Dechadilok and Deen2006), and adhere to the wall at a rate proportional to their local concentration. The full advection–diffusion model, with dependence on radial coordinate within the pore as well as axial coordinate
$X$
down the pore, is non-trivial, and the details of its analysis are relegated to appendix A. Here we present just the result that emerges after an asymptotic analysis based on a distinguished limit of the particle Péclet number, and averaging over the pore cross-section: the averaged concentration
$C(X,T)$
of small particles satisfies a simple advection model,
$$\begin{eqnarray}U_{p}\frac{\unicode[STIX]{x2202}C}{\unicode[STIX]{x2202}X}=-\unicode[STIX]{x1D6EC}\frac{C}{A},\end{eqnarray}$$
to be solved subject to specified particle concentration at the inlet,
$$\begin{eqnarray}C(0,T)=C_{0}.\end{eqnarray}$$
The (dimensional) constant
$\unicode[STIX]{x1D6EC}$
is intended to capture the physics of the attraction between particles and pore wall that leads to deposition. More details are provided in appendix A, but (2.14) models, in a crude way, effects such as van der Waals’ interactions between suspended particles and the membrane material, and attractive forces due to electrostatic charge. Inherent in (2.14) is an assumption that all small particles are identical with regard to their deposition dynamics, which may well not be true in practice. In addition, the model assumes that particle adherence to a clean membrane is the same as for a pre-fouled membrane, again a questionable assumption. Nonetheless, for appropriate choices of parameters, we expect our model to provide a reasonable approximation to a real system with sufficiently homogeneous feed solution, and to be quite broadly applicable. The pore radius
$A(X,T)$
shrinks in response to the deposition: we propose
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}A}{\unicode[STIX]{x2202}T}=-\unicode[STIX]{x1D6EC}\unicode[STIX]{x1D6FC}C,\end{eqnarray}$$
for some constant
$\unicode[STIX]{x1D6FC}$
(related to the particle size), which simply assumes that the pore cross-sectional area shrinks at a rate determined by the total area of particles deposited locally at depth
$X$
. The initial pore radius is specified throughout the membrane,
$$\begin{eqnarray}A(X,0)=A_{0}(X).\end{eqnarray}$$
Note that particle deposition in our model is permanent and irreversible. It is possible, however, that in the later stages of filtration, as the pores narrow, shear forces become significant enough to lead to some re-dispersion of particles. Such re-dispersal is beyond the scope of our model, though could be incorporated in a more sophisticated treatment. We observe that purely adsorptive fouling can be retrieved by setting
$N\equiv N_{0}$
in (2.9).
2.2 Scaling and non-dimensionalization
When filtration is driven by a constant pressure drop
$P_{0}$
across the upstream and downstream membrane surfaces, we non-dimensionalize the modified Darcy model (2.2), (2.3), (2.5)–(2.17), using the scalings
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle P=P_{0}p,\quad X=Dx,\quad C=C_{0}c,\quad N=N_{0}n,\quad (A,S)=W(a,s),\quad B=\frac{b}{W},\\ \displaystyle (U,U_{p})=\frac{\unicode[STIX]{x03C0}W^{2}P_{0}}{32\unicode[STIX]{x1D707}D}(u,u_{p}),\quad G=G_{\infty }g(s),\quad T=\frac{8\unicode[STIX]{x1D707}D}{\unicode[STIX]{x03C0}P_{0}W^{4}G_{\infty }}t\end{array}\right\}\end{eqnarray}$$
(time here is non-dimensionalized on the blocking time scale), giving the following dimensionless model for
$u(x,t)$
,
$u_{p}(x,t)$
,
$p(x,t)$
,
$a(x,t)$
,
$c(x,t)$
,
$n(t)$
(dimensionless Darcy velocity, cross-sectionally averaged pore velocity, pressure, pore radius, cross-sectionally averaged particle concentration, and number density of unblocked pores, respectively):
$$\begin{eqnarray}4u=\unicode[STIX]{x03C0}a^{2}u_{p},\end{eqnarray}$$
$$\begin{eqnarray}u=-a^{4}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x}\left(\frac{1-n}{1+\unicode[STIX]{x1D70C}_{b}a^{4}}+n\right),\quad \frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}x}=0,\end{eqnarray}$$
$$\begin{eqnarray}u_{p}\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}x}=-\hat{\unicode[STIX]{x1D706}}\frac{c}{a},\quad \hat{\unicode[STIX]{x1D706}}=\frac{32\unicode[STIX]{x1D6EC}\unicode[STIX]{x1D707}D^{2}}{\unicode[STIX]{x03C0}P_{0}W^{3}},\end{eqnarray}$$
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}a}{\unicode[STIX]{x2202}t}=-\unicode[STIX]{x1D6FD}c,\quad \unicode[STIX]{x1D6FD}=\frac{8\unicode[STIX]{x1D707}D\unicode[STIX]{x1D6EC}\unicode[STIX]{x1D6FC}C_{0}}{\unicode[STIX]{x03C0}P_{0}W^{5}G_{\infty }},\end{eqnarray}$$
$$\begin{eqnarray}\frac{\text{d}n}{\text{d}t}=\left.na^{4}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x}(1-g(a))\right|_{x=0},\end{eqnarray}$$
with boundary and initial conditions
$$\begin{eqnarray}p(0,t)=1,\quad p(1,t)=0,\quad c(0,t)=1,\quad a(x,0)=a_{0}(x),\quad n(0)=1.\end{eqnarray}$$
Equation (2.20) then gives the pressure
$p$
as
$$\begin{eqnarray}p=u\int _{x}^{1}\frac{\text{d}x^{\prime }}{a^{4}\left({\displaystyle \frac{1-n}{1+\unicode[STIX]{x1D70C}_{b}a^{4}}}+n\right)},\end{eqnarray}$$
and the superficial Darcy velocity in terms of the pore radius
$a(x,t)$
as
$$\begin{eqnarray}u=\left(\int _{0}^{1}\frac{\text{d}x^{\prime }}{a^{4}\displaystyle \left(\frac{1-n}{1+\unicode[STIX]{x1D70C}_{b}a^{4}}+n\right)}\right)^{-1},\quad \text{hence}\quad u_{p}=4\left(\unicode[STIX]{x03C0}a^{2}\int _{0}^{1}\frac{\text{d}x^{\prime }}{a^{4}\displaystyle \left(\frac{1-n}{1+\unicode[STIX]{x1D70C}_{b}a^{4}}+n\right)}\right)^{-1}.\end{eqnarray}$$
Substituting in (2.21) and (2.23) our final system reduces to (2.26), plus
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}n}{\text{d}t}=-\!\left.n\left(\int _{0}^{1}\frac{\text{d}x^{\prime }}{a^{4}\displaystyle \left(\frac{1-n}{1+\unicode[STIX]{x1D70C}_{b}a^{4}}+n\right)}\right)^{-1}\left(\frac{1-n}{1+\unicode[STIX]{x1D70C}_{b}a^{4}}+n\right)^{-1}(1-g(a))\right|_{x=0}\!\!, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}x}=-\unicode[STIX]{x1D706}ca\int _{0}^{1}\frac{\text{d}x^{\prime }}{a^{4}\displaystyle \left(\frac{1-n}{1+\unicode[STIX]{x1D70C}_{b}a^{4}}+n\right)},\quad \unicode[STIX]{x1D706}=\frac{8\unicode[STIX]{x1D6EC}\unicode[STIX]{x1D707}D^{2}}{P_{0}W^{3}},\quad \frac{\unicode[STIX]{x2202}a}{\unicode[STIX]{x2202}t}=-\unicode[STIX]{x1D6FD}c,\qquad & \displaystyle\end{eqnarray}$$
with
$\unicode[STIX]{x1D6FD}$
as defined in (2.22) and initial/boundary conditions
$$\begin{eqnarray}n(0)=1,\quad c(0,t)=1,\quad a(x,0)=a_{0}(x)<1.\end{eqnarray}$$
In line with (2.12) and (2.13) the two forms considered for the cumulative particle distribution function
$g(s)$
are
$$\begin{eqnarray}g(s)=\left\{\begin{array}{@{}ll@{}}1 & \text{if }s>a_{0}(0),\\ 0 & \text{if }s\leqslant a_{0}(0),\end{array}\right.\end{eqnarray}$$
and
$$\begin{eqnarray}g(s)=1-\text{e}^{-bs},\end{eqnarray}$$
where
$b=BW$
characterizes the ratio of characteristic membrane-pore size to typical particle size. Note that for
$g(s)$
as specified in (2.30), particles are larger than pores throughout, and
$g(a)|_{x=0}=0$
in (2.27).
3 Results
In this section we present some simulations of the model (2.26)–(2.29) described in § 2 above, paying particular attention to how results depend on the pore geometry. Selected results for the case driven by constant flux, equations (B 3)–(B 5), will be presented later in appendix B.
Table 2. Dimensional parameter values (A. Kumar, Private Communication, 2014 and Giglia & Straeffer (Reference Giglia and Straeffer2012)).
Table 3. Dimensionless parameters and approximate values.
Our model contains several dimensionless parameters:
$\unicode[STIX]{x1D706}$
which captures the physics of the attraction between particles and the pore wall; the ratio
$\unicode[STIX]{x1D70C}_{b}$
of the additional resistance due to pore blocking to the original resistance of the unblocked pore; and the dimensionless pore shrinkage rate
$\unicode[STIX]{x1D6FD}$
. For the case in which we consider a distribution of large-particle sizes, with cumulative particle size distribution
$g(s)$
specified by (2.31), we also need to specify the ratio
$b$
of pore size to characteristic particle size in
$g(s)$
. The values of each of these dimensionless quantities depend on physical dimensional parameters that must be measured for the particular system under investigation, and we lack such detailed experimental data; hence we have to make our best estimate as to the most appropriate values to use in our simulations. The parameters are summarized in tables 2 (dimensional parameters) and 3 (dimensionless parameters) along with typical values, where known. Considerable variation is possible from one system to another, as noted in the tables; nonetheless we believe that our simulations illustrate the predictive potential of our model if detailed data are available. Most of the parameters in table 2 (such as
$W$
,
$D$
,
$B$
,
$\unicode[STIX]{x1D6FC}$
,
$G_{\infty }$
,
$N_{0}$
,
$P_{0}$
and
$C_{0}$
) depend on physical characteristics of the filter membrane and the feed fluid, therefore in principle could be measured directly or obtained from the manufacturer. Other parameters are harder to measure directly, but indirect methods can be useful. For example, as noted earlier, the particle–wall attraction coefficient
$\unicode[STIX]{x1D6EC}$
may be estimated by comparing solutions of (2.14) to experiments that reveal the density of particles absorbed within the filter (as obtained by, for example, Jackson et al. (Reference Jackson, Bakhshayeshi, Zydney, Mehta, van Reis and Kuriyel2014) via fluorescence microscopy); but such experiments are non-trivial. Direct estimation of
$\unicode[STIX]{x1D70C}_{b}$
for a real membrane system could be more problematic, but if this is the only model parameter for which no data can be inferred then this may be viewed as an overall fitting parameter.
Given the number of parameters, most of them will be fixed throughout our simulations. The value of the dimensionless attraction coefficient between pore wall and particles,
$\unicode[STIX]{x1D706}$
, is unknown, and could certainly vary widely from one system to another depending on the detailed structure of the filter membrane and on the nature of the feed solution. In the absence of firm data we take
$\unicode[STIX]{x1D706}=2$
for most simulations. The dimensionless pore shrinkage rate,
$\unicode[STIX]{x1D6FD}$
, is unknown but will normally be small (it represents the time scale on which pores close due to adsorption, relative to that on which particles block individual pores from upstream): we set
$\unicode[STIX]{x1D6FD}=0.1$
. We note here that the parameters
$\unicode[STIX]{x1D706}$
and
$\unicode[STIX]{x1D6FD}$
are not independent: if we wish to consider the effect of changing membrane thickness
$D$
for example, we must keep
$\unicode[STIX]{x1D706}\propto \unicode[STIX]{x1D6FD}^{2}$
, while if we wish to investigate the effect of the attraction coefficient
$\unicode[STIX]{x1D6EC}$
we must keep
$\unicode[STIX]{x1D706}\propto \unicode[STIX]{x1D6FD}$
. We expand upon this point later. In the absence of definitive data on
$\unicode[STIX]{x1D70C}_{b}$
, assuming that blocking of a pore by a particle increases its resistance by twice the original resistance of the unblocked pore, we set
$\unicode[STIX]{x1D70C}_{b}=2$
for most simulations. Finally, for those simulations where we allow an exponential distribution of large-particle sizes, with
$g(s)=1-\text{e}^{-bs}$
(see (2.13) and (2.31)), we consider values of
$b$
in the range 0.2–10. We briefly demonstrate the effect of changing certain key parameters in figures 4 and 5, below.
As noted in the Introduction § 1, most prior work considers a uniform initial pore profile. While we cannot consider all possible initial pore profiles, we present results for a selection of profiles that allow us to model uniform, increasing, decreasing and non-monotone membrane resistances, as functions of depth. According to the Darcy model, the local membrane resistance is proportional to
$A(X,0)^{-4}$
. Using the non-dimensionalization of (2.18), one can define a dimensionless averaged membrane resistance,
$r(t)$
, as
$$\begin{eqnarray}r(t)=\int _{0}^{1}\frac{\text{d}x}{a(x,t)^{4}}.\end{eqnarray}$$
In order to make a meaningful comparison, we run simulations for pore shapes that give the same initial membrane resistance
$r_{0}=r(0)$
. This means that we are comparing membranes that perform identically when no fouling occurs – they would give identical throughputs when filtering pure water under the same applied pressure drop. The pore profiles considered are:
$$\begin{eqnarray}a(x,0)=\left\{\begin{array}{@{}ll@{}}a_{1}(x,0)=0.904 & \text{uniform initial pore profile,}\\ a_{2}(x,0)=0.16x+0.83 & \text{linear increasing initial pore profile,}\\ a_{3}(x,0)=0.99-0.16x & \text{linear decreasing initial pore profile,}\\ a_{4}(x,0)=0.874+0.39(x-0.5)^{2} & \text{convex parabolic initial pore profile,}\\ a_{5}(x,0)=0.933-0.33(x-0.5)^{2} & \text{concave parabolic initial pore profile,}\end{array}\right.\end{eqnarray}$$
all of which correspond to the same initial membrane resistance
$r_{0}=1.50$
. We note that these initial profiles also happen to have very similar average porosity or void fraction,
$\unicode[STIX]{x1D719}_{0}=(\unicode[STIX]{x03C0}/4)\int _{0}^{1}a(x,0)^{2}\,\text{d}x$
(though this would not be true of all equal-resistance membranes):
$\unicode[STIX]{x1D719}_{0}=0.64$
for
$a_{1}$
and
$a_{5}$
, and
$\unicode[STIX]{x1D719}_{0}=0.65$
for
$a_{2}$
,
$a_{3}$
and
$a_{4}$
. Three of the profiles (
$a_{1}$
,
$a_{4}$
and
$a_{5}$
) are initially symmetric about the membrane centreline, but we will see that in all cases asymmetry rapidly develops due to the particle adsorption within pores.
We solve the model numerically for each chosen pore profile, until the membrane becomes impermeable and the total flux through it falls to zero at final time
$t=t_{f}$
(when the pore radius
$a\rightarrow 0$
). Our numerical scheme is straightforward, based on first-order accurate finite difference spatial discretization of the equations, with a simple implicit time step in the pore-blocking equation (2.22) and trapezoidal quadrature to find the integrals in (2.26).
Figure 3. Filtration simulations: (a) the cross-sectionally averaged pore velocity
$u_{p}$
with uniform initial pore profile
$a_{1}(x,0)=0.904$
; (b–f) the pore radius and particle concentration at selected times up to the final blocking time (
$t_{f}$
, indicated in the legends) for different initial pore radius profiles: (b)
$a_{1}(x,0)=0.904$
, (c)
$a_{2}(x,0)=0.16x+0.83$
, (d)
$a_{3}(x,0)=0.99-0.16x$
, (e)
$a_{4}(x,0)=0.874+0.39(x-0.5)^{2}$
, (f)
$a_{5}(x,0)=0.933-0.33(x-0.5)^{2}$
; (g) total flux versus throughput for these initial profiles for homogeneous ((2.30), black curves) and exponential ((2.31), red curves) distributions of large particles, with
$\unicode[STIX]{x1D706}=2$
,
$\unicode[STIX]{x1D6FD}=0.1$
,
$\unicode[STIX]{x1D70C}_{b}=2$
and
$b=0.5$
.
3.1 Model simulations
We present results for the model summarized in § 2.2 according to the scenarios discussed above. The main results are shown in figure 3: we simulate the model for each of the initial profiles given in (3.2), with parameters
$\unicode[STIX]{x1D706}=2$
and
$\unicode[STIX]{x1D70C}_{b}=2$
characterizing the effects of fouling by adsorption and pore blocking (sieving). Figure 3(a) shows the cross-sectionally averaged pore velocity
$u_{p}$
for the initially uniform pore profile, and figure 3(b–f) shows the pore radius
$a(x,t)$
and the concentration of small particles
$c(x,t)$
for each of the five different initial pore profiles, at various times throughout the evolution. Other parameter values are given in the figure caption. The cumulative large-particle size distribution function is as given in (2.30), so that all of the large particles in the feed are bigger than the pore inlet size and are therefore sieved out. A striking feature of these plots is that pore closure (accompanied by cessation of filtration) occurs first at the upstream membrane surface, even for pores that are initially widest on that side. This is consistent with the graph of the pore velocity
$u_{p}$
(figure 3
a), which is initially uniform (for the initially uniform channel) but rapidly becomes non-uniform, becoming much higher at the narrowing pore inlet; and also with the particle concentration graphs, which show that most of the particle deposition occurs at the pore inlet. This effect becomes more pronounced at later times as the pore radius shrinks near the inlet, further enhancing the deposition there. The graphs of
$c(x,t)$
in figure 3(b–f) demonstrate that the filter membranes are initially capturing more than 90 % of small particles (by adsorption) in all cases, with this proportion increasing to nearly 100 % at later times. The capture proportion could be adjusted by varying the parameter
$\unicode[STIX]{x1D706}$
: increasing (decreasing)
$\unicode[STIX]{x1D706}$
will increase (decrease) the proportion of particles captured. The effect of
$\unicode[STIX]{x1D706}$
is discussed further below.
Figure 4. Flux–throughput graphs for the uniform initial pore profile
$a_{1}(x,0)=0.904$
, with
$\unicode[STIX]{x1D706}=2$
and
$\unicode[STIX]{x1D6FD}=0.1$
, (a) for several different values of
$\unicode[STIX]{x1D70C}_{b}$
(relative increase in pore resistance when blocked by a large particle) with homogeneous particle size distribution
$g(s)$
given by (2.30); and (b) for both homogeneous large-particle distribution (2.30), and for several non-homogeneous particle distributions given by (2.31) with different values of
$b$
(a measure of the relative sizes of pores and particles), with
$\unicode[STIX]{x1D70C}_{b}=2$
.
Figure 5. (a) Particle concentration at the downstream membrane surface
$x=1$
(pore exit) versus throughput, with
$\unicode[STIX]{x1D706}=1$
,
$\unicode[STIX]{x1D6FD}=0.05$
,
$\unicode[STIX]{x1D70C}_{b}=2$
and
$g(s)$
given by (2.30), for several initial pore profiles given in (3.2). (b) Particle concentration versus throughput graph for several different values of
$\unicode[STIX]{x1D706}$
, for the uniform initial pore profile
$a_{1}(x,0)=0.904$
, with
$\unicode[STIX]{x1D70C}_{b}=2$
and
$g(s)=1$
. For the black curves we set
$\unicode[STIX]{x1D6FD}\propto \unicode[STIX]{x1D706}$
(corresponding to varying
$\unicode[STIX]{x1D6EC}$
in (2.14)) and for the red curves
$\unicode[STIX]{x1D6FD}\propto \sqrt{\unicode[STIX]{x1D706}}$
(corresponding to varying membrane thickness
$D$
).
A common experimental characterization of membrane filtration performance is the graph of total flux through the membrane at any given time versus the total volume of filtrate processed at that time (throughput); the so-called flux–throughput graph for the membrane. Since the flux is directly proportional to the averaged Darcy velocity, we define our dimensionless flux by
$q(t)=u(0,t)$
; throughput is then defined by
$\int _{0}^{t}q(t^{\prime })\,\text{d}t^{\prime }$
. We plot these curves, for each of the five pore profiles considered, in figure 3(g). This plot also shows the equivalent flux–throughput graphs for the exponential large-particle distribution function of (2.31), for comparison with the homogeneous case (2.30). The graphs collectively demonstrate that, although all pore profiles give the same initial average membrane resistance (and, to a good approximation, the same initial porosity or void fraction), they exhibit significant differences in performance over time. In particular, membranes whose pores are widest on the upstream side give notably better performance overall according to this performance measure, with more filtrate processed under the same conditions. The membrane with least total throughput is that whose pores are initially narrowest on the upstream side, exhibiting rapid pore closure (pore profile
$a_{2}(x)$
in (3.2)). Furthermore, we see that the flux–throughput curves are initially concave, becoming convex only as total system blockage is approached. This change in curvature has been observed in experimental systems (e.g. Giglia & Straeffer (Reference Giglia and Straeffer2012)) but rarely in model simulations: it seems that only models that incorporate multiple blocking mechanisms simultaneously can exhibit such behaviour. It appears to be indicative of the different blocking regimes: in the early stages the pore blocking is the dominant mechanism responsible for the decrease in flux, while in the latter stages adsorptive blocking dominates (at least for the choice of parameters used here). The differences in performance noted here for homogeneous (2.30) and exponential (2.31) distributions of large-particle sizes are in part due to the fact that for the homogeneous distribution all large particles are sieved, while in the exponential case some are smaller than the pore and pass through the membrane. System resistance therefore increases more rapidly in the former case.
Figure 4(a) demonstrates the effect of varying the parameter
$\unicode[STIX]{x1D70C}_{b}$
, which measures the relative increase in pore resistance when a pore is blocked by a large particle. All of these simulations are performed for the homogeneous large-particle size distribution (2.30). We note that as
$\unicode[STIX]{x1D70C}_{b}$
varies from large to small there is a clear qualitative change in the shape of the flux–throughput performance curves, as the model transitions from blocking-dominated to adsorption-dominated. A large value of
$\unicode[STIX]{x1D70C}_{b}$
means that blocking of a pore by a large particle leads to a tight seal at the inlet and a significant increase in resistance (blocking-dominated); while a small value means that a blocked pore is very loosely sealed, and offers only marginally more resistance to flow than an unblocked pore (hence fouling will be adsorption-dominated). Again, see the experimental data in, for example Giglia & Straeffer (Reference Giglia and Straeffer2012), which reveal similar features as the membrane type and/or filtrate is varied. Since increasing the value of
$\unicode[STIX]{x1D70C}_{b}$
adds more total system resistance, the total throughput decreases monotonically as
$\unicode[STIX]{x1D70C}_{b}$
increases.
Similar qualitative changes may be observed in the case where we model a feed solution with a distribution of particle sizes, described by the exponential size distribution function
$g(s)=1-\text{e}^{-bs}$
(from (2.31)). In this case, whether or not particles are sieved (and hence block pores) depends on their size relative to the instantaneous pore radius. As the parameter
$b$
, measuring the relative characteristic sizes of pores and particles, is varied (figure 4
b), we again see the transition from blocking-dominated to adsorption-dominated behaviour: when
$b$
is small, pores are smaller than particles, so pores rapidly become blocked and this is the dominant fouling mode; whereas when
$b$
is large pores are larger than particles, hence little sieving occurs, and adsorption is the dominant fouling mode. A secondary consequence of increasing
$b$
(‘large’ particles becoming increasingly small relative to pores) is that the fouling is slower and net system resistance increases more slowly (large particles that are not sieved pass straight through the membrane), so total throughput increases. It must be remembered, however, that this increased throughput is achieved at the expense of decreased particle removal. For comparison, figure 4(b) also shows the result for the case when the feed solution contains a homogeneous distribution of large particles that are larger than the pores (
$g(s)=0$
throughout). As anticipated, this case is close to the small-
$b$
simulations, where most particles are larger than pores.
Another key consideration in evaluating membrane performance is the concentration of particles that remain in the filtrate as it exits the membrane,
$c(1,t)$
: in general a lower particle concentration at the outflow side of the membrane indicates superior separation efficiency for the filter membrane. Figure 5(a) plots
$c(1,t)$
versus throughput for each of the initial profiles given in (3.2). The results here are qualitatively consistent with those of the flux–throughput graphs of figure 3(g), in particular, for a given ‘tolerance’ value of the particle concentration at the outlet, membranes with narrow pores on the upstream side always give less total throughput than those whose pores are wide on the upstream side. Note, however, that in figure 5(a), in order to obtain sufficiently distinct graphs, we set the dimensionless membrane–pore attraction coefficient
$\unicode[STIX]{x1D706}=1$
and the dimensionless pore shrinkage rate
$\unicode[STIX]{x1D6FD}=0.05$
, while values
$\unicode[STIX]{x1D706}=2$
and
$\unicode[STIX]{x1D6FD}=0.1$
were used in figure 3.
It is also of interest to study the influence of the membrane–pore attraction coefficient
$\unicode[STIX]{x1D6EC}$
, and of characteristics such as the membrane thickness
$D$
. These each appear in two of our dimensionless parameters:
$\unicode[STIX]{x1D706}=8\unicode[STIX]{x1D6EC}\unicode[STIX]{x1D707}D^{2}/(P_{0}W)$
, and
$\unicode[STIX]{x1D6FD}=8\unicode[STIX]{x1D707}D\unicode[STIX]{x1D6EC}\unicode[STIX]{x1D6FC}C_{0}/(\unicode[STIX]{x03C0}P_{0}W^{5}G_{\infty })$
(see table 3); hence we cannot vary the parameters
$\unicode[STIX]{x1D706}$
and
$\unicode[STIX]{x1D6FD}$
in isolation. To study the effect of the attraction coefficient
$\unicode[STIX]{x1D6EC}$
(see (2.14); this would be changed by, for example, changing the membrane material, or the type of particles in the feed solution) we consider different values of
$\unicode[STIX]{x1D706}$
, with
$\unicode[STIX]{x1D6FD}$
changed in proportion to
$\unicode[STIX]{x1D706}$
, consistent with the way
$\unicode[STIX]{x1D6EC}$
appears in the definitions of these two parameters. To study the effect of changing the membrane thickness
$D$
, we again consider different values of
$\unicode[STIX]{x1D706}$
, but now take
$\unicode[STIX]{x1D6FD}\propto \sqrt{\unicode[STIX]{x1D706}}$
. In figure 5(b) we plot the particle concentration at the pore outlet,
$c(1,t)$
, versus throughput for several different values of
$\unicode[STIX]{x1D706}$
, with
$\unicode[STIX]{x1D6FD}\propto \unicode[STIX]{x1D706}$
or
$\unicode[STIX]{x1D6FD}\propto \sqrt{\unicode[STIX]{x1D706}}$
, corresponding to these two distinct system changes.
In the former case (
$\unicode[STIX]{x1D6FD}\propto \unicode[STIX]{x1D706}$
) we associate small values of
$\unicode[STIX]{x1D706}$
with weaker membrane–particle attraction. As anticipated, this is observed to give rise to poor separation of particles from feed, with a significant fraction of the small particles remaining suspended in the flow at the pore outlet. Large values of
$\unicode[STIX]{x1D706}$
, corresponding to strong membrane–particle attraction, give uniformly low particle concentrations at the outlet. Such strong attraction is, of course, associated with faster total blocking of the membrane pores: if all particles adhere to the pore wall then the pore will close sooner. In the latter case (
$\unicode[STIX]{x1D6FD}\propto \sqrt{\unicode[STIX]{x1D706}}$
), we associate smaller/larger values of
$\unicode[STIX]{x1D706}$
with thinner/thicker membranes. As above, we expect worse/better separation in this case (as measured by the particle concentration at the pore outlet), and this is borne out in the simulations. With a thin membrane the feed solution transits too quickly to deposit all of its particles (though the flux–throughput characteristics would look favourable, since a thin membrane affords little resistance to the flow, and poor particle removal corresponds to a slow fouling process). With a thick membrane the feed remains within the membrane long enough to deposit nearly all its particles, so separation is good; but of course the tradeoff is poorer flux–throughput characteristics, the thick membrane providing higher flow resistance and the good separation leading to more rapid fouling. Interestingly though, in comparing these two scenarios, changing the attraction coefficient appears to have a much larger effect on the overall behaviour than does changing the membrane thickness, suggesting that membrane chemistry could be a very important consideration affecting overall performance.
In the context of these observations we emphasize that our model considers all small particles to be identical. It may be the case that the feed contains several populations of small particles, each with a different deposition coefficient. While we do not explicitly model such a scenario, it would be a fairly straightforward extension to our model, and could be useful for an application in which only certain species are to be removed from a feed solution. We note also that
$\unicode[STIX]{x1D706}$
and
$\unicode[STIX]{x1D6FD}$
are assumed to be constant throughout the duration of filtration. In reality it may well be the case that particles adhere differently to the clean membrane than to the fouled membrane, so that the values of
$\unicode[STIX]{x1D706}$
and
$\unicode[STIX]{x1D6FD}$
should change as filtration progresses. Such considerations are beyond the scope of this paper.
Figure 6. (a) Total throughput
$\int _{0}^{t}q(t^{\prime })\,\text{d}t^{\prime }$
versus initial pore gradient
$b_{0}$
, keeping dimensionless initial net resistance,
$r_{0}=\int _{0}^{1}(\text{d}x^{\prime }/(a_{0}+b_{0}x^{\prime })^{4})=5$
, fixed, for several different values of
$\unicode[STIX]{x1D706}$
, with homogeneous distribution of large-particle sizes (2.30) and
$\unicode[STIX]{x1D70C}_{b}=2$
. For the thin curves we set
$\unicode[STIX]{x1D6FD}\propto \unicode[STIX]{x1D706}$
(corresponding to varying
$\unicode[STIX]{x1D6EC}$
; maximum throughput marked by red dots) and for the thick curves
$\unicode[STIX]{x1D6FD}\propto \sqrt{\unicode[STIX]{x1D706}}$
(corresponding to varying membrane thickness
$D$
; maximum throughput marked by blue diamonds). (b) Total throughput versus initial pore gradient
$b_{0}$
, for several different values of dimensionless initial resistance
$r_{0}$
with homogeneous distribution of large-particle sizes (2.30),
$\unicode[STIX]{x1D706}=2$
,
$\unicode[STIX]{x1D6FD}=0.1$
and
$\unicode[STIX]{x1D70C}_{b}=2$
. The red dots are maximum throughput for each given initial resistance
$r_{0}$
.
3.2 Optimal initial membrane-pore profile
One question of interest to manufacturers is: for a membrane of given net (average) resistance, what is the optimum porosity profile as a function of depth through the membrane? For our model this translates to: what is the optimal shape of the filter pores among all filters with the same initial average resistance? To answer this question, we must first decide how to define filtration performance. This definition will vary depending on the user requirements (we have already seen above the tradeoff inherent between maximizing throughput and simultaneously removing as many particles as possible from the feed), but for purposes of illustration we will use the common experimental characterization of performance as the total throughput over the filter lifetime, as introduced earlier, noting that our methods can easily be adapted to give predictions for any other chosen efficiency measure. In our exposition below we consider optimizing performance while fixing the initial average membrane resistance.
Since the general optimization problem is very challenging, requiring consideration of pores of all possible shapes, we simplify by restricting attention to the class of membranes with pores whose initial radius
$A(X,0)$
varies linearly with membrane depth
$X$
. In order to make a meaningful comparison, we consider members of the family of all linear initial pore profiles,
$a(x,0)=a_{0}+b_{0}x$
, with the same initial resistance
$r_{0}=r(0)$
(as defined by (3.1)). For a given value of the pore profile gradient
$b_{0}$
, the intercept value
$a_{0}$
is then fixed. Note that, depending on the chosen value of
$r_{0}$
, not all values of
$b_{0}$
may be possible: for a low-resistance, highly permeable membrane, the pore occupies a large fraction of the period box (within which it must be entirely confined), and hence the range of values of
$b_{0}$
will be limited in such cases to small absolute values.
Figure 7. Initial total particle concentration at the pore outlet versus initial pore gradient
$b_{0}$
, for several different (small) values of dimensionless initial resistance
$r_{0}$
with homogeneous distribution of large-particle sizes (2.30),
$\unicode[STIX]{x1D706}=2$
,
$\unicode[STIX]{x1D6FD}=0.1$
and
$\unicode[STIX]{x1D70C}_{b}=2$
.
Figure 6(a,b) illustrates our results, plotting throughput versus pore gradient for several different scenarios. In figure 6(a), the dimensionless initial resistance is fixed at
$r_{0}=5$
(a value chosen large enough such that a wide range of pore gradients are available), and throughput is plotted as a function of pore gradient for several different values of the deposition coefficient
$\unicode[STIX]{x1D706}$
. Recalling the discussion at the end of § 3.1 above, we cannot change
$\unicode[STIX]{x1D706}$
in isolation; here we consider the two cases discussed there: (i) we change
$\unicode[STIX]{x1D6FD}$
proportionally to
$\unicode[STIX]{x1D706}$
, modelling changes in the dimensional particle–membrane attraction coefficient
$\unicode[STIX]{x1D6EC}$
(the thin curves), and (ii) we change
$\unicode[STIX]{x1D6FD}$
proportionally to
$\sqrt{\unicode[STIX]{x1D706}}$
, modelling changes in the membrane thickness
$D$
(the thick curves). (We note that there are many ways in which this study could be extended to give a more general optimization: we could for example allow both pore gradient and membrane thickness, and/or the size of the period box, to vary simultaneously while keeping initial membrane resistance fixed, which would involve a sweep through a larger parameter space. However, viewing the present work as a preliminary study, we defer a more general investigation to a future publication.) In figure 6(b) total throughput is again plotted versus slope of the initial pore profile, for several different values of the membrane resistance
$r_{0}$
. As noted above, only a limited range of pore gradients are realizable at low resistances.
In all cases shown in figure 6, the optimum (as measured by maximal total throughput) is achieved at the most negative value of the pore profile gradient – in other words, the pore profile giving maximal total throughput is always that which is as wide as possible at the upstream membrane surface. This result is perhaps unsurprising given our previous simulation results of figure 3 showing the rapid pore closure at the upstream surface: maximum throughput will be achieved by delaying this closure as long as possible.
Note that figure 6 tells us nothing about the proportion of small particles captured in each filtration scenario (though they assume capture by sieving of all large particles). In figure 7, we plot the (initial) total small-particle concentration at the pore outlet, for each case shown in figure 6(b). It is evident from figure 7 that (i) the net capture of small particles depends only weakly on the pore gradient; and (ii) in any case, the most negative pore gradient is favourable to improved net particle capture, leading to a slightly lower concentration of particles at the outlet compared with most other pore gradients.
If we set such concerns aside then we might suspect that a different optimal result would be obtained for very low values of the deposition coefficient, where pore closure might in fact occur at an internal point for pores of decreasing radius. Figure 8 confirms this expectation: for
$\unicode[STIX]{x1D706}=0.01$
and
$\unicode[STIX]{x1D6FD}=0.1$
the optimum profile is no longer the widest possible at the upstream side. A more uniform profile is now favoured, but it should be noted that (i) gains in total throughput are only marginal in this situation; and (ii) in such small-
$\unicode[STIX]{x1D706}$
simulations, only a very small fraction of the small particles is removed by the membrane.
Though we show results here only for the class of linear pore profiles, we note that preliminary investigations of other classes of pore shape (quadratic, cubic and simple exponential profiles) suggest very similar findings. In particular, although the actual optimal pore shapes obtained are somewhat different in each case, the maximal throughput in all cases is very similar, with only marginal improvements over the linear case.
Figure 8. Total throughput versus initial pore gradient
$b_{0}$
, for several different (small) values of dimensionless initial resistance
$r_{0}$
with homogeneous distribution of large-particle sizes (2.30),
$\unicode[STIX]{x1D706}=0.01$
,
$\unicode[STIX]{x1D6FD}=0.1$
and
$\unicode[STIX]{x1D70C}_{b}=2$
. The red dots are maximum throughput for each case: note that in contrast to figure 6 maximal throughput now occurs at some intermediate value of the pore gradient, and not the most negative value.
4 Conclusions
We have presented a model that can describe the key effects of membrane morphology on separation efficiency and fouling of a membrane filter. Our model accounts for Darcy flow through the membrane, and for fouling by two distinct mechanisms: pore blocking (sieving) by large particles, and adsorption of small particles within pores. While essentially predictive, our model contains several parameters that may be difficult to measure for a given system – most notably, the relative increase in pore resistance due to a blocking event,
$\unicode[STIX]{x1D70C}_{b}$
; the dimensionless attraction coefficient between the membrane-pore wall and particles,
$\unicode[STIX]{x1D706}$
; and the dimensionless pore shrinkage rate,
$\unicode[STIX]{x1D6FD}$
. In practice such parameters could be inferred by fitting to a reliable dataset; but even so these parameters will vary from one membrane-feed system to another, since they depend on membrane structure, on the type and size of the particles carried by the feed, and on the chemical interactions between the particles in the feed and the membrane material. The model as presented here implicitly assumes that all of the ‘small’ particles comprising the concentration
$C$
are identical, but it would not (we think) be difficult to extend the modelling to account for several different types of small particles, each with its own concentration and its own sticking parameter.
In the absence of firm data on model parameters we have chosen what we believe to be plausible parameter values (summarized at the start of § 3) for most of our simulations. The focus in this paper is on development of a model that can be used to quantify the effects of membrane morphology on separation efficiency, in terms of the performance (flux–throughput) curve of a membrane filter with known characteristics under given operating conditions, and by the graphs of particle concentration at pore outlet.
Our model accounts, in the simplest possible way, for variations in membrane-pore profiles. The pore profile variation in real membranes is undoubtedly highly complex: here we restrict attention to simple axisymmetric pore profiles characterized by depth-dependent initial radius
$a(x,0)$
, which span the entire membrane depth, and we investigate how filtration performance varies as these initial pore profiles change. Our results simulating filtration at constant pressure drop indicate that such variations in pore profile lead to different fouling patterns within the membrane. More importantly, if the initial pore radius at the top of membrane is small (pore profile
$a_{2}(x)$
in (3.2)), it can give rise to a marked decrease in filter performance, as quantified by the total amount of filtrate processed under the same operating conditions, as shown by figure 3(g). This figure, which summarizes results for five distinct (equal resistance) initial pore profiles, shows that the case where the initial pore profile is linear decreasing across the membrane, given by
$a(x,0)=a_{3}(x)$
(see (3.2)), gives significantly higher total throughput when compared with the other cases considered (initial pore profile uniform, linear increasing, concave/convex parabolic across the membrane).
Similar differences in performance, though less pronounced, are observed for the case where the total flux is prescribed (rather than the pressure drop), as shown in figure 9(f) in appendix B. Maintaining the same flux requires a significantly higher pressure drop for the linear increasing case
$a(x,0)=a_{2}(x)$
than for all other cases, while the case
$a(x,0)=a_{3}(x)$
, linear decreasing pore profile in (3.2), shows the best performance, requiring the lowest pressure drop to sustain the desired flux.
When studying the influence of the deposition (or particle–pore attraction) parameter
$\unicode[STIX]{x1D706}=8\unicode[STIX]{x1D6EC}\unicode[STIX]{x1D707}D^{2}/(P_{0}W)$
, care must be taken to track the concentration of small particles, both within the membrane and in particular at the downstream edge of the membrane (figures 5 and 6). A naive interpretation of figure 6(a) would suggest that smaller values of
$\unicode[STIX]{x1D706}$
are preferable, since these lead to greater total throughput. However, a glance at figure 5(b) confirms that if
$\unicode[STIX]{x1D706}$
is too small then a large concentration of small particles remain in the filtrate, which is likely undesirable. From the definition of
$\unicode[STIX]{x1D706}$
this could stem from several causes; for example, there may be insufficiently strong attractive forces between particles and membrane (
$\unicode[STIX]{x1D6EC}$
); the filter membrane may be too thin (
$D$
); or the pressure drop may be too high (
$P_{0}$
) so filtration is too fast to give good deposition.
This brings us to another major performance requirement of filtration: to achieve the desired level of particle separation from the feed solution. This separation level may vary from one application to another (sometimes a filtrate should be as clean as possible, with all impurities removed; at other times a threshold level of impurities may be tolerated, or it may be desirable to remove only a certain type of particles from the feed) and hence the best choice of filter may depend on the application. If we consider the simplest scenario in which the filtrate should be as clean as possible, while simultaneously maximizing throughput, then for a given tolerance level of impurities (maximum allowable concentration
$c(1,t)$
at pore outlet), the more throughput the filter gives, the better performance it has. Figure 5(a) shows that among all the initial pore profiles given in (3.2),
$a_{2}(x)$
and
$a_{3}(x)$
have the worst and the best performance respectively, under these criteria.
We note that the flux–throughput curves generated by our model are in good qualitative agreement with experimental data from the literature, as seen in, for example, Griffiths et al. (Reference Griffiths, Kumar and Stewart2014), Kumar, Martin & Kuriyel (Reference Kumar, Martin and Kuriyel2015) and many other works. This consistency between our model results and the experimental data gives us further confidence that our model, based as it is on first-principles assumptions about how fouling occurs, is sound, and provides a good basis for predictive simulations. While a more complicated model could perhaps provide more accurate predictions, our model has the advantage that it is simple and quick to simulate, offering a useful tool for investigating filter design characteristics.
Our predictive model leads naturally to questions of how membrane structure may be optimized. Defining optimal performance is very application-dependent, and in this paper we consider only a simple optimization, maximizing total throughput of filtrate over the filter lifetime as the pore profile is varied. We do not explicitly optimize also for particle removal efficiency, which is of course important; but we compare performance only of membranes with the same initial net resistance, assuming that particle removal is comparable for such membranes (as is borne out by the simulations shown in figures 5(a) and 7). We optimize only within a restricted class of pore profiles, presenting detailed results for linear pore profiles (preliminary work suggests that considering a larger class of profiles yields only marginal improvements). We note that our findings here are remarkably consistent with those of Griffiths et al. (Reference Griffiths, Kumar and Stewart2016), despite the many differences in approach. Collaborative work is ongoing to reach quantitative agreement between these different models.
There are of course many alternative approaches to optimization that could be considered, and a full investigation is beyond the scope of this paper. One could, as noted in § 3.2 above, extend the investigation to sweep through a larger parameter space in which the membrane thickness
$D$
and the size of the period-box
$2W$
are varied while keeping membrane resistance fixed. Another approach that could be interesting is to reverse time in the problem: if one assumes that the ‘optimal’ filtration scenario is that in which pore closure occurs uniformly along the length of the pore (such a scenario would maximize the time for which the pore is open), then one could run the model backwards to simulate the opening of an infinitesimally thin, parallel-sided pore. Stopping the simulation when the net resistance reaches a chosen value would then provide the optimal pore shape for that chosen resistance (of course, there are still other optimization questions here relating to varying the membrane thickness and the period-box size).
Finally, though our model represents an important first step in systematically accounting for internal membrane complexity, it must be emphasized that real membranes are much more complicated in structure than our simple assumptions allow, as is evident from a glance at figure 1. They may consist of many randomly oriented pores, which branch and reconnect, so that the feed solution takes a winding and tortuous path through the membrane rather than the simple flow assumed here. In future work, we will describe more sophisticated models, with branching and reconnecting pores, to better account for such internal membrane complexity with multiple membrane fouling modes (including so-called cake filtration, which occurs normally at a late stage of the filtration process). This model is able to capture blocking by large particles not only at the top surface of the membrane, but also at interior points of the membrane where pores branch or reconnect.
Acknowledgements
The authors acknowledge financial support from the NSF under grants DMS-1261596 and DMS-1615719, and from Pall Corporation. Several very helpful conversations with Drs I. Griffiths (University of Oxford), A. Kumar (Pall Corporation) and M. Hurwitz (IFFRISS, Cornell Engineering) are gratefully acknowledged, as are many helpful comments made by the reviewers. The work carried out in this paper arose from a problem presented at the 2014 Mathematical Problems in Industry workshop, held at NJIT under the support listed above.
Appendix A. Derivation of the advection model for small particles
In this appendix we justify the advection model for the concentration of small particles presented in (2.21), which in dimensionless form is (2.14). We begin from first principles, considering the flow of a suspension of small particles through a periodic array of identical channels (see figure 2). We seek solutions for the flow and particle concentration within a single channel of radius
$A(X,T)$
, in cylindrical coordinates
$(R,\unicode[STIX]{x1D703},X)$
in which properties vary only in the axial direction
$X$
, the radial direction
$R$
, and time
$T$
. Throughout this appendix, unstarred dependent variables are radially averaged over the pore cross-section, while the starred equivalents have radial dependence. Though the pore geometry changes due to the particle deposition, this occurs on a time scale much longer than that of the flow, hence we use a quasi-static model in which the flow domain
$0\leqslant R\leqslant A$
may be regarded as fixed. The pore velocity vector
$\boldsymbol{V}_{p}^{\ast }=(V_{p}^{\ast },0,U_{p}^{\ast })$
and the pressure
$P^{\ast }$
satisfy the Stokes equations, inertia being negligible in all scenarios of interest:
$$\begin{eqnarray}\unicode[STIX]{x1D735}P^{\ast }=\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{V}_{p}^{\ast },\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{V}_{p}^{\ast }=0,\quad 0\leqslant X\leqslant D,\quad 0\leqslant R\leqslant A,\end{eqnarray}$$
subject to
$$\begin{eqnarray}P^{\ast }|_{X=0}=P_{0},\quad P^{\ast }|_{X=D}=0,\quad \boldsymbol{V}_{p}^{\ast }=0\quad \text{at }R=A.\end{eqnarray}$$
The full advection–diffusion equation satisfied by the small-particle concentration
$C^{\ast }(R,X,T)$
is:
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}C^{\ast }}{\unicode[STIX]{x2202}T}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{Q}_{c}^{\ast },\quad \boldsymbol{Q}_{c}^{\ast }=-\unicode[STIX]{x1D6EF}\unicode[STIX]{x1D735}C^{\ast }+\boldsymbol{V}_{p}^{\ast }C^{\ast }+F^{\ast }C^{\ast }\boldsymbol{e}_{\boldsymbol{ r}},\end{eqnarray}$$
where
$\boldsymbol{Q}_{c}^{\ast }$
is the total particle flux,
$\unicode[STIX]{x1D6EF}$
is the diffusion coefficient of the small particles in the feed solution and
$F^{\ast }$
is the radial particle drift speed induced by interaction with the pore wall. This term is intended to describe any forces, such as electrostatic interactions and perhaps van der Waals’ forces, that act to attract particles towards the wall where they can adhere. It could be argued that this attraction force should act in the direction locally perpendicular to the wall rather than radially, requiring a more careful analysis using the normal vector; but since the pore is slender the normal direction is very close to the radial direction and, to the order we consider in our asymptotics, the end result may be shown to be the same with the purely radial term in (A 3). The boundary conditions are
$$\begin{eqnarray}C^{\ast }(R,0,T)=C_{0},\quad \left.\frac{\unicode[STIX]{x2202}C^{\ast }}{\unicode[STIX]{x2202}R}\right|_{R=0}=0,\quad \boldsymbol{Q}_{c}^{\ast }\boldsymbol{\cdot }\boldsymbol{n}=\frac{\unicode[STIX]{x1D6EC}}{2}C^{\ast }\quad \text{at }R=A(\text{wall deposition}),\end{eqnarray}$$
for some constant
$\unicode[STIX]{x1D6EC}$
. Consistent with our quasi-static assumption we solve the steady-state version of (A 3),
$$\begin{eqnarray}\unicode[STIX]{x1D6EF}\left(\frac{1}{R}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}R}\left(R\frac{\unicode[STIX]{x2202}C^{\ast }}{\unicode[STIX]{x2202}R}\right)+\frac{\unicode[STIX]{x2202}^{2}C^{\ast }}{\unicode[STIX]{x2202}X^{2}}\right)=V_{p}^{\ast }\frac{\unicode[STIX]{x2202}C^{\ast }}{\unicode[STIX]{x2202}R}+U_{p}^{\ast }\frac{\unicode[STIX]{x2202}C^{\ast }}{\unicode[STIX]{x2202}X}+\frac{1}{R}\frac{\unicode[STIX]{x2202}(RF^{\ast }C^{\ast })}{\unicode[STIX]{x2202}R},\end{eqnarray}$$
where we have used the continuity equation,
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{V}_{p}^{\ast }=0$
. Furthermore the wall deposition boundary condition in (A 4), by use of the zero-velocity boundary condition in (A 2) and calculation of
$\boldsymbol{n}\Vert \unicode[STIX]{x1D735}(R-A)$
, becomes
$$\begin{eqnarray}-\unicode[STIX]{x1D6EF}\frac{\unicode[STIX]{x2202}C^{\ast }}{\unicode[STIX]{x2202}R}+\unicode[STIX]{x1D6EF}\frac{\unicode[STIX]{x2202}C^{\ast }}{\unicode[STIX]{x2202}X}\frac{\unicode[STIX]{x2202}A}{\unicode[STIX]{x2202}X}+F^{\ast }C^{\ast }=\frac{\unicode[STIX]{x1D6EC}}{2}C^{\ast }\quad \text{at }R=A\text{ (wall deposition)}.\end{eqnarray}$$
We non-dimensionalize using the scalings
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \boldsymbol{V}_{p}^{\ast }=\frac{\unicode[STIX]{x03C0}P_{0}W^{2}}{32\unicode[STIX]{x1D707}D}(\unicode[STIX]{x1D716}v_{p}^{\ast },0,u_{p}^{\ast }),\quad X=Dx,\quad (A,R)=W(a,r),\\ \displaystyle P^{\ast }=P_{0}p^{\ast },\quad C^{\ast }=C_{0}c^{\ast },\quad F^{\ast }=\frac{\unicode[STIX]{x03C0}P_{0}W^{2}}{32\unicode[STIX]{x1D707}D}\unicode[STIX]{x1D716}f^{\ast },\end{array}\!\right\}\end{eqnarray}$$
where
$\unicode[STIX]{x1D716}=W/D\ll 1$
. This gives the dimensionless form of (A 1), (A 2), (A 4), (A 5) and (A 6) as
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{32}{\unicode[STIX]{x03C0}\unicode[STIX]{x1D716}^{2}}\frac{\unicode[STIX]{x2202}p^{\ast }}{\unicode[STIX]{x2202}r}=\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(r\frac{\unicode[STIX]{x2202}v_{p}^{\ast }}{\unicode[STIX]{x2202}r}\right)-\frac{v_{p}^{\ast }}{r^{2}}+\unicode[STIX]{x1D716}^{2}\frac{\unicode[STIX]{x2202}^{2}v_{p}^{\ast }}{\unicode[STIX]{x2202}x^{2}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{32}{\unicode[STIX]{x03C0}}\frac{\unicode[STIX]{x2202}p^{\ast }}{\unicode[STIX]{x2202}x}=\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(r\frac{\unicode[STIX]{x2202}u_{p}^{\ast }}{\unicode[STIX]{x2202}r}\right)+\unicode[STIX]{x1D716}^{2}\frac{\unicode[STIX]{x2202}^{2}u_{p}^{\ast }}{\unicode[STIX]{x2202}x^{2}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{1}{r}\frac{\unicode[STIX]{x2202}(rv_{p}^{\ast })}{\unicode[STIX]{x2202}r}+\frac{\unicode[STIX]{x2202}u_{p}^{\ast }}{\unicode[STIX]{x2202}x}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}p^{\ast }|_{x=0}=1,\quad p^{\ast }|_{x=1}=0,\quad u_{p}^{\ast }=v_{p}^{\ast }=0\quad \text{at }r=a(x),\end{eqnarray}$$
$$\begin{eqnarray}Pe\left(\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(r\frac{\unicode[STIX]{x2202}c^{\ast }}{\unicode[STIX]{x2202}r}\right)+\unicode[STIX]{x1D716}^{2}\frac{\unicode[STIX]{x2202}^{2}c^{\ast }}{\unicode[STIX]{x2202}x^{2}}\right)=v_{p}^{\ast }\frac{\unicode[STIX]{x2202}c^{\ast }}{\unicode[STIX]{x2202}r}+u_{p}^{\ast }\frac{\unicode[STIX]{x2202}c^{\ast }}{\unicode[STIX]{x2202}x}+\frac{1}{r}\frac{\unicode[STIX]{x2202}(rf^{\ast }c^{\ast })}{\unicode[STIX]{x2202}r},\quad \text{where }Pe=\frac{32\unicode[STIX]{x1D6EF}\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x03C0}P_{0}W^{2}},\end{eqnarray}$$
$$\begin{eqnarray}c^{\ast }(r,0,t)=1,\quad \left.\frac{\unicode[STIX]{x2202}c^{\ast }}{\unicode[STIX]{x2202}r}\right|_{r=0}=0,\quad -Pe\frac{\unicode[STIX]{x2202}c^{\ast }}{\unicode[STIX]{x2202}r}+\unicode[STIX]{x1D716}^{2}Pe\frac{\unicode[STIX]{x2202}c^{\ast }}{\unicode[STIX]{x2202}x}\frac{\unicode[STIX]{x2202}a}{\unicode[STIX]{x2202}x}+f^{\ast }c^{\ast }=\unicode[STIX]{x1D706}_{1}c^{\ast },\quad \text{at }r=a,\end{eqnarray}$$
where
$\unicode[STIX]{x1D706}_{1}=16\unicode[STIX]{x1D6EC}\unicode[STIX]{x1D707}D^{2}/\unicode[STIX]{x03C0}P_{0}W^{3}$
. We exploit asymptotic analysis of this model, based on the assumption
$\unicode[STIX]{x1D716}=W/D\ll 1$
and in the distinguished limit
$\hat{Pe}=\unicode[STIX]{x1D716}Pe=O(1)$
, representing a specific balance between advective and diffusive particle transport. We expand the dependent variables in powers of
$\unicode[STIX]{x1D716}$
, e.g.
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}c^{\ast }(r,x,t)=c_{0}^{\ast }(r,x,t)+\unicode[STIX]{x1D716}c_{1}^{\ast }(r,x,t)+\unicode[STIX]{x1D716}^{2}c_{2}^{\ast }(r,x,t)+\cdots \,,\\ f^{\ast }(r)=f_{0}^{\ast }(r)+\unicode[STIX]{x1D716}f_{1}^{\ast }(r)+\unicode[STIX]{x1D716}^{2}f_{2}^{\ast }(r)+\cdots \,,\end{array}\right\} & & \displaystyle\end{eqnarray}$$
etc. Solving (A 8)–(A 11) gives the pore velocity and pressure at leading order (
$u_{p0}^{\ast }$
,
$v_{p0}^{\ast }$
and
$p_{0}$
, respectively) as
$$\begin{eqnarray}\displaystyle u_{p0}^{\ast }=\frac{8}{\unicode[STIX]{x03C0}}\frac{\unicode[STIX]{x2202}p_{0}^{\ast }}{\unicode[STIX]{x2202}x}(r^{2}-a^{2}),\quad v_{p0}^{\ast }=\frac{2}{\unicode[STIX]{x03C0}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left(\frac{\unicode[STIX]{x2202}p_{0}^{\ast }}{\unicode[STIX]{x2202}x}(r^{3}-2a^{2}r)\right),\quad p_{0}^{\ast }=\frac{\displaystyle \int _{x}^{1}\frac{\text{d}x^{\prime }}{a^{4}}}{\displaystyle \int _{0}^{1}\frac{\text{d}x^{\prime }}{a^{4}}}.\quad & & \displaystyle\end{eqnarray}$$
At leading-order equation (A 12) reduces to
$c_{0}^{\ast }(r,x,t)=c(x,t)$
, consistent with leading order in the boundary conditions (A 13). At order one, we obtain
$$\begin{eqnarray}\frac{\hat{Pe}}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(r\frac{\unicode[STIX]{x2202}c_{1}^{\ast }}{\unicode[STIX]{x2202}r}\right)=u_{p0}^{\ast }\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}x}+\frac{c}{r}\frac{\unicode[STIX]{x2202}(rf_{0}^{\ast })}{\unicode[STIX]{x2202}r},\end{eqnarray}$$
where
$u_{p0}^{\ast }$
is the leading-order dimensionless axial pore velocity. Using the operator
$(1/\unicode[STIX]{x03C0}a^{2})\int _{0}^{2\unicode[STIX]{x03C0}}\int _{0}^{a}\,\cdot \,r\,\text{d}r\,\text{d}\unicode[STIX]{x1D703}$
to take the cross-sectional average of (A 16), we obtain
$$\begin{eqnarray}\left.\frac{2\hat{Pe}}{a}\frac{\unicode[STIX]{x2202}c_{1}^{\ast }}{\unicode[STIX]{x2202}r}\right|_{r=a}=u_{p}\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}x}+\frac{2c}{a}f_{0}^{\ast }(a),\end{eqnarray}$$
where
$u_{p}$
denotes the cross-sectional average of
$u_{p0}^{\ast }$
, obtained from (A 15)
$$\begin{eqnarray}u_{p}=-\frac{4}{\unicode[STIX]{x03C0}}a^{2}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x},\end{eqnarray}$$
where, in line with our notation,
$p$
is the cross-sectional average of
$p_{0}^{\ast }$
(note that here
$p\equiv p_{0}^{\ast }$
). In addition the second boundary condition in (A 13) gives, at
$O(1)$
$$\begin{eqnarray}-\hat{Pe}\frac{\unicode[STIX]{x2202}c_{1}^{\ast }}{\unicode[STIX]{x2202}r}+cf_{0}^{\ast }=\unicode[STIX]{x1D706}_{1}c,\quad \text{at }r=a.\end{eqnarray}$$
Combining (A 17) and (A 19) gives the desired result (2.21)
$$\begin{eqnarray}u_{p}\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}x}=-\hat{\unicode[STIX]{x1D706}}\frac{c}{a},\quad \hat{\unicode[STIX]{x1D706}}=\frac{32\unicode[STIX]{x1D6EC}\unicode[STIX]{x1D707}D^{2}}{\unicode[STIX]{x03C0}P_{0}W^{3}},\quad c_{0}(0,t)=1.\end{eqnarray}$$
Appendix B. Constant flux case
Here we briefly outline how the results change if conditions of constant flux, rather than constant pressure drop, are applied. We first outline the modifications to the model, then present a sample simulation.
Figure 9. Simulations at constant flux: the pore radius at several different times at different final blocking times (
$t_{f}$
, indicated in the legends) for different initial pore radius profiles: (a)
$a_{1}(x,0)=0.904$
, (b)
$a_{2}(x,0)=0.16x+0.83$
, (c)
$a_{3}(x,0)=0.99-0.16x$
, (d)
$a_{4}(x,0)=0.874+0.39(x-0.5)^{2}$
, (e)
$a_{5}(x,0)=0.933-0.33(x-0.5)^{2}$
, and (f) inverse pressure drop versus throughput for those initial pore radius profiles with homogeneous distributions of large-particle sizes, (2.30) (
$g(s)=0$
),
$\tilde{\unicode[STIX]{x1D706}}=2$
,
$\unicode[STIX]{x1D6FD}=0.1$
and
$\unicode[STIX]{x1D70C}_{b}=2$
.
B.1 The model for specified flux
The original model (2.2)–(2.17) still holds, but now
$P_{0}$
in (2.2) must be considered a function of time,
$P_{0}(T)$
, while (2.3) integrates directly to give
$$\begin{eqnarray}\unicode[STIX]{x03C0}A^{2}U_{p}=Q_{pore},\end{eqnarray}$$
where
$Q_{pore}$
is the constant flux per pore. We non-dimensionalize the model using the same scalings as in (2.18), except for
$$\begin{eqnarray}P=\frac{8\unicode[STIX]{x1D707}DQ_{pore}}{\unicode[STIX]{x03C0}W^{4}}p,\quad (U,U_{p})=\frac{Q_{pore}}{4W^{2}}(1,u_{p}).\end{eqnarray}$$
The resulting model is easily reduced to:
$$\begin{eqnarray}\frac{\text{d}n}{\text{d}t}=-\!\left.n\left(\frac{1-n}{1+\unicode[STIX]{x1D70C}_{b}a^{4}}+n\right)^{-1}(1-g(a))\right|_{x=0},\quad \frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}x}=-\tilde{\unicode[STIX]{x1D706}}ca,\quad \frac{\unicode[STIX]{x2202}a}{\unicode[STIX]{x2202}t}=-\tilde{\unicode[STIX]{x1D6FD}}c,\end{eqnarray}$$
where
$$\begin{eqnarray}\tilde{\unicode[STIX]{x1D6FD}}=\frac{1}{Q_{pore}G_{\infty }},\quad \tilde{\unicode[STIX]{x1D706}}=\frac{\unicode[STIX]{x03C0}\unicode[STIX]{x1D6EC}WD}{Q_{pore}},\end{eqnarray}$$
with modified Darcy pressure
$p$
within the membrane given by
$$\begin{eqnarray}p=\int _{x}^{1}\frac{\text{d}x^{\prime }}{\displaystyle a^{4}\left(\frac{1-n}{1+\unicode[STIX]{x1D70C}_{b}a^{4}}+n\right)}.\end{eqnarray}$$
The model (B 3) must be solved subject to the same boundary and initial conditions (2.29). Note that (B 5) allows the pressure
$p(0,t)$
at the membrane inlet (the dimensionless pressure drop in this constant flux case) to be evaluated.
B.2 Results
Figure 9(a–e) shows results for the same initial pore profiles given in (3.2). Figure 9(f) shows the inverse pressure drop versus throughput for each of those pore profiles.
The results differ quite significantly from those for the constant pressure case. In contrast to those simulations, the pore radius evolution is now much more uniform along the pore length. Pore closure still always occurs first at the upstream end of the pores for the cases shown here, but with a smaller value of
$\tilde{\unicode[STIX]{x1D706}}$
this is not inevitable. Since the total flux through the system is held constant, the flux–throughput graph gives no characterization of the system in these simulations, hence we instead plot how the pressure drop rises over time as blocking occurs in order to maintain the specified flow rate (figure 9
e). Mathematically, the pressure must go to infinity within finite time to sustain the same flux as total blockage is reached, but this is of course not practical. In reality the driving pressure is increased until some specified maximum value (based on practical constraints of the system under consideration) is reached. Following this, the system then reverts to the fixed pressure drop operation, with the pressure fixed at this maximal value, and the subsequent behaviour can then be modelled as discussed in § 3.1. As with the constant pressure simulations, the best overall performance (in terms of efficiency) is provided by the pores of monotone decreasing radius (profile
$a_{3}(x)$
in (3.2)), and the worst performance by pores of monotone increasing radius (profile
$a_{2}(x)$
). However, the differences in performance are less pronounced than for the constant pressure case.
