Abstract

Background

The flow of suspensions through bifurcations is encountered in several applications. It is known that the partitioning of particles at a bifurcation is different from the partitioning of the suspending fluid, which allows particle separation and fractionation. Previous works have mainly investigated the dynamics of particles suspended in Newtonian liquids.

Methods

In this work, we study through 2D direct numerical simulations the partitioning of particles suspended in non-Newtonian fluids flowing in a T-junction. We adopt a flow configuration such that the two outlets are orthogonal, and their flow rates can be tuned. A fictitious domain method combined with a grid deformation procedure is used. The effect of fluid rheology on the partitioning of particles between the two outlets is investigated by selecting different constitutive equations to model the suspending liquid. Specifically, an inelastic shear-thinning (Bird-Carreau) and a viscoelastic shear-thinning (Giesekus) models have been chosen; the results are also compared with the case of a Newtonian suspending liquid.

Results

Simulations are carried out by varying the confinement, the inlet flow rate and the relative weight of the two outlet flow rates. For each condition, the fluxes of particles through the two outflow channels are computed. The results show that shear-thinning does not have a relevant effect as compared to the equivalent Newtonian case, i.e., with the same choice of the relative outlet flow rates. On the other hand, fluid elasticity strongly alters the fraction of particles exiting the two outlets as compared to the inlet. Such effect is more pronounced for larger particles and inlet flow rates.

Conclusions

The results illustrated here show the feasibility to efficiently separate/fractionate particles by size, through the use of viscoelastic suspending liquids.

Keywords

Background

The flow of suspensions through channel bifurcations is a relevant topic in several applications such as microfluidics and biology. Indeed, a number of microfluidic devices consisting of a main channel with several side branches, aimed at separating particles with different size, has been proposed [1–4]. From the biological side, the microcirculatory network is made of many capillary bifurcations. The partitioning of red blood cells, white blood cells and platelets at vessel junctions plays an important role in determining the microvessel hematocrits and has important physiological consequences [5–9]. It is not surprising, then, that a vast literature addressing the problem of microcirculation in capillary networks is available [10–12].

Experimental and numerical studies have demonstrated that, when a suspension of rigid, non-Brownian particles under inertialess conditions flows through an asymmetric bifurcation, the partitioning of particles between the two downstream branches is different from the partitioning of the total suspension between these branches [13–15]. Volume excluded effects in the inlet channel, due to the existence of a zone near the channel wall unaccessible to the particles, and particle-wall hydrodynamic interactions are responsible for this different partitioning at bifurcations [16]. Several studies have addressed the role of particle deformability [17–20], concentration [21–24] and shape of bifurcation [16, 24] on the fluid-particle partitioning. Finally, the enrichment/depletion of particles in the two outlet streams has been proposed as an effective mechanism to design simple microfluidic devices for particle separation [16], fractionation [16] and filtering [15].

The available theoretical, experimental and numerical studies dealing with suspensions flowing through bifurcations have mainly considered Newtonian liquids as suspending medium. However, an increasing interest in using non-Newtonian fluids in microfluidic devices has been recently observed [25]. Indeed, the elastic properties of these fluids can lead to ‘internal’ forces that push the suspended particles towards specific regions of the channel [26, 27]. By properly exploiting these forces, non-Newtonian fluids can be employed to achieve operations like alignment [28–32] and separation [33, 34] in simple microdevices.

The flow of non-Newtonian fluids (without particles) in T-shaped channels with one inlet and two orthogonal outlets has been the subject of a vast literature [35–39], mainly focused on investigating the stress and velocity fields, and the characteristics of recirculation regions that may appear as inertia, shear-thinning and elasticity are varied. For what concerns the flow behavior of solid viscoelastic suspensions in bifurcating channels, we are aware of one experimental work where the effect of non-Newtonian properties on the separation of particles with different size is investigated [40]. Four different fluids, namely a Newtonian and three non-Newtonian fluids with different degree of shear-thinning and elasticity, have been considered. The experimental measurements show that fluid shear-thinning leads to the highest separation efficiency, which is attributed to the migration mechanism [40]. In the strong shear-thinning fluid, the lateral motion is towards the walls and it is fastest for the largest particles. Hence, small and large particles travel on different streamlines and can be separated. This is not the case when the elastic, constant-viscosity fluid is considered since the migration is towards the channel centerline and does not contribute to the separation.

In this work, we carry out a detailed 2D numerical investigation on the partitioning of particles suspended in Newtonian and non-Newtonian fluids flowing through a microfluidic T-shaped bifurcation. The adopted geometrical configuration is such that the microchannel has one inlet and two orthogonal outlets. Inertia is neglected and the suspension is assumed to be dilute, thus, simulations with a single particle suffice to describe the suspension dynamics. The finite element method (FEM) is used to solve the governing equations and a fictitious domain method (FDM) is employed to handle the particle motion. To improve the accuracy of the FDM around the particle, a grid deformation method (GDM) is also used whereby the mesh is dynamically adapted in order to distribute the elements around the particle boundary. To highlight the effect of the fluid rheological properties, different constitutive equations are considered, namely the Newtonian fluid (constant-viscosity, inelastic), the Bird-Carreau model (shear-thinning, inelastic) and the Giesekus model (shear-thinning, elastic). The partitioning of the particles between the two downstream branches is investigated by varying the inlet flow rate, the relative weight of the two outlet flow rates and the ratio of the particle and channel size.

The paper is organized as follows. In Section “Methods”, the governing equations are presented along with the constitutive models used in this work. The predictions of the constitutive equations in shear and Poiseuille flows are briefly reported. In the same Section, the relevant dimensionless parameters, the adopted numerical method and the details of the computational domain are also discussed. In Section “Results and discussion”, the simulation results are presented in terms of deviation between the particle and suspension fluxes through one outlet, both normalized by the corresponding inlet values. Finally, in Section “Conclusions”, final considerations are drawn.

Methods

Governing equations

Figure 1

Schematic representation of the T-shaped channel. a Global view of the channel. The blue trajectory, at a distance \(y_\text {f}\) from the inlet channel centerline, divides the domain in two regions corresponding to fractions of the fluid without particles that end up in the two outlets. The green trajectory, at a distance \(y_\text {c}\) from the inlet channel centerline, divides the domain in two regions corresponding to fractions of particles that end up in the two outlets. These two limiting trajectories have been drawn assuming that \(Q_1 < Q_2\). b Notation adopted for the domain boundaries. c Close view of the smoothed corners between the main and the side channels.

