1 Introduction
Lateral migration and focusing of neutrally buoyant solid particles at moderate Reynolds number (
$Re$
) in a confined pressure-driven flow is a well-known phenomenon first documented by Segré and Silberberg in 1961. Specifically, hydrodynamic inertial stresses cause particles to laterally migrate across streamlines and ultimately focus into distinct locations in the channel (Ho & Leal Reference Ho and Leal1974). In addition to the seminal work on inertial migration and focusing (Saffman Reference Saffman1965; Cox & Brenner Reference Cox and Brenner1968; Ho & Leal Reference Ho and Leal1974), there have also been a few recent reviews highlighting progress (Amini, Lee & Di Carlo Reference Amini, Lee and Di Carlo2014; Martel & Toner Reference Martel and Toner2014; Zhang et al.
Reference Zhang, Yan, Yuan, Alici, Nguyen, Warkiani and Li2016). However, comparatively, there have only been a few studies on the motion of inertial particles in the presence of a permeate flow. Systems with such flows are ubiquitous in applications related to pressure-driven membrane filtration for the separation of particles and cells from liquid suspensions as well as in a multitude of other areas, including wastewater treatment (Chang & Kim Reference Chang and Kim2005), food (Fernández García, Álvarez Blanco & Riera Rodríguez Reference Fernández García, Álvarez Blanco and Riera Rodríguez2013) and beverage (Pall Corporation 2015) processing and biotechnology (Charcosset Reference Charcosset2006; Palmer, Sun & Harris Reference Palmer, Sun and Harris2009).
In general, to precisely solve for the forces on a particle in inertial migration problems, one relies on simulation of the combined fluid–particle interaction problem. For example, Chun & Ladd (Reference Chun and Ladd2006) used the lattice-Boltzmann method for
$Re$
ranging from 100 to 1000 to show that spherical particles migrate to one of a finite number of equilibrium locations in the cross-section of a non-porous square channel, where only the face-centred equilibrium locations are stable for lower
$Re$
. Di Carlo et al. (Reference Di Carlo, Edd, Humphry, Stone and Toner2009) performed simulations to show how particle equilibrium locations vary with particle sizes (a trend not captured by previous asymptotic approaches). Similarly, Liu et al. (Reference Liu, Hu, Jiang and Sun2015) performed exhaustive simulations to determine the non-trivial relationship between the focusing location of particles and their size, channel aspect ratio, and
$Re$
. These numerical studies – fitted into generalised formulae – shed light on how the forces on a migrating particle are distributed within a channel and how this distribution depends on particle size. Such findings are crucial for practical device design.
Figure 1. (a) Schematic illustration of our system. We model a square channel with average velocity
$U$
. In the channel a spherical particle of diameter
$a$
is migrating within the confines of the bounding walls where two of the walls are permeable (
$y{-}z$
plane) and allow flow to penetrate at a constant rate of
$U_{W}$
. For each
$x$
–
$z$
location of the particle in the channel cross-section, the lateral lift forces (
$F_{z}$
and
$F_{x}$
) are calculated. (b) Mesh sensitivity analysis showing that the calculated lift forces have converged and are, thus, insensitive to the DOF in the model. The inset depicts the error relative to a model with
$1.6\times 10^{6}$
DOF. (
$Re=100$
,
$a/W=0.1$
.)
With the addition of permeate flow, the interplay between the effects of axial and permeate flow (see figure 1 a) results in a much richer particle behaviour than the case of pure inertial migration of a particle in axial flow (Altena & Belfort Reference Altena and Belfort1984; Drew, Schonberg & Belfort Reference Drew, Schonberg and Belfort1991; Lebedeva & Asmolov Reference Lebedeva and Asmolov2011). Earlier theoretical work by Belfort and co-workers accounts for the effect of wall porosity on particle motion in a two-dimensional geometry (Altena & Belfort Reference Altena and Belfort1984; Drew et al. Reference Drew, Schonberg and Belfort1991) and has been validated experimentally (Otis et al. Reference Otis, Altena, Mahar and Belfort1986). In these theoretical studies, researchers employed asymptotic analysis to derive expressions for the forces on a migrating particle as a sum of inertially induced forces and permeate drag. While useful from a theoretical perspective, the researchers were unable to precisely predict the lateral lift forces of highly confined particles in a porous channel. This limitation is particularly critical since precise force values across the cross-section are required for efficient device design that exploits the lateral migration of particles to separate and concentrate particles with high specificity. More recently, Garcia & Pennathur (Reference Garcia and Pennathur2017) experimentally demonstrated that a permeate flow can drastically alter the migration and focusing of confined particles in inertial flow, with particle size significantly impacting the resulting migration behaviour. However, no model to date has accurately described this behaviour. The ability to construct a simple model that reliably predicts these forces across a wide range of conditions serves as the overarching motivation for the present work.
Based on observations from a finite number of full-scale simulations, we develop a linear model to efficiently predict the lateral forces acting on a neutrally buoyant spherical particle migrating in a pressure-driven porous channel. We first employ a previously utilised full-physics simulation approach (Di Carlo et al.
Reference Di Carlo, Edd, Humphry, Stone and Toner2009; Gossett et al.
Reference Gossett, Tse, Dudani, Goda, Woods, Graves and Di Carlo2012; Hood, Lee & Roper Reference Hood, Lee and Roper2015; Kim et al.
Reference Kim, Lee, Wu, Nam, Di Carlo and Lee2016) to numerically simulate a spherical particle translating in a microchannel of square cross-section. We additionally validate this model against previously published work and experiments, then exploit the model to study the lateral forces in a system with permeate flow of varying magnitude and direction. We observe that the relative permeate flow,
$\unicode[STIX]{x1D6FE}$
, can generalise the data over many particle sizes and
$Re$
, which provides a conceptual basis for developing a new linear model. We finally validate our linear model using the full-physics simulation, and determine that the limits of applicability are well within those of most experimental systems. To prove as much, we design a microfluidic system to experimentally validate the results of the linear model.
2 Full-physics simulations of particle–fluid interaction
The present work focuses on the combined inertial and permeate migration of a single neutrally buoyant spherical particle in Poiseuille flow within the confining geometry of a porous channel of square cross-section (
$W\times W$
) and length
$L_{C}=10W$
. Here the particle of diameter
$a$
is translating with velocity
$\boldsymbol{U}_{P}=[0,U_{P},0]$
and angular velocity
$\unicode[STIX]{x1D734}$
in a flow of average (local) velocity
$U=(U_{in}+U_{out})/2$
, where
$U_{in}$
and
$U_{out}$
are the average axial flow velocities at the inlet and outlet of the porous domain respectively. The porous channel allows for a lateral flow to exist that emanates from only two parallel walls at a constant velocity
$U_{W}$
. Here we define the channel Reynolds number
$Re=\unicode[STIX]{x1D70C}UW/\unicode[STIX]{x1D707}$
and the relative permeate flow
$\unicode[STIX]{x1D6FE}=U_{W}/U$
, where
$\unicode[STIX]{x1D70C}$
and
$\unicode[STIX]{x1D707}$
are the fluid density and viscosity respectively.
We consider a moving frame of reference translating with the accelerating particle. (In this moving reference frame, the lateral flow extends a length
$3W$
as depicted in figure 1.) In this reference frame, we solve the steady Navier–Stokes equations:
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}+\boldsymbol{U}_{P}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\bar{\boldsymbol{u}})=\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}-\unicode[STIX]{x1D735}p, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$
where
$p$
is the fluid pressure,
$\boldsymbol{u}$
is the fluid velocity in the reference frame of the translating particle and
$\bar{\boldsymbol{u}}$
is the undisturbed flow. Note that the second inertial term arises due to the acceleration of our chosen reference frame as shown in more detail in appendix A.
The velocity of the suspended particle (
$U_{P}$
and
$\unicode[STIX]{x1D734}$
) can be self-consistently determined by setting conditions such that the axial motion satisfies a drag constraint
$F_{y}=0$
(2.3) and its rotational motion satisfies a torque constraint
$\unicode[STIX]{x1D70F}_{x}=\unicode[STIX]{x1D70F}_{y}=\unicode[STIX]{x1D70F}_{z}=0$
(2.4):
$$\begin{eqnarray}\displaystyle & \displaystyle F_{y}\equiv \boldsymbol{e}_{y}\boldsymbol{\cdot }\int _{s}\boldsymbol{n}_{\boldsymbol{r}}\boldsymbol{\cdot }\unicode[STIX]{x1D64F}\,\text{d}s=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70F}_{i}\equiv \boldsymbol{e}_{i}\boldsymbol{\cdot }\int _{s}(\boldsymbol{r}-\boldsymbol{r}_{p})\times \boldsymbol{n}_{\boldsymbol{r}}\boldsymbol{\cdot }\unicode[STIX]{x1D64F}\,\text{d}s=0\quad \text{where}~i=x,y,z, & \displaystyle\end{eqnarray}$$
where the integrals are over the surface of the sphere,
$\boldsymbol{r}_{p}$
is the particle position vector that points from the centre of the channel to the centre of the particle and
$\boldsymbol{n}_{\boldsymbol{r}}=(\boldsymbol{r}-\boldsymbol{r}_{p})/|\boldsymbol{r}-\boldsymbol{r}_{p}|$
is the unit normal at each point on the surface of the sphere. Here
$\unicode[STIX]{x1D64F}$
is the total stress tensor and, for an incompressible Newtonian fluid, is given by
$\unicode[STIX]{x1D64F}=-p\unicode[STIX]{x1D644}+2\unicode[STIX]{x1D707}\unicode[STIX]{x1D640}$
, where
$\unicode[STIX]{x1D640}=(1/2)(\unicode[STIX]{x1D735}\boldsymbol{u}+\unicode[STIX]{x1D735}\boldsymbol{u}^{\text{T}})$
is the rate of strain tensor.
The boundary conditions of this problem are in the reference frame of the translating particle. Therefore, the no slip condition on the wall is
$$\begin{eqnarray}\boldsymbol{u}=-U_{p}\boldsymbol{e}_{y}\quad \text{on all walls}.\end{eqnarray}$$
Additionally, the channel has two parallel porous walls that can advect fluid into or out of the channel at a constant rate of
$U_{W}$
. Thus the condition on these walls also must satisfy
$$\begin{eqnarray}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}=U_{W},\end{eqnarray}$$
where
$\boldsymbol{n}$
is the wall unit normal pointing out of the channel and a positive value of
$U_{W}$
indicates flow out of the channel. The no-slip condition on the rotating particle is enforced by assigning a velocity to the surface of the sphere corresponding to that of rigid body rotation at angular velocity
$\unicode[STIX]{x1D734}$
:
$$\begin{eqnarray}\boldsymbol{u}_{surface}=\unicode[STIX]{x1D734}\times (\boldsymbol{r}-\boldsymbol{r}_{p}).\end{eqnarray}$$
Finally, far from the particle the flow field satisfies
$$\begin{eqnarray}\boldsymbol{u}=\bar{\boldsymbol{u}}-U_{p}\boldsymbol{e}_{y}.\end{eqnarray}$$
To solve for the unknowns (i.e.
$\boldsymbol{u}$
,
$p$
,
$U_{p}$
and
$\unicode[STIX]{x1D734}$
) we couple equations (2.3) and (2.4) to the fluid flow equations and solve directly using COMSOL Multiphysics software. This procedure is performed for a lattice of discrete positions of the particle within the cross-section of the channel (using only a quarter of the domain via symmetry arguments to minimise computational effort). To calculate the lift force on the particle, we integrate the surface stresses on the particle in the appropriate direction (
$x$
or
$z$
):
$$\begin{eqnarray}F_{i}=\boldsymbol{e}_{i}\boldsymbol{\cdot }\int _{s}\boldsymbol{n}_{r}\boldsymbol{\cdot }\unicode[STIX]{x1D64F}\,\text{d}s\quad \text{where}~i=x,z.\end{eqnarray}$$
3 Numerical results
3.1 Mesh convergence analysis
We discretise the fluid domain using tetrahedral elements in COMSOL Multiphysics. The surface of the sphere is discretised into triangular boundary elements with six boundary layer elements at the surface of the sphere to accommodate large gradients. We use quadratic basis functions for representing the velocity solution and linear basis functions for the pressure. The resulting discretisation has approximately
$3.6\times 10^{5}$
degrees of freedom (DOF). To show that the calculated results were independent of mesh density we increased the number of DOF in our model by varying the boundary elements, boundary layer and fluid domain tetrahedral density. The resulting calculated lift forces are shown in figure 1(b), and show no significant differences (
${\leqslant}1\,\%$
change) in the spatial lift force amongst the different cases.
Figure 2. (a) Inertial forces calculated in a single quadrant of the square cross-section for a particle of
$a/W=0.10$
at
$Re=100$
. Here the black small squares indicate the locations of the stable equilibrium points, i.e. locations where the lift forces go to zero. (b) A comparison of the
$x$
-component inertial forces along the
$z/W=0$
axis between the present study and that of Di Carlo et al. (Reference Di Carlo, Edd, Humphry, Stone and Toner2009), showing good agreement at
$Re=38$
for two particles
$a/W=0.30$
and
$a/W=0.38$
. (c) Experimental measurements (discrete points) from Garcia & Pennathur (Reference Garcia and Pennathur2017) of the inertial focusing locations (
$x_{eq}/W$
) as a function of
$Re$
for three particle sizes (
$a/W$
) with overlaid corresponding numerical simulations (solid lines).
3.2 Model validation
To demonstrate the validity of the model, we performed simulations of a particle within the channel with no permeate flow (
$U_{W}=0$
), since this reduces to the well-studied case of an impermeable channel (Di Carlo et al.
Reference Di Carlo, Edd, Humphry, Stone and Toner2009; Miura, Itano & Sugihara-Seki Reference Miura, Itano and Sugihara-Seki2014; Shichi et al.
Reference Shichi, Yamashita, Seki, Itano and Sugihara-Seki2017). Typically in experimentally relevant flows for inertial microfluidic systems, particle diameters are a significant fraction of the channel width (
$a/W\geqslant 0.10$
) and generally
$Re$
ranges from
${\sim}10$
to 100. In the case of an impermeable channel under these conditions, particles focus to four symmetric equilibrium positions near the centre of each wall face and approximately
$0.3W$
away from the channel centre (Di Carlo et al.
Reference Di Carlo, Edd, Humphry, Stone and Toner2009). Figure 2(a) provides a detailed map of the spatially varying inertial lift forces in a single quadrant of a channel. With this map, we can identify the location of zero lateral lift force, which agrees with previous experimental studies (Miura et al.
Reference Miura, Itano and Sugihara-Seki2014; Shichi et al.
Reference Shichi, Yamashita, Seki, Itano and Sugihara-Seki2017). In figure 2(b), we compare the spatial variation of the lift forces on the particle which show good agreement with data from Di Carlo et al. (Reference Di Carlo, Edd, Humphry, Stone and Toner2009). Finally, we show that our model can correctly predict changes in the locations of equilibrium with increasing
$Re$
(Asmolov Reference Asmolov1999) and increasing particle size (Di Carlo et al.
Reference Di Carlo, Edd, Humphry, Stone and Toner2009) in figure 2(c), where we compare our previous experimental data (Garcia & Pennathur Reference Garcia and Pennathur2017) with current simulation results.
3.3 Permeate flow results
We next consider the effects of permeate flow, for a range of
$Re$
,
$\unicode[STIX]{x1D6FE}$
(
$\unicode[STIX]{x1D6FE}=U_{W}/U$
) and particle sizes. We consider
$Re=25$
, 50, 75, 100,
$\unicode[STIX]{x1D6FE}=\pm 0.003,\pm 0.002,\pm 0.001,0$
and particle sizes
$a/W=0.05,0.1,0.15$
. Figure 3(a,b) plots the lateral lift force vectors for a subset of the simulations. Figure 3(c–h) shows the equilibrium locations for the full range of
$\unicode[STIX]{x1D6FE}$
and
$Re$
, with figure 3(a,b) showing full force vector fields for the extreme cases. At large
$Re$
and large
$a/W$
, the vector fields are only slightly disturbed by the presence of the permeate flow when compared to the case with no permeate flow (
$\unicode[STIX]{x1D6FE}=0$
) (e.g. figure 3
a), presumably because inertial lift forces are more dominant (since the lift force
$F_{L}\sim a^{4}$
(Asmolov Reference Asmolov1999)). Conversely, for flows with less inertia (i.e. smaller
$Re$
and
$a/W$
) (e.g. figure 3
b), the vector fields are highly disturbed due to the substantial effects of permeate flow. Specifically, for
$\unicode[STIX]{x1D6FE}=0.003$
(e.g. figure 3
b right), we observe a complete suction where the force equilibrium coincides with the location where wall contact occurs and is when all ‘streamlines’ appear to be guided towards the porous wall (
$x/W=0.5$
). (We do not actually model the physics of wall contact. When referring to wall contact we do so in the context of the particle equilibrium location, i.e.
$x_{eq}/W\geqslant 0.5-a/2W$
.) Conversely, when
$\unicode[STIX]{x1D6FE}=-0.003$
(e.g. figure 3
b left), there is a reversal of ‘streamlines’, now directed towards the centreline (
$x/W=0$
).
Figure 3. Lateral force vector fields spanning the extremities of the parameter space
$\unicode[STIX]{x1D6FE}=[-0.003:0.001:0.003]$
and
$Re=[25:25:100]$
for two particle sizes: (a)
$a/W=0.15$
(b) and
$a/W=0.05$
. The blue ‘streamlines’ are for visualisation of the vector fields which are bounded by a
$0.45W$
$\times$
$0.45W$
box. (c–e) Locations of
$x$
-equilibrium for all direct simulations in this study along the
$z/W=0$
axis for (c)
$a/W=0.15$
, (d)
$a/W=0.10$
and (e)
$a/W=0.05$
. Note that the equilibrium shift is in the direction of the permeate flow. (f–h) Locations of
$z$
-equilibrium for all direct simulations in this study along the
$x/W=0$
axis for (f)
$a/W=0.15$
, (g)
$a/W=0.10$
and (h)
$a/W=0.05$
.
Figure 3(c–h) captures all these observations as we plot the
$x$
- and
$z$
-force equilibrium location (
$x_{eq}/W$
and
$z_{eq}/W$
) for all
$Re$
and
$a/W$
. Through this representation we can clearly see that as permeate flow increases, the equilibrium location of each vector field shifts in the direction of the permeate flow (figure 3
c–e). Furthermore, in agreement with the qualitative observation in figure 3(a,b), the shift in equilibrium location becomes less sensitive with increasing
$a/W$
and
$Re$
. We also see a slight shift in the
$z$
-equilibrium with the addition of a permeate flow, but in general the
$z$
-equilibrium shift is less sensitive than the
$x$
-equilibrium counterpart. One thing to note is that the
$z$
-equilibrium shifts either towards or away from the wall as a function of
$\unicode[STIX]{x1D6FE}$
, depending on the flow
$Re$
. In general, there is a clear dependence on two factors, the first being the magnitude of permeate forces (set by
$U_{W}$
) and the second being the magnitude of inertial forces (set by
$U$
). The combination of these two factors dictates the behaviour of the particles in our geometry, ultimately leading to non-trivial force fields.
Figure 4. (a) Residual force plot. Each curve represents a different value of
$\unicode[STIX]{x1D6FE}$
, where the residual is the difference in the force between a particle in the presence of permeate flow and a particle in a channel with no permeate flow (
$\unicode[STIX]{x1D6FE}=0$
) at different
$x$
locations. (b) The force residuals for all
$\unicode[STIX]{x1D6FE}$
normalised by a characteristic drag force resulting in three distinct master curves corresponding to each particle size at
$Re=100$
. (c) Normalised force residuals for four distinct
$Re$
at a constant particle size (
$a/W=0.10$
). Note that normalising in this manner results in the collapse of the curves into a single ‘master curve’. Due to noise in our simulations, we averaged values in both (b) and (c) as shown in the coloured lines. Black curves underneath represent the raw normalisations.
3.4 Linearised model
Solving for the effect of combined inertial and permeate forces with full-physics numerical simulations can quickly become computationally infeasible. For instance, the results in the previous sub-section involved performing
$4(Re)\times 3(a/W)\times 9(x/W)\times 9(z/W)\times 7(\unicode[STIX]{x1D6FE})=6804$
simulations, which required expending substantial computational resources. The computational requirements become even larger when modelling dynamic processes where a particle may be migrating in a continuously varying flow field, or in which
$\unicode[STIX]{x1D6FE}$
is spatially and/or temporally changing. Therefore, we develop a linearised model (based on observations from the previous sub-section) to produce quantitatively similar results to the full simulation, but with significantly less computational requirements.
We consider a linearisation of the lateral lift force
$\boldsymbol{F}(\unicode[STIX]{x1D6FE},Re,a/W,x/W,z/W)$
about
$\unicode[STIX]{x1D6FE}=0$
, i.e. impermeable channel case. There are multiple reasons for choosing such a linearisation strategy: (i) this builds upon available lift force results for the impermeable channel case, which enables easy generalisation to other cross-sections, (ii) this takes advantage of the fact that parameter
$\unicode[STIX]{x1D6FE}$
is naturally very small for the purposes of linearisation (since
$|\unicode[STIX]{x1D6FE}|\leqslant 0.01$
in typical microfiltration processes (Belfort, Davis & Zydney Reference Belfort, Davis and Zydney1994)), and (iii) such a linearisation has prior analytical precedent that suggests an additive decomposition of the lift force into an inertial lift force and a permeate drag force (Altena & Belfort Reference Altena and Belfort1984).
We write the linearisation as
$$\begin{eqnarray}\displaystyle \boldsymbol{F}(\unicode[STIX]{x1D6FE},Re,a/W,x/W,z/W) & {\sim} & \displaystyle \boldsymbol{F}_{0}(Re,a/W,x/W,z/W)|_{\unicode[STIX]{x1D6FE}=0}\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D6FE}\boldsymbol{F}_{1}(Re,a/W,x/W,z/W)+O(\unicode[STIX]{x1D6FE}^{2}),\end{eqnarray}$$
where
$\boldsymbol{F}_{0}$
is the lift force calculated for a particle at a given location in the absence of any permeate flow (i.e.
$\unicode[STIX]{x1D6FE}=0$
, which is the standard inertial migration in an impermeable channel) and
$\boldsymbol{F}_{1}$
represents the first-order linearisation effect. As indicated earlier, for small
$\unicode[STIX]{x1D6FE}$
, we speculate that
$\boldsymbol{F}_{1}$
is proportional to the drag force,
$F_{d}=3\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}aU$
:
$$\begin{eqnarray}\boldsymbol{F}_{1}=\boldsymbol{g}(Re,a/W,x/W,z/W)\times F_{d},\end{eqnarray}$$
where
$\boldsymbol{g}(Re,a/W,x/W,z/W)$
encodes the spatial variation across the cross-section. In this work, we do not try to analytically identify the form of the scaling function
$\boldsymbol{g}$
, but explore how it can be constructed by using a minimal set of full-physics simulations. In fact, we show below that
$\boldsymbol{g}(Re,a/W,x/W,z/W)$
(and hence,
$\boldsymbol{F}_{1}$
) can be reliably constructed using just two full-physics simulations. We rewrite (3.1) (for a fixed
$Re$
and
$a/W$
) to compute
$\boldsymbol{F}_{1}$
as
$$\begin{eqnarray}\boldsymbol{F}_{1}(Re,a/W,x/W,z/W)=\frac{\boldsymbol{F}(Re,a/W,\unicode[STIX]{x1D6FE},x/W,z/W)-\boldsymbol{F}_{0}(Re,a/W,x/W,z/W)}{\unicode[STIX]{x1D6FE}}.\end{eqnarray}$$
Figure 5. (a) Force fields for various values of
$\unicode[STIX]{x1D6FE}$
for both the DNS and the LM with a streamline trace overlaid in blue for visualisation (
$a/W=0.05$
and
$Re=100$
). The bounding box for each field is
$0.45W\times 0.45W$
. (b) The
$x$
-equilibrium plotted against the relative permeate flow
$\unicode[STIX]{x1D6FE}$
. (c) The
$z$
-equilibrium plotted against the relative permeate flow
$\unicode[STIX]{x1D6FE}$
. The discrete points represent the results from the DNS and the continuous black line is the result of the LM.
Here
$\boldsymbol{F}-\boldsymbol{F}_{0}$
is the difference between the calculated forces maps, where
$\boldsymbol{F}$
represents the force vector field of a particle in the presence of a permeate flow at some fixed
$Re$
and
$\unicode[STIX]{x1D6FE}$
and
$\boldsymbol{F}_{0}$
is the force vector field for the same particle at the same
$Re$
but in the absence of permeate flow (
$\unicode[STIX]{x1D6FE}=0$
). Figure 4(a) shows the spatial distribution of this residual for various values of
$\unicode[STIX]{x1D6FE}$
. Interestingly, it appears that these curves are self-similar when scaled by
$F_{P}=3\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}a\unicode[STIX]{x1D6FE}U$
, suggesting the utility of the proposed linear model. In figure 4(b,c) we construct
$\boldsymbol{g}$
by dividing the residual
$\boldsymbol{F}-\boldsymbol{F}_{0}$
by the characteristic permeate force
$F_{P}$
. In figure 4(b), keeping
$Re$
constant, we see three distinct curves representing our three particle sizes, showing that all normalised residual curves fall on top of each other when divided by
$\unicode[STIX]{x1D6FE}$
. Similarly, the curves collapse as shown in figure 4(c) corresponding to each value of
$Re$
(at a fixed
$a/W$
). We note that
$\boldsymbol{g}$
is insensitive to the choice of
$Re$
(appendix B).
We compare the force maps constructed using
$\boldsymbol{F}_{1}$
(i.e. linearised model, LM) with that from the full-physics simulations (i.e. direct numerical simulation, DNS) in figure 5(a). We chose the interesting case of
$Re=100$
and
$a/W=0.05$
, and seven discrete values of
$\unicode[STIX]{x1D6FE}$
, where fields yield drastically different ‘streamline’ morphologies, from a complete wall suction at
$\unicode[STIX]{x1D6FE}=0.003$
to centre-plane focusing at
$\unicode[STIX]{x1D6FE}<-0.003$
.
First we compute
$\boldsymbol{F}_{1}$
according to (3.3) for each of the six permeate solutions (i.e.
$\unicode[STIX]{x1D6FE}=\pm 0.001,\pm 0.002,\pm 0.003$
). With this result, we average all six
$\boldsymbol{F}_{1}$
fields to produce a master solution, and thus eliminate any numerical noise that is inherent in this type of model construction. From figure 5, it is apparent that the LM reconstructs the lateral lift force maps qualitatively well with little discernible error; the advantage of the LM is that it only requires knowledge of a finite number of simulations. More importantly, the LM is not limited to discrete values of
$\unicode[STIX]{x1D6FE}$
, but can be used to evaluate forces for continuous values of
$\unicode[STIX]{x1D6FE}$
. We next compare the predicted location of the focusing points (force equilibrium points) in figure 5(b,c) for different
$\unicode[STIX]{x1D6FE}$
. In figure 5(b) we see that the
$x$
-equilibrium diagram spans a wide range of the channel where the limits are set by the centreline (
$x/W=0$
) and by contact with the porous wall (
$x/W=0.475$
). Similarly, figure 5(c) shows the
$z$
-equilibrium diagram for the same system, where the equilibrium shift is seemingly independent of
$\unicode[STIX]{x1D6FE}$
. The LM provides the ability to tune the permeate flow rate
$\unicode[STIX]{x1D6FE}$
to precisely select a desired equilibrium location. This is potentially very useful for particle separation applications.
Figure 6. (a) The
$x$
-equilibrium as a function of the relative permeate force (
$F_{P}/F_{L}$
) for three particle diameters (
$a/W=0.05,0.10$
and 0.15) at
$Re=100$
. Of note is that the data are limited by either the centreline (
$x/W=0$
) or the confining walls (for large values of
$x/W=0.50$
). (b) The
$z$
-equilibrium diagram for three particle diameters (
$a/W=0.05,0.10$
and
$0.15$
). Here the equilibrium shift is less sensitive than the
$x_{eq}$
counterpart, that is, the particle deviates only slightly from the zero permeate equilibrium. (c) The
$x$
-equilibrium as a function of the relative permeate force (
$F_{P}/F_{L}$
) for four
$Re$
(25, 50, 75 and 100) for a single particle of diameter
$a/W=0.05$
. (d) The
$z$
-equilibrium as a function of the relative permeate force (
$F_{P}/F_{L}$
) for four
$Re$
(25, 50, 75 and 100) for a single particle of diameter
$a/W=0.05$
.
We next investigate the scaling interplay between the inertial force
$\boldsymbol{F}_{0}$
and the permeate drag
$\unicode[STIX]{x1D6FE}\boldsymbol{F}_{1}$
with increasing
$Re$
and particle size. We note that
$\unicode[STIX]{x1D6FE}\boldsymbol{F}_{1}$
scales as the Stokes drag force, i.e.
$F_{P}\sim 3\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}a\unicode[STIX]{x1D6FE}U$
, and the inertial lift forces scale as
$F_{L}\sim \unicode[STIX]{x1D70C}U^{2}a^{4}/W^{2}$
. Hence, when
$|F_{P}|\gg F_{L}$
we would expect that permeate forces are dominant and the force fields would appear as if inertial forces were negligible. Conversely, when
$|F_{P}|\ll F_{L}$
we would expect that inertial forces would dictate the migration and focusing of a particle in the channel. In figure 6 we plot the equilibrium location against the relative permeate force
$F_{P}/F_{L}$
for all particle sizes (figure 6
a,b) and
$Re$
(figure 6
c,d). The continuous diagrams seen in figure 6 show very clear and specific trends that would be difficult to interpret from discrete DNS data such as that seen in figure 3. The equilibrium diagrams shown in figure 6(a,c) clearly demonstrate the behaviour of a particle at the two extremes as predicted. That is, for large values of magnitude of
$|F_{P}|\gg F_{L}$
we see a dominance of the permeate flow and the force equilibrium is shifted either to the walls (
$x_{eq}/W=(W-a)/2W$
) or to the centreline (
$x_{eq}=0$
); on the other extreme, when
$|F_{P}|\ll F_{L}$
, we see that the force equilibrium is only slightly perturbed from the case of no permeate flow. Finally, figure 6(b,d) shows the change in
$z$
-equilibrium location as a function of the relative permeate force (
$F_{P}/F_{L}$
) for
$a/W$
and
$Re$
. The
$z_{eq}$
of large particles are more affected when compared to small particles under this representation of
$F_{P}/F_{L}$
(figure 6
b). Further, in figure 6(d) we see that as
$F_{P}/F_{L}$
increases, the
$z_{eq}$
migrates towards the channel wall for low
$Re$
and towards the centreline for high
$Re$
. This behaviour is indicative of secondary effects that cannot be captured by such a simple scaling argument (Di Carlo et al.
Reference Di Carlo, Edd, Humphry, Stone and Toner2009).
3.5 Limits of the linear model
We next explore when the additive decomposition of the force into an inertial component and a linear viscous component breaks down. We examine the error in the LM by comparing the lateral forces of the DNS with those constructed with our LM to determine when nonlinear effects due to
$\unicode[STIX]{x1D6FE}$
,
$Re$
and
$a/W$
cause the model to become unreliable. In figure 7(a) we show a comparison of the inertial force for three distinct
$\unicode[STIX]{x1D6FE}$
at
$Re=100$
and
$a/W=0.10$
and in figure 7(b) for three distinct
$Re$
for
$\unicode[STIX]{x1D6FE}=0.2$
and
$a/W=0.10$
. As expected, for increasing values of both
$\unicode[STIX]{x1D6FE}$
and
$Re$
, mismatch between the two models increases. However, the increasing mismatch associated with a change in
$Re$
is difficult to discern and so we plot in the inset to figure 7(b) the relative spatial error for each case. We define the error in the LM as
$$\begin{eqnarray}\text{error}=\frac{\Vert \boldsymbol{F}_{DNS}-\boldsymbol{F}_{LM}\Vert _{2}}{\Vert \boldsymbol{F}_{DNS}\Vert _{2}},\end{eqnarray}$$
where
$\boldsymbol{F}_{DNS}$
and
$\boldsymbol{F}_{LM}$
are the force distribution calculated using the DNS and the LM respectively. Figure 7(b) shows how the error in the LM increases with
$a/W$
,
$\unicode[STIX]{x1D6FE}$
and
$Re$
. In general, as the permeate Reynolds number
$Re_{U_{W}}$
(
$Re_{U_{W}}=|\unicode[STIX]{x1D6FE}|(a/W)Re$
) increases, so does error. This confirms our hypothesis that as both inertial and permeate flow increase, nonlinear effects become more prominent and cannot be captured in our linear model.
Over all sets of data for
$Re_{U_{W}}\leqslant 1$
the error in the LM is less than 5 %, represented by the dashed lines in figure 7(b). Therefore, in figure 7(d) we determine the maximum
$|\unicode[STIX]{x1D6FE}|$
for a given
$(a/W)Re$
defined as where error will be 5 %. Superimposed on this figure we show the
$\unicode[STIX]{x1D6FE}$
(at a given
$(a/W)Re$
) at which the particle equilibrium location coincides with the confining walls (i.e. wall contact) and therefore an increase in
$\unicode[STIX]{x1D6FE}$
will no longer result in a change in particle equilibrium. It is clear that for the particles studied in this work, we will never reach the 5 % error limit of the LM, suggesting that the LM is a viable, useful and applicable model for fast experimental exploration of the effect inertial migration under permeate flow conditions.
Figure 7. A direct comparison of the forces on a particle as computed by a DNS and the LM showing how error in the LM increases with increasing
$\unicode[STIX]{x1D6FE}$
(a) and
$Re$
(b). It is difficult to see from a direct comparison that the error is increasing with
$Re$
; therefore, we plot in the inset the local relative error. (c) A plot illustrating how error in the LM increases with the permeate Reynolds number
$Re_{U_{W}}=|\unicode[STIX]{x1D6FE}|Re(a/W)$
. The dashed lines show that for
$Re_{U_{W}}>1$
, error in the LM is greater than 5 %. (d) The
$\unicode[STIX]{x1D6FE}-Re(a/W)$
operating parameter space for the LM, where using the LM in the space beneath the black dashed line should yield less than 5 % error in the model. Superposed onto this operating parameter space are data points for the three particle sizes representing the value of
$|\unicode[STIX]{x1D6FE}|$
at which complete wall suction occurs for the particles studied.
Figure 8. (a) Schematic showing the microfluidic device composed of a long straight region (
${\approx}1.9$
cm) followed by a region (
$L=1.0$
cm) where there is an array of perpendicular permeate channels (inset) that allow permeate flow to enter or exit the channel. (b) Long-exposure image of
$10~\unicode[STIX]{x03BC}\text{m}$
fluorescent polystyrene particles in a
$100\times 100~\unicode[STIX]{x03BC}\text{m}$
(
$W\times W$
) channel. The in-plane particles are measured at a distance
$x$
relative to the centreline. (c) Comparison of computed (solid lines) and experimentally measured (shapes) trajectories of a particle
$a/W=0.10$
for various operating conditions of the microfluidic device.
4 Experiments
In this section, we illustrate the utility of the LM by comparing experiments with predicted trajectories of inertial particles within a porous channel. We microfabricated a model porous system as shown in figure 8(a). The microfluidic device is composed of a primary channel that is 3.2 cm long with a square cross-sectional area of
$100\times 100~\unicode[STIX]{x03BC}\text{m}$
. The particles enter the channel through a long straight region intended to prefocus the particles. Following this region, there is a permeate region
$L=1.0$
cm where flow can enter and exit the channel through an array of permeate channels of
$W_{P}=5~\unicode[STIX]{x03BC}\text{m}$
width spaced
$\unicode[STIX]{x1D6FF}=50~\unicode[STIX]{x03BC}\text{m}$
apart, with a length of
$L_{P}=4.95$
mm. A two-syringe pump system (Harvard Apparatus) provides both inlet flow and permeate flow: the first pump infuses the inlet flow into the primary channel at a constant volumetric flow rate (
$Q_{F}$
), and the second pump has two possible configurations depending on the ratio
$\unicode[STIX]{x1D6FD}=Q_{R}/Q_{F}$
. If
$\unicode[STIX]{x1D6FD}$
is less than one the second pump is placed at the exit of the primary channel and limits the flow rate to a rate of
$Q_{R}$
through the primary outlet. If
$\unicode[STIX]{x1D6FD}$
is greater than one the second pump infuses flow through the permeate channels at a rate
$Q_{P}$
where
$Q_{F}+Q_{P}=Q_{R}$
. In general the flow in a device like this is highly dependent on the relative hydrodynamic resistance of the permeate channels to that of the main channel. In our geometry, the permeate resistance is large compared to that of the main channel and it can be shown that in this limit the volumetric flow rate decreases linearly with axial location, so
$Q(y)=Q_{F}[1+(y/L)(\unicode[STIX]{x1D6FD}-1)]$
(appendix C). Using conservation of mass we find the constant permeate velocity to be
$U_{W}=(Q_{F}(1-\unicode[STIX]{x1D6FD}))/2WL$
, and thus
$\unicode[STIX]{x1D6FE}(y)=W^{2}U_{W}/Q(y)$
.
We performed migration experiments with a suspension of
$10~\unicode[STIX]{x03BC}\text{m}$
fluorescent polystyrene particles in deionised water at a concentration of
$10^{4}$
particles per millilitre, and we add 0.5 % (
$v/v$
) Tween 20 (Sigma-Aldrich) to reduce particle aggregation. To find the particle locations, we record streak images with a CCD camera (Andor Luca) by accumulating approximately 25 s of image data at each downstream location and the post-processing using image-processing software (MATLAB). Figure 8(b) shows an example of a streak image. We post-process images like this to locate three peaks: the centre peak corresponds to the out-of-plane equilibrium while the remaining two correspond to the in-plane equilibrium. Since the in-plane equilibria are symmetric about the centreline, we build a trajectory from one in-plane equilibrium at various axial locations along the length of the channel to compare with our LM. More details of the experiments can be found in Garcia & Pennathur (Reference Garcia and Pennathur2017).
Figure 8(c) shows four experimental trajectories of
$10~\unicode[STIX]{x03BC}\text{m}$
particles. The prefocused particles enter the permeate region at
$y/L=0$
and begin to migrate either towards the wall or away, depending on the direction of the permeate flow. Here values of
$\unicode[STIX]{x1D6FD}>1$
indicate the permeate flow is directed towards the centre of the channel and conversely
$\unicode[STIX]{x1D6FD}<1$
indicates that the flow is directed out of the channel. Even though the permeate flow is constant, the trajectories can be non-monotonic, where the particles begin migrating in one direction and subsequently reverse directions due to the evolving nature of inertial and permeate forces in the spatially varying flow field (figure 8
c). We can calculate the particle trajectories with our linear model by knowing how the flow parameters (i.e.
$\unicode[STIX]{x1D6FE}$
and
$Re$
) change with axial location in the channel.
We calculated the theoretical trajectories with a first-order time-stepping approximation (i.e. Euler method), evaluating the force at
$z=0$
. We do not expect much motion in the
$z$
-direction because the lift forces in the
$z$
-direction are much smaller than the
$x$
-direction lift force for the prefocused stream of particles we are modelling. Thus the equations of motion become
$$\begin{eqnarray}\displaystyle & \displaystyle y_{n+1}=y_{n}+u_{y}(x_{n},z=0)\frac{Q(y)}{Q_{F}}\unicode[STIX]{x0394}t, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle x_{n+1}=x_{n}+\frac{F_{x}(x_{n},y_{n},z=0)}{3\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}a\unicode[STIX]{x1D706}}\unicode[STIX]{x0394}t, & \displaystyle\end{eqnarray}$$
where
$u_{y}(x,z)$
is the flow field in a square channel (White Reference White2006) and
$F_{x}(x,y,z)$
is the predicted force in the lateral direction calculated using our linear model. The axial dependence of the lateral force
$F_{x}$
is calculated by mapping the axial location of the particle to a corresponding local value of
$Re$
. With this information, we interpolate between a pre-calculated dataset and generate the local zero permeate lift force. Using the observation that
$g_{x}$
is invariant to
$Re$
(figure 4
c), we construct the linear model using (3.1). Finally,
$\unicode[STIX]{x1D706}$
is the correction factor for the Stokes drag near a confining wall (appendix B). Using the linear model we are able to reproduce the experimental trajectories reasonably well (figure 8
c), further proving the viability of our linear model.
5 Conclusion
Our findings describe the spatially inhomogeneous forces on confined inertial particles in the presence of a permeate flow. Our numerical simulations suggest that the relative permeate force (
$F_{P}/F_{L}$
) is an important parameter in the characterisation of behaviour of these particles. For very small magnitudes of the relative permeate force, the location of force equilibrium remains unchanged and reminiscent of flow in a non-porous duct. As the magnitude of the relative permeate force is increased, the equilibrium position deviates further from the non-porous case until it is limited by either the wall or centreline. Using the results from our numerical simulations we are able to construct a model which superposes the linear viscous effects of the permeate flow to that of the underlying inertial forces. This linear model shows excellent agreement over a continuous span of permeate flow with both full simulations and experimental observations when
$Re_{U_{W}}<1$
with no added computational penalty. This is especially noteworthy because the flow field in a channel with permeable walls is not trivial and would normally require much effort to simulate with another approach. We speculate that this model can help in the rapid design of microfluidic devices that can precisely manipulate particle streams. Furthermore, the framework for our linear model presented in this work can also be implemented in other systems where external forcing of inertial particles exists, such as with magnetic or electric forces. Our model greatly reduces the complexity of a well-studied and ubiquitous flow. We emphasise the fact that the linearisation, while useful, does not account for potentially important phenomena like particle deformability and non-Newtonian suspending fluids, which are typical of real-world systems. Developing similar simple and computationally efficient methods for accounting for such interactions is a future area of research work.
Acknowledgements
We would like to acknowledge Professor B. Bamieh for enlightening discussions on nonlinear bifurcating systems. M.G. was supported by the Institute for Collaborative Biotechnologies through grant nos. W911NF-09-D-0001 and W911NF-12-1-0031 through the US Army Research Office. The content of the information does not necessarily reflect the position of the policy of the government, and no official endorsement should be inferred.
Appendix A. Additional momentum term
The momentum equation can be transformed from a velocity field in the laboratory frame to that of a frame translating with a particle moving at a velocity
$U_{p}$
using the following transformations:
$$\begin{eqnarray}\displaystyle & t^{\prime }=t, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & x_{i}^{\prime }=x_{i}+x_{p,i}(t), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & u_{i}^{\prime }=u_{i}+U_{p,i}(t), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & p^{\prime }=p. & \displaystyle\end{eqnarray}$$
Here the variables associated with a prime denote the laboratory frame and those without a prime denote the moving frame. We can then relate the derivatives with respect to time and space for either frame by
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}(\,)}{\unicode[STIX]{x2202}x_{j}}=\frac{\unicode[STIX]{x2202}x_{k}^{\prime }}{\unicode[STIX]{x2202}x_{j}}\frac{\unicode[STIX]{x2202}(\,)}{\unicode[STIX]{x2202}x_{k}^{\prime }}=\frac{\unicode[STIX]{x2202}(x_{k}+x_{p,k})}{\unicode[STIX]{x2202}x_{j}}\frac{\unicode[STIX]{x2202}(\,)}{\unicode[STIX]{x2202}x_{k}^{\prime }}=\unicode[STIX]{x1D6FF}_{jk}\frac{\unicode[STIX]{x2202}(\,)}{\unicode[STIX]{x2202}x_{k}^{\prime }}=\frac{\unicode[STIX]{x2202}(\,)}{\unicode[STIX]{x2202}x_{j}^{\prime }}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}(\,)}{\unicode[STIX]{x2202}t}=\frac{\unicode[STIX]{x2202}t^{\prime }}{\unicode[STIX]{x2202}t}\frac{\unicode[STIX]{x2202}(\,)}{\unicode[STIX]{x2202}t^{\prime }}+\frac{\unicode[STIX]{x2202}x_{k}^{\prime }}{\unicode[STIX]{x2202}t}\frac{\unicode[STIX]{x2202}(\,)}{\unicode[STIX]{x2202}x_{k}^{\prime }}=\frac{\unicode[STIX]{x2202}(\,)}{\unicode[STIX]{x2202}t^{\prime }}+U_{p,k}\frac{\unicode[STIX]{x2202}(\,)}{\unicode[STIX]{x2202}x_{k}^{\prime }}. & \displaystyle\end{eqnarray}$$
We can then relate the momentum equation in the laboratory frame to that of a moving frame by using (A 5) and (A 6):
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}\left[\frac{\unicode[STIX]{x2202}u_{i}^{\prime }}{\unicode[STIX]{x2202}t^{\prime }}+u_{j}^{\prime }\frac{\unicode[STIX]{x2202}u_{i}^{\prime }}{\unicode[STIX]{x2202}x_{j}^{\prime }}\right]=-\frac{\unicode[STIX]{x2202}p^{\prime }}{\unicode[STIX]{x2202}x_{i}^{\prime }}+\unicode[STIX]{x1D707}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}^{\prime }}\left(\frac{\unicode[STIX]{x2202}u_{i}^{\prime }}{\unicode[STIX]{x2202}x_{j}^{\prime }}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}\left[\frac{\unicode[STIX]{x2202}(u_{i}+U_{p,i})}{\unicode[STIX]{x2202}t^{\prime }}+(u_{j}+U_{p,j})\frac{\unicode[STIX]{x2202}(u_{i}+U_{p,i})}{\unicode[STIX]{x2202}x_{j}^{\prime }}\right]=-\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x_{i}^{\prime }}+\unicode[STIX]{x1D707}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}^{\prime }}\left(\frac{\unicode[STIX]{x2202}(u_{i}+U_{p,i})}{\unicode[STIX]{x2202}x_{j}^{\prime }}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle \unicode[STIX]{x1D70C}\left[\frac{\unicode[STIX]{x2202}(u_{i}+U_{p,i})}{\unicode[STIX]{x2202}t}-U_{p,j}\frac{\unicode[STIX]{x2202}(u_{i}+U_{p,i})}{\unicode[STIX]{x2202}x_{j}}+(u_{j}+U_{p,j})\frac{\unicode[STIX]{x2202}(u_{i}+U_{p,i})}{\unicode[STIX]{x2202}x_{j}}\right]\nonumber\\ \displaystyle & & \displaystyle \quad =-\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x_{i}}+\unicode[STIX]{x1D707}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}}\left(\frac{\unicode[STIX]{x2202}(u_{i}+U_{p,i})}{\unicode[STIX]{x2202}x_{j}}\right).\end{eqnarray}$$
Now cancelling like terms and knowing that
$U_{p,i}$
has no spatial gradients because it is only a linear translation of the laboratory frame and not a continuum value, we can write the momentum equation as
$$\begin{eqnarray}\unicode[STIX]{x1D70C}\left[\frac{\unicode[STIX]{x2202}(u_{i}+U_{p,i})}{\unicode[STIX]{x2202}t}+u_{j}\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{j}}\right]=-\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x_{i}}+\unicode[STIX]{x1D707}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}}\left(\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{j}}\right).\end{eqnarray}$$
We assume that in the moving reference frame the flow is quasi-steady and thus the time derivative of the flow is zero, yielding
$$\begin{eqnarray}\unicode[STIX]{x1D70C}\left[\frac{\unicode[STIX]{x2202}U_{p,i}}{\unicode[STIX]{x2202}t}+u_{j}\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{j}}\right]=-\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x_{i}}+\unicode[STIX]{x1D707}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}}\left(\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{j}}\right).\end{eqnarray}$$
At a given moment in time the acceleration of the particle that is being tracked by the moving reference frame is dictated by the underlying flow (
$\bar{u}_{i}$
) and its gradient and thus the acceleration is given by
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}U_{p,i}}{\unicode[STIX]{x2202}t}=U_{p,j}\frac{\unicode[STIX]{x2202}\bar{u}_{i}}{\unicode[STIX]{x2202}x_{j}}.\end{eqnarray}$$
Finally, substituting this term back into the momentum equation yields
$$\begin{eqnarray}\unicode[STIX]{x1D70C}\left[U_{p,j}\frac{\unicode[STIX]{x2202}\bar{u}_{i}}{\unicode[STIX]{x2202}x_{j}}+u_{j}\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{j}}\right]=-\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x_{i}}+\unicode[STIX]{x1D707}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}}\left(\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{j}}\right),\end{eqnarray}$$
where
$\bar{u}$
is the undisturbed velocity field. However, the only non-zero component of (A 12) in our model acts in the axial direction and scales with
$\unicode[STIX]{x1D6FE}$
as
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}U_{p,i}}{\unicode[STIX]{x2202}t}=U_{p,j}\frac{\unicode[STIX]{x2202}\bar{u}_{i}}{\unicode[STIX]{x2202}x_{j}}\sim U_{p}\frac{\unicode[STIX]{x1D6FE}U}{W}.\end{eqnarray}$$
Note that the term is negligible for small
$\unicode[STIX]{x1D6FE}$
and near the channel walls where
$U_{p}$
approaches zero.
Appendix B. Correction to the Stokes drag and the relationship to g
When the size of a flowing particle is comparable to the dimensions of its confining channel, it is important to consider the hydrodynamic effects from the walls on the lateral motion of the particle. To account for this, we introduce a spatially varying retardation factor
$\unicode[STIX]{x1D706}$
that relates force
$F$
on a particle of diameter
$a$
and its migration velocity
$U$
(Brenner Reference Brenner1961):
$$\begin{eqnarray}F=3\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}aU\unicode[STIX]{x1D706}.\end{eqnarray}$$
Figure 9. (a) The relationship between the translational velocity of a sphere (
$U$
) and the force (
$F$
) required to produce such a motion can be significantly modified by the effects of confining walls. We show the effect of this spatial retardation (
$\unicode[STIX]{x1D706}$
) on spheres of diameter
$a/W=0.05,0.10,0.15$
(at
$z/W=0$
) and compare with the analytic results of Brenner (Reference Brenner1961). (b) A plot of the spatially varying permeate velocity (at
$z/W=0$
) within the channel modelled in this work. At the walls the permeate velocity is maximal (
$u_{x}/U_{W}=1$
) and decays to zero at the centre of the channel. (c) Normalised residual curves (
$g_{x}$
) as a function of
$x/W$
at
$z/W=0$
for a particle of diameter
$a/W=0.05$
. (d) Normalised residual curves (
$g_{x}$
) as a function of
$x/W$
at
$z/W=0$
for a particle of diameter
$a/W=0.10$
. (e) Normalised residual curves (
$g_{x}$
) as a function of
$x/W$
at
$z/W=0$
for a particle of diameter
$a/W=0.15$
. The solid black lines in (c–e) represent
$g_{x}$
modelled using the retardation factor and permeate velocity depicted in (a) and (b) respectively.
Here,
$\unicode[STIX]{x1D706}$
is calculated with a numerical simulation where a particle is modelled in a square channel of cross-section
$W\times W$
and is assigned a velocity
$U$
(
$Re\ll 1$
) in the direction of a wall. Then we calculate the required force to produce such a motion by integrating the fluid stress on the surface of the particle and using (B 1) we solve for
$\unicode[STIX]{x1D706}$
. We repeat this process for discrete locations that span the width of the channel. In figure 9(a) we plot
$\unicode[STIX]{x1D706}$
for three particles (
$a/W=0.05,0.10,0.15$
at
$z/W=0$
) and show how our calculations compare to the existing analytic model (Brenner Reference Brenner1961). Interestingly, even though in our simulations the particle is confined by four walls, the semi-infinite domain of the Brenner model seems to match well at
$z/W=0$
. In our work we find that the normalised residual
$\boldsymbol{g}$
(figure 4
b,c) is weakly dependent on
$Re$
. A quick approximation of
$\boldsymbol{g}$
(for a particle of any size) may useful for any researcher who may be interested in applying this linear model. Here we provide a simple method for approximating
$\boldsymbol{g}$
. We begin by interpreting
$\boldsymbol{g}$
to represent the effects of the permeate flow (
$U_{W}$
) on the force (
$\boldsymbol{F}_{P}$
) experienced by a particle:
$$\begin{eqnarray}\boldsymbol{F}_{P}=\unicode[STIX]{x1D6FE}\boldsymbol{F}_{1}=3\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}aU_{W}\boldsymbol{g}.\end{eqnarray}$$
Likewise, if we interpret (B 1) as the force experienced by particle with a flow moving past it (i.e. particle reference frame), then we can construct a force field everywhere in the channel according to the local permeate velocity given by the
$x$
- and
$z$
-components of undisturbed flow
$\bar{\boldsymbol{u}}$
that should approximate the permeate force field:
$$\begin{eqnarray}\boldsymbol{F}_{p}\approx 3\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}a\bar{\boldsymbol{u}}\unicode[STIX]{x1D706}.\end{eqnarray}$$
Without loss of generality we consider only the
$x$
-component of the permeate force. If we equate (B 2) and (B 3) we can then use this relationship to obtain
$g_{x}$
:
$$\begin{eqnarray}g_{x}\approx \unicode[STIX]{x1D706}\frac{\bar{u}_{x}}{U_{W}}.\end{eqnarray}$$
Here
$\bar{u}_{x}$
is the
$x$
-component of the undisturbed flow and can be seen in figure 9(b). The results of implementing (B 4) can be seen in figure 9(c–e). In general the approximation works best for smaller particles at a lower
$Re$
. However, the approximation is remarkably accurate and only requires knowledge of flow field, which is relatively simple to determine.
Appendix C. Local volumetric flow rate in a porous channel
The local volumetric flow rate in a porous duct
$Q(y)$
is typically axially varying due to a non-zero fluid flux. Thus the local volumetric flow rate can be modelled as
$$\begin{eqnarray}Q(y)=Q_{F}-\int _{0}^{y}U_{W}(y^{\prime })W\,\text{d}y^{\prime },\end{eqnarray}$$
where
$Q_{F}$
is the feed flow rate, i.e.
$Q(y=0)=Q_{F}$
,
$U_{W}$
is the average permeate flow velocity,
$W$
is the channel height and
$y$
is the axial coordinate. Now assuming that the wall can be treated as a continuously porous wall of a constant permeability, we can use Darcy’s law to write
$$\begin{eqnarray}Q(y)=Q_{F}-\int _{0}^{y}\unicode[STIX]{x1D705}\frac{P(y^{\prime })-P_{0}}{\unicode[STIX]{x1D707}L_{P}}W\,\text{d}y^{\prime },\end{eqnarray}$$
where
$P(y^{\prime })$
is the local pressure in the channel,
$P_{0}$
is the reservoir pressure on the other side of the porous wall,
$L_{p}$
is the wall thickness and
$\unicode[STIX]{x1D705}$
is the wall permeability. It is convenient to define the constant
$R_{P}=\unicode[STIX]{x1D707}L_{P}/\unicode[STIX]{x1D705}W$
so that we can rewrite (C 2) as
$$\begin{eqnarray}Q(y)=Q_{F}-\int _{0}^{y}\frac{P(y^{\prime })-P_{0}}{R_{P}}\,\text{d}y^{\prime }.\end{eqnarray}$$
Equation (C 3) can be rewritten in terms of pressure by using the Hagen–Poiseuille equation:
$$\begin{eqnarray}Q=-\frac{1}{R_{C}}\frac{\text{d}P}{\text{d}y}.\end{eqnarray}$$
Here
$R_{C}$
represents the channel resistance per unit length in the axial direction. This equation can now be substituted into (C 3):
$$\begin{eqnarray}-\frac{1}{R_{C}}\frac{\text{d}P}{\text{d}y}=Q_{F}-\int _{0}^{y}\frac{P(y^{\prime })-P_{0}}{R_{P}}\,\text{d}y^{\prime }.\end{eqnarray}$$
If we differentiate (C 5) once and use the fundamental theorem of calculus we obtain an ordinary differential equation for the axial pressure distribution:
$$\begin{eqnarray}-\frac{1}{R_{C}}\frac{\text{d}^{2}P}{\text{d}y^{2}}=-\frac{P(y)-P_{0}}{R_{P}}.\end{eqnarray}$$
Now rearranging the equation to a standard form
$$\begin{eqnarray}\frac{R_{P}}{R_{C}}\frac{\text{d}^{2}P}{\text{d}y^{2}}-P(y)+P_{0}=0,\end{eqnarray}$$
say that
$\unicode[STIX]{x1D6FC}^{2}=R_{P}/R_{C}$
and solve the ordinary differential equation:
$$\begin{eqnarray}P(y)=C_{1}\exp \left(\frac{y}{\unicode[STIX]{x1D6FC}}\right)+C_{2}\exp \left(\frac{-y}{\unicode[STIX]{x1D6FC}}\right)+P_{0}.\end{eqnarray}$$
This solution is valid; however, the known boundary conditions are in terms of volumetric flow rate. Therefore, it is convenient to differentiate the solution once:
$$\begin{eqnarray}\frac{\text{d}P}{\text{d}y}=\frac{C_{2}}{\unicode[STIX]{x1D6FC}}\exp \left(\frac{-y}{\unicode[STIX]{x1D6FC}}\right)-\frac{C_{1}}{\unicode[STIX]{x1D6FC}}\exp \left(\frac{y}{\unicode[STIX]{x1D6FC}}\right)\end{eqnarray}$$
and then again apply the Hagen–Poiseuille relationship to obtain the axial flow rate distribution:
$$\begin{eqnarray}Q(y)=\frac{C_{1}}{\unicode[STIX]{x1D6FC}R_{C}}\exp \left(\frac{y}{\unicode[STIX]{x1D6FC}}\right)-\frac{C_{2}}{\unicode[STIX]{x1D6FC}R_{C}}\exp \left(\frac{-y}{\unicode[STIX]{x1D6FC}}\right).\end{eqnarray}$$
Here the boundary conditions are
$Q(0)=Q_{F}$
and
$Q(L)=Q_{R}$
and are used to solve for
$C_{1}$
and
$C_{2}$
:
$$\begin{eqnarray}C_{1}=\frac{-\unicode[STIX]{x1D6FC}R_{C}\exp \left({\displaystyle \frac{L}{\unicode[STIX]{x1D6FC}}}\right)\left[Q_{R}-Q_{F}\exp \left({\displaystyle \frac{-L}{\unicode[STIX]{x1D6FC}}}\right)\right]}{\exp \left({\displaystyle \frac{2L}{\unicode[STIX]{x1D6FC}}}\right)-1}\end{eqnarray}$$
and
$$\begin{eqnarray}C_{2}=\frac{-\unicode[STIX]{x1D6FC}R_{C}\exp \left({\displaystyle \frac{L}{\unicode[STIX]{x1D6FC}}}\right)\left[Q_{R}-Q_{F}\exp \left({\displaystyle \frac{L}{\unicode[STIX]{x1D6FC}}}\right)\right]}{\exp \left({\displaystyle \frac{2L}{\unicode[STIX]{x1D6FC}}}\right)-1}.\end{eqnarray}$$
Upon simplification,
$$\begin{eqnarray}Q(y)=\frac{Q_{R}\sinh \left({\displaystyle \frac{y}{\unicode[STIX]{x1D6FC}}}\right)+Q_{F}\sinh \left({\displaystyle \frac{L-y}{\unicode[STIX]{x1D6FC}}}\right)}{\sinh \left({\displaystyle \frac{L}{\unicode[STIX]{x1D6FC}}}\right)}.\end{eqnarray}$$
This model is sufficient to model the flow rate in a porous channel, but it can be further simplified in the limit where the flow resistance through the channel wall is much greater than the resistance through the primary channel
$L/\unicode[STIX]{x1D6FC}\ll 1$
and
$\sinh (L/\unicode[STIX]{x1D6FC})\approx L/\unicode[STIX]{x1D6FC}$
, and in this limit (C 13) reduces down to
$$\begin{eqnarray}Q(y)=Q_{F}\left[1+\frac{y}{L}(\unicode[STIX]{x1D6FD}-1)\right]\quad \text{and}\quad \unicode[STIX]{x1D6FD}=Q_{R}/Q_{F}.\end{eqnarray}$$
The resistance of a rectangular section of channel of height
$2b$
and width
$2a$
(
$a>b$
) is given by (White Reference White2006)
$$\begin{eqnarray}R=\frac{3\unicode[STIX]{x1D707}}{4ba^{3}\left(1-{\displaystyle \frac{192a}{\unicode[STIX]{x03C0}^{5}b}}\displaystyle \mathop{\sum }_{n=1,3,5,\ldots }^{\infty }{\displaystyle \frac{\tanh (n\unicode[STIX]{x03C0}b/2a)}{n^{5}}}\right)}.\end{eqnarray}$$
For a square channel of width
$W$
and height
$W$
, the axial channel resistance is given by
$$\begin{eqnarray}R_{C}\approx \frac{28\unicode[STIX]{x1D707}}{W^{4}}.\end{eqnarray}$$
The permeability of a porous material is defined as
$$\begin{eqnarray}\unicode[STIX]{x1D705}=\unicode[STIX]{x1D707}\bar{U}\frac{L_{p}}{\unicode[STIX]{x0394}P}.\end{eqnarray}$$
Choosing the right formulation for average permeate velocity
$\bar{U}_{W}$
is not obvious, but one reasonable approximation is to treat the wall as continuously permeable (rather than discretely) but then use the resistance of a single permeate channel to calculate the volumetric flow rate for a section of the channel:
$$\begin{eqnarray}Q=2\bar{U}_{W}W(W_{P}+\unicode[STIX]{x1D6FF})=\left[1-\frac{192W}{\unicode[STIX]{x03C0}^{5}W_{p}}\mathop{\sum }_{n=1,3,5,\ldots }^{\infty }\frac{\tanh (n\unicode[STIX]{x03C0}W_{p}/2W)}{n^{5}}\right]\frac{W^{3}W_{P}}{12\unicode[STIX]{x1D707}}\frac{\unicode[STIX]{x0394}P}{L_{p}}.\end{eqnarray}$$
Figure 10. (a) Schematic showing the microfluidic device used in the experiments. The channel has a height
$W$
into the page. (b) Volumetric flow rate curves for channels of varying channel resistance
$L/\unicode[STIX]{x1D6FC}$
(
$\unicode[STIX]{x1D6FC}^{2}=R_{P}/R_{C}$
) for
$\unicode[STIX]{x1D6FD}=0.1$
(
$\unicode[STIX]{x1D6FD}=Q_{R}/Q_{F}$
). As
$L/\unicode[STIX]{x1D6FC}$
approaches 0, the axial flow rate distribution becomes linear.
Equation (C 17) can be used to solve for an expression giving
$\bar{U}_{W}$
. Here
$\unicode[STIX]{x0394}P=P(y)-P_{0}$
and
$W_{P}$
is the width of the permeate channels (figure 10
a). Knowing
$\bar{U}_{W}$
we can substitute that into the definition of
$\unicode[STIX]{x1D705}$
(C 16) then substitute
$\unicode[STIX]{x1D705}$
into the definition of
$R_{P}$
to arrive at
$$\begin{eqnarray}R_{P}\approx \frac{24\unicode[STIX]{x1D707}L_{P}(W_{p}+\unicode[STIX]{x1D6FF})}{W^{3}W_{P}[1-0.63(W/W_{P})\tanh (\unicode[STIX]{x03C0}W_{p}/2W)]}.\end{eqnarray}$$
If we use values that are representative of our microchannel we can calculate a value of
$\unicode[STIX]{x1D6FC}\approx 0.0194$
m and the resulting ratio
$L/\unicode[STIX]{x1D6FC}=0.52$
. This ratio is of
$O(10^{-1})$
which is sufficiently small to use the small-angle approximation that results in the distribution of (C 14). Figure 10(b) shows that as
$L/\unicode[STIX]{x1D6FC}$
decreases we observe a more linear trend in the axial flow rate distribution, and for
$L/\unicode[STIX]{x1D6FC}=0.01$
the distribution is for all practical purposes linear which is approximated well by
$L/\unicode[STIX]{x1D6FC}=0.52$
.
