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

(3.1) $$\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$ :

(3.2) $$\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

(3.3) $$\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

(3.4) $$\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.