The geometry considered in this work is a T-shaped channel with one inlet and two outlets, as depicted in Figure 1a. The two channel outlets are orthogonal and all the channels have the same height H. The outlet channel having the same axis of symmetry of the inlet one is denoted as ‘main branch’ or ‘branch 1’. The other outlet, having the axis of symmetry orthogonal to that of the inlet one, is denoted as ‘side branch’ or ‘branch 2’. A suspension of circular particles with diameter \(d_\text {p}=2r_\text {p}\) is injected in the channel with a flow rate \(Q_\text {IN}\). The suspension is assumed to be very dilute which implies that particle-particle hydrodynamic interactions can be neglected. Thus, the dynamics of a single particle suffices to describe the problem. We denote by \(\Omega \) the fluid and particle domain. The external boundaries are denoted as shown in Figure 1b: \(\Gamma _\text {IN}\) is the inflow section, \(\Gamma _\text {OUT1}\) and \(\Gamma _\text {OUT2}\) are the outflow sections of the main and side branches, respectively, \(\Gamma _\text {W}\) collects the remaining boundaries. Notice that, as illustrated in Figure 1c, the corners between the main and side channels are rounded to avoid numerical problems due to point singularity. The space occupied by the particle is denoted by P(t) with boundary \(\partial P(t)\). Notice that, since the particle is transported by the fluid, the particle domain depends on time. A Cartesian reference frame is selected with the origin on the midpoint of the inlet section (see Figure 1a). The x-axis coincides with the direction of the fluid in the inflow channel and the \(y-\)axis identifies the direction of the inflow boundary \(\Gamma _\text {IN}\).

Due to the small length scales adopted in microfluidic devices, the Reynolds number is generally very small [41], thus inertia can be neglected. Assuming also incompressible flow conditions, the equations governing the dynamics of the fluid domain \(\Omega -P(t)\) are:

where \(\varvec{u}\) is the fluid velocity and \(\varvec{\sigma }\) is the stress tensor. These equations are the mass and momentum balance equations, respectively. A general expression for the stress tensor is:

where p is the pressure, \(\varvec{I}\) is the unit tensor, \(\eta \) is the viscosity, \(\varvec{D}=(\nabla \varvec{u}+(\nabla \varvec{u})^{T})/2\) is the rate-of-deformation tensor and \(\varvec{\tau }\) is the extra viscoelastic stress contribution. Notice that the viscosity is, in general, not constant, but it is assumed to be a function of the magnitude of the rate-of-deformation tensor \(\dot{\gamma }=\sqrt{2\varvec{D}:\varvec{D}}\) [42]. By selecting the constitutive equations for the viscosity \(\eta (\dot{\gamma })\) and for the viscoelastic stress tensor \(\varvec{\tau }\), fluids with different rheological properties can be modeled. Details on the constitutive models adopted in this work are given in the next section.

Due to the inertialess assumption and in absence of external forces and torques, the particle is force- and torque-free, i.e. the total force \(\varvec{F}\) and torque \(\varvec{T}\) acting on its boundary are zero:

In these equations, \(\varvec{x}\) is the position vector of a point on the boundary \(\partial P(t)\), \(\varvec{X}\) is the position of the particle center, \(\varvec{n}\) is the outwardly directed unit normal vector on \(\partial P(t)\) and ds is the local boundary length.

To solve the governing equations, boundary conditions need to be specified on the external fluid boundaries and on the fluid-particle interface. No-slip conditions are applied on the channel walls:

The function \(\varvec{u}_\text {IN}(y)\) is the well-known parabolic profile for a Newtonian suspending fluid. However, no analytical expression is available for a general non-Newtonian fluid. In this case, the velocity profile is precalculated by simulating the fluid without particles in a periodic channel with height H so that a fully developed flow is obtained. For a viscoelastic fluid (\(\varvec{\tau }\ne \varvec{0}\)), the viscoelastic stress tensor profile needs to be specified on the inflow section as well:

Like the velocity field, \(\varvec{\tau }_\text {IN}(y)\) is taken from a preliminary simulation of the fluid without particles in a periodic channel. In the viscoelastic case, the steady-state velocity and stress profiles are applied on the inflow section.

where \(\varvec{n}\) is the outwardly directed unit normal vector on the two outflow sections. These conditions assume that the fluid, after crossing the outlet sections of the main and side branches, ends up in an ambient with pressure \(p_\text {OUT1}\) and \(p_\text {OUT2}\), respectively. Notice that the pressure level in a point of the domain does not need to be set because it is specified through the outflow equations. Notice also that the relevant quantity is the pressure difference \(\Delta p=p_\text {OUT1}-p_\text {OUT2}\) rather than the single pressure values. Indeed, the pressure difference \(\Delta p\) sets the relative weight of the two outlet flow rates. For a fixed \(\Delta p\), different values of the outflow pressures only change the level of the fluid pressure field in the channel.

Finally, the rigid-body motion condition is assumed on the particle boundary:

where \(\varvec{U}\) and \(\varvec{\omega }\) are the translational and angular velocities of the particle.

Concerning initial conditions, no initial velocity is needed due to the inertialess assumption. For the viscoelastic fluid, we need to specify the initial stress tensor profile in the fluid domain. We select an initial stress-free state condition in the whole domain:

$$\begin{aligned} \varvec{\tau }|_{t=0}=\varvec{0} \end{aligned}$$

(12)

Once the fluid velocity, pressure and, for the viscoelastic case, the stress fields are calculated along with the particle kinematic quantities, the particle position is updated by integrating the following equation:

Constitutive equations

The stress tensor in Eq. (3) needs to be specified by selecting the constitutive equations for the viscosity \(\eta (\dot{\gamma })\) and for the viscoelastic stress tensor \(\varvec{\tau }\). In this work, we consider three suspending fluids: the Newtonian fluid, a generalized Newtonian fluid modeled by the Bird-Carreau constitutive equation [42, 43] and a viscoelastic fluid modeled by the Giesekus constitutive equation [42, 43]. The equations of these three models, their rheological features in shear flow and the velocity profiles obtained in an infinite slit channel (like the one corresponding to the inlet channel of the geometry investigated in this work) are reported below.

which is characterized by a viscosity that depends on the shear rate and no elasticity. In Eq. (16), the parameter \(n<1\) sets the fluid shear-thinning, \(\eta _0\) and \(\eta _{\infty }\) are the viscosities at low and high shear rates, and 1/K is the shear rate value corresponding to the start of viscosity thinning in simple shear flow. For \(n=1\), a Newtonian fluid with \(\eta _\text {s}=\eta _0\) is recovered.

where \(\eta _\text {s}\) is the (constant) viscosity of a Newtonian solvent, \(\lambda \) is the fluid relaxation time, \(\eta _\text {p}\) is the polymer contribution to the viscosity and \(\alpha \) is the so-called ‘mobility’ parameter. The symbol \((^{\nabla })\) denotes the upper-convected time derivative defined in Eq. (20). In shear flow, the Giesekus constitutive equation predicts a non-zero first normal stress difference \(N_1=\sigma _{xx}-\sigma _{yy}\). For \(\alpha > 0\), the model also predicts shear-thinning and a non-zero second normal stress difference \(N_2=\sigma _{yy}-\sigma _{zz}\). The Newtonian case is recovered for \(\lambda =0\).

Figure 2a shows the viscosity trends predicted by the above constitutive equations in shear flow with the constitutive parameters reported in the top part of Table 1. First of all, notice that all the models have the same zero-shear viscosity \(\eta _0\). For the Giesekus model, \(\alpha \) is set to 0.2 denoting a shear-thinning fluid. The parameters of Bird-Carreau equation are selected in order to match the viscosity trend of the Giesekus fluid. This model is referred as ‘Bird-Carreau A’, characterized by \(n=0.25\). As the Bird-Carreau model does not predict elasticity but only shear-thinning, by matching the viscosity curve of this model with the one obtained from the Giesekus equation, we will be able to separately investigate the effects of shear-thinning and elasticity on the particle partitioning in the T-shaped channel. Furthermore, we will also perform simulations with the Bird-Carreau model for a lower value of the parameter n (without varying the other parameters) to investigate the influence of a more pronounced shear-thinning. This model is referred as ‘Bird-Carreau B’, where \(n=0.05\). Finally, in the inset of Figure 2a, the first normal stress difference coefficient \(\Psi _1=N_1/\dot{\gamma }^2\) is reported for the Giesekus model.

Figure 2

Rheology and velocity predictions of the fluids. a Viscosity trends in shear flow predicted by the constitutive equations used in this work. The model parameters are reported in Table 1. The inset shows the trend of the first normal stress difference coefficient for the Giesekus fluid. b Velocity profiles obtained in a wide-slit channel for the same constitutive equations as in panel (a). The velocity is normalized by its average value whereas the position along the gap is normalized by the channel height.

Table 1

Parameters of the constitutive equations used in this work in dimensional (top) and dimensionless (bottom) form

Newtonian

Bird-Carreau A

Bird-Carreau B

Giesekus

\(\eta _\text {s}\) [Pa]

1.1

–

–

0.1

\(\eta _0\) [Pa]

–

1.1

1.1

–

\(\eta _{\infty }\) [Pa]

–

0.1

0.1

–

n [−]

–

0.25

0.05

–

K [s]

–

1.0

1.0

–

\(\lambda \) [s]

–

–

–

1.0

\(\eta _\text {p}\) [Pa]

–

–

–

1.0

\(\alpha \) [−]

–

–

–

0.2

\(\Lambda \)

–

0.33–1.33

0.33–1.33

0.33–1.33

\(\eta _\text {r}\)

–

0.091

0.091

0.091

n

–

0.25

0.05

–

\(\alpha \)

–

–

–

0.2

In Figure 2b, the velocity profiles obtained in a 2D Poiseuille flow (i.e., in a channel made by two infinite parallel walls) are shown for the aforementioned constitutive equations with the parameters reported in the top part of Table 1. The channel height is set to \(H=1\) mm and the flow rate to \(Q=1\) mm/s, corresponding, for the Newtonian fluid, to a shear rate at wall of \(\dot{\gamma }_\text {w}=6\) s−1. It is readily noted that the shear-thinning predicted by the Giesekus and Bird-Carreau models alters the parabolic velocity profile observed in the constant viscosity fluid. (The shear rate at wall estimated from the Giesekus velocity profile is around 8.5 s−1, so we are well within the shear-thinning region, see Figure 2a). Specifically, a flatten profile is observed around the channel centerline. Notice also that the profile corresponding to the Giesekus fluid (blue line) essentially coincides with that obtained from the Bird-Carreau A model (red solid line) because of the very similar viscosity trend shown in Figure 2a. The dashed red line corresponding to the Bird-Carreau B model is slightly more flatten due to the more pronounced viscosity thinning. Of course, at lower flow rates, the non-Newtonian velocity profiles tend to the Newtonian one as the maximum shear rate value (i.e., at wall) becomes smaller and, hence, the viscosity remains essentially constant along the gap. In contrast, at higher flow rates, more pronounced deviations from the parabolic profile are expected as the viscosity strongly changes through the channel gap due to the fluid shear-thinning. As a final comment, we would like to remark that the velocity profiles in Figure 2b refer to a fluid without particles. The presence of a particle, of course, alters the velocity profile and, in turn, the shear rate distribution. In particular, since the presence of the particle strongly deforms the flow field, the maximum shear rate can be much higher than the shear rate at wall, especially in highly confined situations. Therefore, the local viscosity around a particle suspended in a shear-thinning fluid might be much lower the the zero-shear value.

Dimensionless equations

The governing equations together with the constitutive equations can be made dimensionless by selecting appropriate characteristic quantities for length, velocity and stress. We choose the channel height H as characteristic length, the average velocity in the inflow channel \(\bar{u}_\text {IN}=Q_\text {IN}/H\) as characteristic velocity and \(\sigma _\text {c}=\eta _0 \bar{u}_\text {IN}/H\) as characteristic stress. By using these characteristic scales, a number of dimensionless parameters appear that are listed in the bottom part of Table 1 for the different constitutive equations.

First of all, we notice that no dimensionless quantity appears for the case of a Newtonian suspending fluid. Indeed, the fluid viscosity \(\eta _\text {s}\), that is the only dimensional parameter of this constitutive equation, is used in the definition of the characteristic stress and disappears from the dimensionless form of the equations. This implies that the value of the flow rate \(Q_\text {IN}\) has no influence on the particle dynamics in dimensionless terms. Of course, this is valid as long as inertial terms are neglected.

For the Bird-Carreau and Giesekus models, together with the parameters n and \(\alpha \) that are already dimensionless in Eqs. (17) and (19), other two parameters appear: the viscosity ratio \(\eta _\text {r}\) and the parameter \(\Lambda \). The viscosity ratio is defined as \(\eta _\text {r}=\eta _\infty /\eta _0\) and \(\eta _\text {r}=\eta _\text {s}/\eta _0\) for the BC and GSK models, respectively. The parameter \(\Lambda \) is defined as \(\Lambda =K Q_\text {IN}/H^2\) and \(\Lambda =\lambda Q_\text {IN}/H^2\) for the BC and GSK models, respectively. Notice that, both K and \(\lambda \) have the dimension of time. However, whereas \(\lambda \) is the fluid characteristic time related to the stress build-up, K has nothing to do with memory effects. Since \(\lambda \), K and H are kept constant throughout this work, the parameter \(\Lambda \) is changed by varying the inlet flow rate \(Q_\text {IN}\). Notice that the parmeter \(\Lambda \) for the viscoelastic fluid coincides with the Deborah number [25].

By making the boundary conditions dimensionless, the dimensionless relative pressure difference \(\Delta p^{*}=(p_\text {OUT1}-p_\text {OUT2})/\sigma _\text {c}\) appears. Finally, the last dimensionless parameter to be specified is the confinement ratio \(\beta =d_\text {p}/H\), i.e. the ratio between the particle diameter and the channel height.

Numerical method

The governing equations are solved by the finite element method (FEM). To handle the particle motion, a fictitious domain method (FDM) is used [44]. In this method, the particles are filled by the same fluid of the external phase. The fluid and solid domains are embedded into a single weak form. The balance of forces and torques at fluid-particle boundaries is satisfied, but the forces and torques do not explicitly end up in the variational form. To release the constraint on the variational space, the rigid-body motion on the particle surface is imposed through Lagrange multipliers. In this work, we adopt a weak implementation of the constraints that has been proven to give higher accuracies and better conditioned linear systems as compared to the collocation method [45]. Therefore, the kinematic quantities of the particle are unknowns of the full system of equations and are directly computed by solving the system. The main advantage of the fictitious domain method is the possibility to use a fixed, time-independent, regular mesh for both fluid and solid. Concerning the particle domain discretization, due to the inertialess assumption, we can adopt a rigid-ring description of the particle [46, 47], i.e. the particle is described as a rigid-shell and the rigid-body motion is imposed only on the particle boundary.

For the viscoelastic fluid case, a two-step decoupled formulation is adopted [48], whereby the momentum and continuity equations are decoupled from the constitutive equation. At each time level, the stress field is first updated by solving Eq. (19) where the velocity field is taken from the previous step. Then, the viscoelastic stress tensor just computed is used as known force term in the momentum balance which is solved together with the continuity equation and the force- and torque-free conditions to get the velocity and pressure fields, and the particle kinematic quantities. Finally, the particle position is updated by integrating Eq. (13). A DEVSS-G/SUPG formulation is adopted [49–51] with a log-conformation representation [52, 53] to improve the numerical stability.

For the discretization of the weak form, we use quadrilateral elements with continuous biquadratic interpolation (Q2) for the velocity and bilinear continuous interpolation (Q1) for the pressure. In the viscoelastic case, bilinear continuous interpolation (Q1) for the velocity gradient (coming from the DEVSS-G scheme) and bilinear continuous interpolation (Q1) for the log-conformation tensor are used. More details on the weak form and the implementation are given in Refs. [48, 54].

Figure 3

Grid deformation method. Close view of the mesh around a circular particle near the right corner of the T-shaped channel without (a) and with (b) the application of the grid deformation algorithm.

One major drawback of the fictitious domain method is a relatively low accuracy around the particle due to the fact that the field variables at the particle boundary are interpolated between the values possessed by the physical fluid (outside the particle) and the fictitious fluid (inside the particle). Consequently, discontinuities required at the interface are lost. To reduce such a problem, we adopt a grid deformation algorithm (GDA) to locally refine the mesh around the particle boundary [55]. The method redistributes the mesh nodes while preserving the mesh topology. An example of the mesh around a particle near the right corner of the T-channel is reported in Figure 3 without (left) and with (right) applying the grid deformation method. The details of the algorithm are given in Ref. [55]. Therefore, by combining the FDM and the GDA, the mesh is not fixed and time-independent anymore. However, its topology is preserved and, at each time step, it can be readily obtained for a given particle position from an initial regular mesh. As a final remark, since the mesh moves each time step, an issue arises when solving the constitutive equation for the viscoelastic case. Indeed, the time integration of Eq. (19) requires the knowledge of the stress tensor field at previous time levels corresponding to the grid node positions of the current time level. We solve this problem by adopting an Arbitrary Lagrangian-Eulerian (ALE) procedure [56] whereby the convective term \(\varvec{u}\cdot \nabla \varvec{\tau }\) in Eq. (20) is replaced by \((\varvec{u}-\varvec{u}_\text {m})\cdot \nabla \varvec{\tau }\) where \(\varvec{u}_\text {m}\) is the node mesh velocity. The grid velocity is calculated as: \(\varvec{u}_\text {m}=(\varvec{x}_\text {m}^{n+1}-\varvec{x}_\text {m}^{n})/\Delta t\) where \(\varvec{x}_\text {m}^{n+1}\) and \(\varvec{x}_\text {m}^{n}\) are the positions of the mesh nodes at the current and previous time level, and \(\Delta t\) is the time step size.

Selection of the computational domain

The selection of the domain geometrical parameters is a critical issue. Indeed, the domain size should be chosen as small as possible to reduce the computational cost. On the other hand, the simulated geometry must be representative of a real device where the three branches are generally much longer than the channel height, resulting in a fully developed flow far from the bifurcation. In this section, we discuss the procedure adopted to select the geometrical parameters of the computational domain.

As mentioned above, we impose the velocity (and the stress for the viscoelastic case) profile on the inflow section \(\Gamma _\text {IN}\) and outflow conditions on the outlet sections \(\Gamma _\text {OUT1}\) and \(\Gamma _\text {OUT2}\). To avoid that the imposed boundary conditions affect the particle motion near the bifurcation, the lengths of the three branches must be chosen sufficiently long. Furthermore, the particle is initially placed at a distance \(X_0\) from the inflow section (see Figure 1a). Also this distance needs to be selected sufficiently large such that the inflow boundary condition does not influence the particle motion. Preliminary simulations have been performed by changing the branch lengths as well as the distance \(X_0\). For the Newtonian and Bird-Carreau models, we found that a good choice of these geometrical quantities is \(L_\text {IN}/H=L_1/H=L_2/H=3\) and \(X_0/H=0.5\). Indeed, for these values, we observe that, as soon as the simulation starts, the particle follows a straight trajectory within the inlet channel at a constant translational velocity (recall that, under inertialess conditions, there is no time-derivative and the velocity and pressure fields instantaneously develop), i.e., its initial dynamics is unaffected by both the inflow profile and the downstream bifurcation. As the particle approaches the T-junction, the path deviates and the particle enters one of the two branches. Finally, when sufficiently far from the bifurcation, the particle travels again along a straight streamline at a constant velocity, i.e. its dynamics is unaffected by the outflow conditions. The simulation is, then, stopped before reaching the outlet section.

Concerning the viscoelastic case, a time-derivative appears in the governing equation through Eq. (20). Therefore, the stress field develops from the stress-free initial condition given in Eq. (12) to a ‘steady-state’ (actually, it is not a true steady-state since the particle moves and, in turn, the surrounding stress field changes in time). Of course, the velocity field and the kinematic particle quantities undergo the same transient dynamics. However, as discussed above, our simulations must be representative of fully developed flow conditions. To satisfy this requirement, we adopt the following procedure. We select a value \(L_\text {IN}\) for the inlet channel length and we place the particle at a distance \(X_0\) from the inflow section. We start the simulation and check the evolution of the particle position and its translational velocity. The velocity trend shows an initial start-up corresponding to the viscoelastic stress development. The choice made for \(L_\text {IN}\) is correct whether, after the initial transient, (1) the translational velocity becomes constant and (2) the particle travels within the inlet channel following a straight path. If both conditions are satisfied, the situation just after the start-up is equivalent to that obtained by running the same simulation in a much longer inflow channel. Hence, we can conclude that \(L_\text {IN}\) has been chosen sufficiently large to assure a fully developed flow far from the bifurcation. If this is not the case, the inflow length \(L_\text {IN}\) is increased and the simulation is run again. In this respect, it should be mentioned that, as the stresses develop, the particle may change its y-position due to the migration mechanism [57]. However, we found that the lateral particle displacement is negligible in all the simulations carried out in this work, due to the short duration of the start-up phase. This is also confirmed by previous simulations in similar conditions [57].

Of course, the space covered by the particle during the start-up increases as the flow rate is higher. As a consequence, \(L_\text {IN}\) depends on the inlet flow rate. An inlet channel length of \(L_\text {IN}/H=9\) is required for the two lowest flow rates considered in this work. The length must be increased to \(L_\text {IN}/H=12\) for the highest flow rate. In Table 2, the values of the geometrical parameters are summarized.

Concerning the adopted mesh, the reference (undeformed) grid is refined along the walls and, in particular, around the corners where the largest field gradients are expected. Spatial and temporal convergence is verified for all the simulations performed in this work. As expected, the most critical case is the viscoelastic one at high \(\Lambda \)-values. For the Newtonian and Bird-Carreau models, about 20,000 quadrilateral elements are sufficient to get grid independence of the solution. The elements must be increased to about 35,000 when the Giesekus model is considered. A time step size of \(\Delta t=0.01\) guarantees time convergence for all the simulations. A typical simulation for a Newtonian fluid takes about 3-4 h of computational time. The CPU time increases by about a factor 3 for the Bird-Carreau case because of the Newton-Rhapson iterations required at each time step. Finally, one day is needed to get the results for the viscoelastic suspending fluid due to the finer mesh and to the increased size of the system of equations as a consequence of the adopted DEVSS-G procedure.

As a final note, we recall that all the simulations are performed by setting the inlet flow rate (through the velocity profile) and the ambient pressures \(p_\text {OUT1}\) and \(p_\text {OUT2}\) on the two outflow sections. Of course, the flow rate partitioning through the outlet branches depends on the pressure difference \(\Delta p=p_\text {OUT1}-p_\text {OUT2}\) as well as on the lengths of the outlet branches \(L_1\) and \(L_2\). It is convenient to present the results in terms of outlet flow rates rather than outlet pressure difference. In this way, in fact, the results are general and independent of the lengths of the outflow channels.

Table 2

Geometrical parameters used in this work. See Figure 1 for the notation

Newtonian/Bird-Carreau

Giesekus

\(L_\text {IN}/H\)

3

9–12

\(L_\text {1}/H\)

3

3

\(L_\text {2}/H\)

3

3

R / H

0.2

0.2

\(X_{0}/H\)

0.5

0.5

Results and discussion

Let us consider a fluid without particles that is pumped through the inlet branch of a T-shaped channel. Near the bifurcation, the fluid is partitioned between the two outlets according to the imposed pressure difference \(\Delta p=p_\text {OUT1}-p_\text {OUT2}\). Specifically, if \(\Delta p=0\) and if the outlet channels have the same height and length, it is \(Q_1 \cong Q_2\), where \(Q_1\) and \(Q_2\) are the flow rates through the main and side outlet branches, respectively. (Actually, these two flow rates are not exactly equal because the bifurcation is asymmetric). In case of \(\Delta p>0\) (\(\Delta p<0\)), the fluid encounters a higher resistance in the main channel (side channel) and prevalently moves through the other branch, i.e., \(Q_2 > Q_1\) (\(Q_1 > Q_2\)). Of course, there will be two limiting values of \(\Delta p\) (one positive and one negative) such that the fluid totally passes through only one outlet. In the general case, a critical streamline dividing the inflow section in two regions corresponding to the fluid portion that ends up in the two outlets can be identified. An example of this streamline is depicted as a blue line in Figure 1a, which is at a coordinate \(y_\text {f}\) along the section of the inlet branch.

Let us now consider a particle suspended in a fluid and initially moving in the inlet branch of the T-shaped channel far from the junction. The particle is transported by the fluid and deviates its path while approaching the bifurcation. Depending on the relative pressure between the two outlets and on the position of the particle center along the section of the inlet branch, the particle can be displaced to the side branch or follow the main flow direction and exit the channel from the main branch. Similarly to the above mentioned critical streamline, for a fixed pressure difference \(\Delta p\), a critical trajectory can be identified that divides the inlet channel section in two parts corresponding to paths that lead the particle to the main (above the critical trajectory) or to the side (below the critical trajectory) outlets. An example of the critical trajectory is reported as green line in Figure 1a. Of course, the critical trajectory changes by varying \(\Delta p\). The key point is that, in general, the critical streamline and the critical trajectory do not coincide. Indeed, the hydrodynamic interaction between the particle and the right corner of the bifurcation alters the path followed by the finite-sized particle as compared to the fluid particle (that has no size). This hydrodynamic interaction is more pronounced for bigger particles and, in fact, larger deviations between the two critical separatrices are expected in this case. If the two separatrices are different, the fraction of fluid flowing through the two outflow branches is, hence, different from the fraction of particles collected from the two outlets. Therefore, the two outlet streams have a different relative amount of suspended particles as compared to the inlet one. In the case of a bi-disperse suspension, this method could be used to fractionate the sample suspension and enrich the outlet streams with bigger or smaller particles.

To quantify the fractionation of particles between the two outlets, we introduce the particle flux ratio through the main outlet [13]:

where \(F_1\) and \(F_\text {IN}\) are the fluxes of particles through the main outlet and the inlet, respectively, \(U_\text {x}(y)\) is the particle translational velocity along the flow direction in the inflow channel far from the bifurcation, and \(y_\text {c}\) is the position of the critical trajectory along the section of the channel inlet (see Figure 1). First of all, we remark that both integrals are computed along one section of the inlet branch corresponding to fully developed flow conditions. Notice that the definition in Eq. (21) assumes a uniform particle distribution in inflow. Comments about such an assumption will be given at the end of this section. Finally, notice also that the integrations are performed over the channel section accessible to the particle (i.e., the zone near the walls with thickness of one radius is excluded). The values \(F_1/F_\text {IN}=0\) or \(F_1/F_\text {IN}=1\) correspond to the limiting conditions such that all the particles end up in the side or the main branches, respectively.

In this work, we carry out simulations by varying \(\Lambda \) (through the inlet flow rate), the relative weight of the outlet flow rates (through the pressure difference \(\Delta p^{*}\)) and the confinement ratio \(\beta \), for the three constitutive equations discussed in the previous section. For each set of these parameters, we evaluate the particle flux ratio from Eq. (21). This requires the knowledge of the particle velocity profile \(U_\text {x}(y)\) along the channel inlet as well as the position \(y_\text {c}\) of the critical trajectory. The computations of these two quantities is done by adopting the following procedure. For a fixed set of parameters, we release the particle at 10 equally distributed positions along the cross-section of the inflow channel, at a distance \(X_0\) from the inflow boundary, selected as discussed in the previous section. (Due to numerical issues, the closest distance between the particle center and the walls is set to 1.1\(r_\text {p}\)). For each initial position, we compute the x-component of the translational velocity \(U_\text {x}\). For the inelastic constitutive equations, a single time step is sufficient to compute the translational velocity since inertia is neglected. For the viscoelastic case, \(U_\text {x}\) is computed after the stress build-up. The simulation data are, then, interpolated and an analytical velocity profile \(U_\text {x}(y)\) is reconstructed. To compute \(y_\text {c}\), we run simulations by initially locating the particle at different y-positions on the section of the inlet branch. We first detect the two closest trajectories that lead the particle to the main and side outlets. Then, we run other simulations by initially locating the particle between the two initial positions corresponding to these trajectories in order to identify the critical path with higher accuracy. The uncertainty in computing \(y_\text {c}\) is set to \(\pm 0.5\%\) of the channel height H, i.e. a particle released at \(y_\text {c}/H+0.005\) ends up in the main outlet and a particle released at \(y_\text {c}/H-0.005\) ends up in the side outlet.

Figure 4

Particle partitioning in a Newtonian fluid. a Particle flux ratio as a function of the flow rate ratio for three different confinement ratios. The inset shows experimental data taken from [16]. b The same data of (a) plotted in terms of difference of particle flux ratio minus flow rate ratio.

We start by presenting the results for particles suspended in the Newtonian fluid. Figure 4a reports the particle flux ratio \(F_1/F_\text {IN}\) as a function of the flow rate ratio \(Q_1/Q_\text {IN}\) for three different confinement ratios \(\beta =0.1\) (red circles), \(\beta =0.2\) (green circles) and \(\beta =0.5\) (blue circles). A point of the diagonal in this diagram (solid black line) indicates that the same fraction of particles is obtained between the inflow and outflow streams, i.e. the inflow particles are fractioned through the two outlets proportionally to the fluid flow rates. On the contrary, a point above (below) the diagonal denotes that the fluid exiting the main outlet (the side outlet) is enriched of particles as compared to the inlet. The enrichment is more and more pronounced as the deviation from the diagonal is large. The data reported in Figure 4a show that fractionation is enhanced for a flow rate ratio above 0.5 and for the largest particle size considered. For lower \(Q_1/Q_\text {IN}\) values, only slight (negative) deviations are noted even for the largest confinement ratio. As the particle becomes smaller, the separation process is less and less efficient regardless of the value of the flow rate ratio (the red and green circles are hardly distinguishable from the diagonal).

Although a quantitative comparison with the available experimental data [15, 16, 21] is not possible due to the 2D assumption of the present simulations, the results in Figure 4 fairly reproduce the experimental trends. To show this, the inset of Figure 4a reports the experimental data taken from Ref. [16]. The two sets refer to confinement ratios of \(\beta =0.5\) (blue triangles) and \(\beta =0.77\) (magenta triangles). In both cases, the volume fraction is \(\phi =0.02\) so that the system can be assumed dilute. In qualitative agreement with our results, the experimental measurements show that larger particles are more easily separated and the process is more efficient for flow rate ratios above 0.5. Notice also that, in spite of the 2D nature of our simulations, the experimental and numerical data at the same \(\beta \) are quantitatively close.

The way to visualize the simulation results shown in Figure 4a has been used in previous works dealing with particle separation processes in branched capillaries [13, 15, 16, 21]. However, to better quantify the fractionation process, in what follows, we prefer to report the data as in Figure 4b. Here, the quantity \(F_1/F_\text {IN}-Q_1/Q_\text {IN}\) is plotted as a function of the flow rate ratio. Therefore, each symbol directly represents the distance from the diagonal that is now the horizontal line at \(y=0\). The data confirm the asymmetric trend around \(Q_1/Q_\text {IN}=0.5\) and the decreasing efficiency of the fractionation process for smaller particle sizes.

Figure 5

Particle partitioning in a Bird-Carreau A fluid \((n=0.25)\). Difference of particle flux ratio minus flow rate ratio as a function of the flow rate ratio for three different confinement ratios: \(\beta =0.1\) (a), \(\beta =0.2\) (b) and \(\beta =0.5\) (c). In each panel, three different \(\Lambda -\)values have been considered. The lines are fits of the simulations data that are added as a guide for eye.

We now investigate the effect of fluid shear-thinning on the particle partitioning. To this aim, we first consider the Bird-Carreau A model where the shear-thinning parameter is set to \(n=0.25\). We recall that such a model does not predict elastic effects but only a viscosity dependent on the shear rate. Figure 5 reports the quantity \(F_1/F_\text {IN}-Q_1/Q_\text {IN}\) as a function of the flow rate ratio \(Q_1/Q_\text {IN}\). The three plots correspond to different confinement ratios: \(\beta =0.1\) (a), \(\beta =0.2\) (b) and \(\beta =0.5\) (c). For each confinement ratio, three different values of \(\Lambda \) have been considered: \(\Lambda =0.33\) (close circles), \(\Lambda =0.66\) (open circles) and \(\Lambda =1.33\) (crossed circles). The lines are fits of the simulation data that are added as a guide for the eye. The trends are qualitatively similar to the Newtonian case. The data, indeed, are negative (positive) for flow rate ratios lower (higher) than a value slightly below 0.5. Furthermore, larger confinement ratios improve the fractionation as clearly visible by the scales of the \(y-\)axis of the three panels. At variance with the Newtonian case, the inlet flow rate \(Q_\text {IN}\) (contained in \(\Lambda \)) now affects the shear rate and the velocity profile along the channel cross-section (see Figure 2b). At the lowest value of \(\Lambda \), the shear rate at wall is around 2.5, corresponding to a moderate viscosity thinning. Nevertheless, by comparing the close circles in Figure 5 with the data for the Newtonian case in Figure 4b, very similar quantitative trends are found for all the confinement ratios considered. An improvement of the fractionation process with respect to the one found in a Newtonian suspending fluid is only observed for the highest confinement ratio at \(Q_1/Q_\text {IN}< 0.5\) (close circles in Figure 5c). In this region, the minimum value of \(F_1/F_\text {IN}-Q_1/Q_\text {IN}\) in the Bird-Carreau A fluid is twice lower than the Newtonian case, and comparable with the maximum value observed for \(Q_1/Q_\text {IN}> 0.5\) in the same fluid. It so appears that shear-thinning does not seem to significantly affect the distribution of particles in a T-shaped bifurcation. Indeed, by increasing \(\Lambda \), i.e. by increasing the inlet flow rate \(Q_\text {IN}\), corresponding to a more pronounced viscosity thinning, the data follow the same qualitative and quantitative trend of the minimum flow rate case just discussed. Specifically, for \(\beta =0.1\) and \(\beta =0.2\), the three sets of data describe a unique trend. For \(\beta =0.5\), the simulation data for the largest inlet flow rate (crossed circles) are slightly higher than the other two flow rates. A possible explanation might be related to the high local shear rate values around a highly confined particle (see the end of section ‘Constitutive equations’) that contribute to reduce the local viscosity. However, the observed deviations from the lower flow rate cases are quite small, denoting a weak effect of shear-thinning on the particle fractionation.

Figure 6

Particle partitioning in a Bird-Carreau B fluid \((n=0.05)\). Difference of particle flux ratio minus flow rate ratio as a function of the flow rate ratio for three different confinement ratios: \(\beta =0.1\) (a), \(\beta =0.2\) (b) and \(\beta =0.5\) (c). In each panel, three different \(\Lambda -\)values have been considered. The lines are fits of the simulations data that are added as a guide for eye.

Figure 7

Particle partitioning in a Giesekus fluid. Difference of particle flux ratio minus flow rate ratio as a function of the flow rate ratio for three different confinement ratios: \(\beta =0.1\) (a), \(\beta =0.2\) (b) and \(\beta =0.5\) (c). In each panel, three different \(\Lambda -\)values have been considered. The lines are fits of the simulations data that are added as a guide for eye. In (a) and (b), the solid lines have been obtained by considering the data at \(Q_\text {IN}=1/3\) and \(Q_\text {IN}=2/3\), whereas the dashed lines refer to the data at \(Q_\text {IN}=4/3\).

To further investigate on the effect of shear-thinning, we repeated the same simulations by decreasing the parameter n to 0.05 (Bird-Carreau B, red dashed curves in Figure 2). The simulation results are shown in Figure 6. The same notation of Figure 5 is used. First of all, we notice that the data for each set of parameters are a bit more scattered than the previous case. This might be ascribed to numerical convergence issues arising when dealing with low n-values. A comparison between the data in Figure 6 and Figure 5 shows that roughly the same particle paritioning is obtained for both \(\beta =0.1\) and \(\beta =0.2\) at any investigated value of the parameter \(\Lambda \). A slightly improved fractionation is, instead, observed for the largest particle size and for the most shear-thinning fluid (compare the scales of the y-axis between Figures 6c and 5c). As mentioned above, the high local shear rate values around the particle that, for \(n=0.05\), lead to even more pronounced viscosity thinning, might be responsible for the enhanced fractionation process. However, as the deviations from the simulation data between the case \(n=0.25\) and \(n=0.05\) are quite small, we can conclude that the viscosity thinning alone does not produce any significant effect on the particle distribution downstream a T-junction as compared to a Newtonian suspending fluid.

Figure 8

Particle partitioning in the Newtonian, Bird-Carreau A and Giesekus fluids. Difference of particle flux ratio minus flow rate ratio as a function of the flow rate ratio for three different confinement ratios: \(\beta =0.1\) (a), \(\beta =0.2\) (b) and \(\beta =0.5\) (c). In each panel, the simulation data for three different suspending fluids are reported with \(\Lambda =1.33\). The lines are fits of the simulations data that are added as a guide for eye. The solid lines have been obtained by considering the data for the Newtonian and the Bird-Carreau A models, whereas the dashed lines refer to the data for Giesekus model.

As final case, we discuss now the partitioning of particles in the T-shaped channel when a viscoelastic fluid (Giesekus model) is used as a suspending medium. The simulation results are shown in Figure 7. First of all, for \(\beta =0.1\) and \(\beta =0.2\), a more efficient separation is found for the highest inlet flow rate (\(\Lambda =1.33\), crossed circles) as compared to the two lower values (close and open circles). (For this reason, we fit separately the two trends that are identified by the solid and dashed curves). Notice that this behavior is at variance with the Bird-Carreau fluid where no significant effect of the parameter \(\Lambda \) was found. By increasing the confinement ratio, the deviations between the data at low and high \(\Lambda \)-values are less and less pronounced and, in fact, the three sets of data follow the same quantitive trend in panel (c) (corresponding to \(\beta =0.5\)). The just discussed results suggest that fluid elasticity, that is enhanced at high flow rates, plays a relevant role on the particle distribution downstream the bifurcation. This is further confirmed by comparing the scales of the y-axis between Figure 7 and the previous figures. The fractionation process obtained in the viscoelastic fluid is from 5 (for \(\beta =0.1\)) to 3 (for \(\beta =0.5\)) times more efficient than the inelastic cases.

A possible explanation of the just mentioned elasticity effect might be related to the normal stress distribution around the particle. It is likely to assume that, as the particle approaches the right corner of the bifurcation, the shear rate gradients in the fluid gap between the particle surface and the wall are very large. The large shear rate gradients lead to large normal stress gradients that deviate the particle from the streamline followed if it is suspended in an inelastic fluid. The shear rate gradients are enhanced whether both the confinement ratio and the flow rate increase. The data reported in Figure 7 show that, for strongly confined conditions, even relatively low flow rates are sufficient to displace the particles (i.e., the large particle size suffices to produce high normal stress gradients even at low flow rates). In contrast, at low confinement ratios, the small particle size is unable to produce relevant shear rate/normal stress gradients and high flow rates are needed to increase the efficiency of the fractionation process.

To summarize the simulation results, we report in Figure 8 the data shown in the previous figures for the three suspending fluids. The three panels differ for the confinement ratio and the data refer to the highest value of the parameter \(\Lambda =1.33\). The trends confirm that the particle partitioning is quantitatively similar for the Newtonian and the shear-thinning inelastic fluids. In contrast, the outlet streams are much more enriched/depleted of particles whether a viscoelastic liquid is used as suspending medium. In particular, the maximum fractionation is observed for \(Q_1/Q_\text {IN}\) between 0.7 and 0.8, for all the three values of the confinement ratio considered in this work.

These results can be used to predict the partitioning of dilute bidisperse suspensions as well, i.e. when particles with different size are simultaneously suspended in the inlet stream. Indeed, the diluteness assumption allows to neglect interparticle hydrodynamic interactions and the partitioning after the T-junction can be still deduced from the single-particle problem. As discussed above, the fractionation is quantitatively different for small or big particles (see the scales of the \(y-\)axis in the panels of Figure 8). Therefore, for a fixed value of \(Q_1/Q_\text {IN}\), the small and big particles will be differently distributed between the two outlet streams (fractionation process). A viscoelastic fluid is preferable as it guarantees a fractionation efficiency in a single device much higher than that obtained in an inelastic fluid. Of course, by changing the relative weight of the outlet flow rates, the partitioning of the particles between the outlet streams can be quantitatively tuned.

As a final comment, it is worthwhile to mention that the definition of the particle flux ratio given in Eq. (21) assumes a uniform particle distribution along the section of the channel inlet. Under inertialess conditions, a single particle suspended in an inelastic fluid flowing in a channel follows the fluid streamlines. Therefore, in case of the particles are randomly distributed when the suspension is injected into the device, the particle distribution remains uniform along the inlet channel (the particle paths are straight lines). This might not be the case when the particles are suspended in a viscoelastic liquid. Indeed, the well-known migration phenomenon, i.e. a motion transversal to the flow direction, may take place [27, 57, 58]. Such a phenomenon is induced by fluid elasticity and, of course, alters the particle distribution along the channel axis. Specifically, it has been shown that the migration is towards the channel axis and, in shear-thinning (and elastic) fluids, also towards the wall [57]. Although the migration velocity depends on several factors (e.g., fluid properties, flow conditions, geometrical parameters), it is generally 2–3 order of magnitude lower than the main flow velocity and a channel much longer than the cross-section characteristic size is required to drive a particle from the wall to the centerline [57]. Therefore, the results presented in this work assume that the length of the inlet channel is relatively short (at most one order of magnitude larger than H) so that migration is negligible. For this reason, a direct comparison with the experimental data reported in Ref. [40] is not possible due to the relevant migration phenomenon that occurs in the experiments.

Conclusions

In this work, we investigated the effect of the fluid rheology on the partitioning of a dilute suspension of particles flowing through a T-shaped bifurcation by 2D direct numerical simulations. A fictitious domain method together with a grid deformation procedure is used to handle the particle motion. Three constitutive equations have been considered: the Newtonian equation (constant-viscosity, inelastic), the Bird-Carreau model (shear-thinning, inelastic) and the Giesekus model (shear-thinning, elastic). Simulations are carried out by changing the inlet flow rate, the relative weight of the two outlet flow rates and the confinement ratio.

The results show that, in agreement with the previous literature, the partitioning of particles at the bifurcation differs from the partitioning of the suspending fluid, regardless of the fluid rheology. The only effect of shear-thinning does not produce any relevant quantitative deviation from the Newtonian suspending fluid case. On the other hand, fluid elasticity greatly enhances the efficiency of the fractionation process. Specifically, the most pronounced deviations between the particle and fluid fractions that end up in the two outlet branches are found for the largest confinement ratio and inlet flow rate.

The results reported in this work show the advantage of using viscoelastic fluids to efficiently separate and fractionate particles by sizes in simple microfluidic devices.

Declarations

Author’s contributions

GDA performed all the computations and drafted the manuscript. MAH and PLM carried out a detailed revision. All authors read and approved the final manuscript.

Compliance with ethical guidelines

Competing interests The authors declare that they have no competing interests.

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.