Particle–liquid transport in curved microchannels: Effect of particle volume fraction and size in Dean flow

Over the past couple of decades, microfluidics has emerged as a popular tool for a multitude of engineering applications and processes (e.g., particle separation and sorting,^1–3 mixing,^4,5 electrochemical energy systems,⁶ sorting of cells,^7,8 cell analysis,^9,10 and biosensing^11,12) owing to factors such as the use of small sample volume, high throughput, and reduced cost. Several of these microfluidic processes intricately depend on the fluid mechanics within the microchannel to impart a desirable motion to the transported moieties (e.g., particles and cells). However, the intrinsically small Reynolds number (Re) associated with the microfluidic transport hinders the ability of the system to produce the desirable flow characteristics necessary to impart the appropriate motion to the transported species (e.g., for low Re flows, it becomes challenging to generate eddies that enable mixing of the transported species). In these micro-scale environments, therefore, researchers have adopted a variety of different approaches to either augment the flow strength or enhance the development of secondary flows with the aim of achieving certain specific applications.

It has been well known that curvature of channels or pipes generates a specific type of secondary flow, known as the Dean flow. Dean flow is attributed to the asymmetry induced within the flow field on account of the curvature-driven centrifugal force experienced by the flow.^13–23 Curvature-induced Dean flows, typically characterized by the presence of two counter-rotating vortices, have been extensively used in curved microchannels for different applications such as cell sorting and isolation^24–29 and inertial particle focusing.^30–34 The extent to which Dean flow impacts the flow path of cells and particles in a curved microchannel stem from the intricate balance of the inertial lift force (e.g., the lift force resulting from the spatial gradient in the flow profile), wall-induced lift force, and the overall drag force experienced by the cells/particles, which in turn determine the focusing position of these cells/particles within the microchannel.

Spiral microchannels generate Dean flow, and there have been extensive studies describing Dean flow in terms of the two induced counter-rotating vortices. These vortices are primarily present for low Re cases (Re < 20),^35–39 whereas vortices are hypothesized to be less prevalent for larger Re cases (Re > 100).⁴⁰ Given the significant implications of using Dean flow effects in microfluidic manipulation, separation, or transport of particles or cells, it is important to understand the possible coupling between the particle and fluid in a spiral microchannel (i.e., microchannel with a curvature) for finitely large volume fraction of the particles. Interestingly, to the best of our knowledge, all the existing numerical and theoretical investigations of microfluidic Dean flow have only quantified the fluid flow (and flow-driven migration of these particles/cells) with no attempt to couple the two motions of fluid flow and the particle/cell transport. While there have been several studies describing the flow-driven particle motion in a framework that appropriately couples these motions,^41–44 there are no such coupled descriptions involving Dean flow in microfluidic curved channels.

In this paper, we describe a particle–liquid two-phase computational fluid dynamics (CFD) model to solve for the coupled fluid-flow-particle-transport in a spiral microchannel giving rise to microfluidic Dean flow. To the best of our knowledge, this is the first computational study that explicitly captures the Dean flow-driven coupled fluid-flow-particle-transport in a spiral microchannel. Our findings explicitly reveal the effect of the volume fraction of the particles as well as the size of the particles in affecting the Dean flow. We start with the analysis of a single-phase fluid flow in a spiral microchannel and identify the distinct Dean flow effect that ensures that the velocity maximum, owing to the presence of the large curvature-induced centrifugal force, gets shifted toward the outer (concave) wall of the microchannel. We validate our CFD simulations by comparing our results against the simulation results of Nivedita et al.¹⁴ With this background, we subsequently probe the transport of the two-phase particle–fluid transport in spiral microchannels. This CFD framework is separately validated against experimental results probing the particle transport in curved channels.⁴⁵ Our simulation results elucidate that the presence of the particles invariably reduces the Dean flow effect and diminishes the tendency of the fluid flow to get skewed toward the outer wall. We further observe that the secondary vortices, which appear in addition to the two counter-rotating vortices for the case of a single-phase Dean flow in the presence of a larger Dean number, disappear for the case of the two-phase flow. We argue that both these scenarios are triggered by the particles exerting distinct drag and lift forces on the fluid flow (these forces are reactions to the drag and the lift forces that the continuous phase exerts on the particles or the dispersed phase), and these forces oppose the large curvature-mediated centrifugal force exerted on the continuous phase. We consider separate cases where we (1) vary the volume fraction of the particles (the discrete phase) for a given size of spherical particles and (2) vary the diameter of spherical particles for a given volume fraction of the particles. We separately study the variation of the drag and lift forces (exerted by the particles on the flow) as functions of these parameters as well as the Dean number and relate the manner in which the velocity profiles affect these forces. Of course, these forces themselves dictate the velocity profile, given the coupled nature of the problem. From these analyses, we are able to better understand the physics of the alteration of the flow and vortex profiles as a function of the presence of the particles under different particle and flow conditions, which results in different Dean numbers.

In this study, we consider the flow of a Newtonian fluid seeded with spherical particles in a spiral microchannel, as shown in Fig. 1. This system is modeled using the mixture theory as a two-phase interacting continuum in an Euler–Euler framework.^42,46–48 The two phases are the fluid phase (or the continuous phase) and the particle phase (or the discrete phase). We investigate the effect of the Reynolds number (or the equivalent Dean number), the initial volume fraction of the loaded particles, and the fluid–solid interacting forces on the fluid flow inside the spiral microchannel. Here, $Re=ρcUcDhμc0$ is the Reynolds number, $De=ReDh2R$ is the Dean number, $ρc$ is the pure fluid density (or the density of the continuous phase), U_c is the average fluid velocity, D_h is the hydraulic diameter of the microchannel, $μc0$ is the pure fluid viscosity (or the dynamic viscosity of the continuous phase), and R is the radius of curvature of the spiral channel.

We study the 3D incompressible laminar flow of two interacting phases in a spiral microchannel with a 2-mm radius of curvature using the Euler–Euler framework. The cross section of a spiral channel is rectangular with dimensions of 250 × 150 $μm$⁠. The Eulerian framework is modeled in COMSOL Multiphysics 5.5 (COMSOL Inc.) using the weak form PDE module.^49,50 We employ a structured hexagonal mesh with 2 864 342 cells throughout the 3D geometry, which is constructed using COMSOL Multiphysics. We performed a grid convergence study to obtain the appropriate mesh size for this study (please see the  Appendix). The mesh size is always maintained to be larger than the size of the particles. We also provide a validation of our numerical model by comparing our results with an experimental study on Dean's flow of fluorescent microparticles in a U-shaped channel (please see the  Appendix).

The multiphase flow field is governed by the continuity (mass conservation) and the momentum conservation equations, as described below:^{47,48,51–53}

• Mass conservation equation

  $divui=0,$

  (1)

• Momentum conservation equations

In Eqs. (1)–(3), u_i is the velocity vector of the phase i [i = c (continuous), d (dispersed)], $ρi$ is the density of the pure phase i, $σi$ is the stress tensor associated with phase i, $ρig$ is the body force experienced by phase i, g is the acceleration due to gravity vector, and f_i is the total interaction force experienced by phase i.

We can further write

where $ϕc$ and $ϕd$ are the fractions of the continuous phase and the dispersed phase (solid particle) in each cell, respectively. The stress tensor for the continuous (liquid) and the dispersed phase (solid) is as follows:

where I is the identity matrix of order 3, and $μc$ and $μd$ are the local dynamic viscosities of the continuous phase and dispersed phase, respectively. The viscosity of the continuous phase is described by the Krieger viscosity model as (see Refs. 42 and 54) $μc=μc01−ϕdϕd,max−2.5 ϕd,max$⁠, with $μc0$ being the pure viscosity of the continuous phase, and $ϕd,max$ being the maximum volume fraction (or packing concentration) (⁠$ϕd,max=0.62$ for Krieger viscosity model⁴²) The viscosity of the dispersed phase is given by the Gidaspow solid viscosity model^42,46,55 as $μd=0.5ϕd$⁠. Gidaspow solid viscosity model is generally used in tandem with a solid pressure force (as explained below). The above-mentioned effective solid viscosity, in general, gives us a scaled estimate of the effective viscosity of the solid, and it also helps to achieve numerical robustness;^55,56 on the other hand, the particle–particle interactions (consisting of the inter-particle friction and collision forces) are modeled with solid pressure.^55,57

The interaction forces consist of the interphase momentum transfer (Schiller–Neumann drag) term F_D^58–60 and the Saffman shear lift force F_L.^47,48,61 The total interaction forces experienced by the continuous and the dispersed phases can be, respectively, expressed as follows:^46–48 

Here, $CD(=24Res(1+0.15(Res)0.687))$ is the drag coefficient for the Schiller–Neumann drag model at low Reynolds number (Re_s < 1000), with $Res=ρduc−uddμc$ being the particle Reynolds number. Also, $ρd$ is the density of dispersed phase, $ui$ is the velocity of phase i, d is the particle diameter, and $μd$ is the viscosity of the dispersed phase. Finally, $Dc=12∇uc+∇ucT$ is the symmetrical part of the velocity gradient of continuous phase and tr(D_c) is the trace of the tensor D_c. For the dispersed phase, in addition to the interphase momentum transfer term and the Saffman lift force, a solid pressure force (the first term) exists that accounts for the collision and friction between the particles.^{42,46,54,62,63}

The Mixture theory or the averaging method that is used in this study uses the form of constitutive equations for the stress tensor of the dispersed phase (based on the dispersed phase volume fraction) that is akin to that of the Newtonian fluid. Such an approach has been widely used in several studies.^42,64–67 However, such an approach has certain limitations: for example, (1) such a stress formulation is only complete when complemented with a solid pressure term, (2) there is an obvious error stemming from the fact that the mixture of fluid–solid phase does not behave like a fluid when volume fraction is close to 1,⁶⁸ etc. These issues have been analyzed and addressed by various researchers by (1) incorporating an appropriate solid pressure term and (2) identifying a maximum mixture packing density for which the formulation is employable. In a fluid phase, the pressure term, described by the isotropic part of the stress tensor, is very well defined. On the other hand, such a term defining the solid phase pressure is not well known. Enwald et al., in their work, have defined the solid pressure as the sum of the component due to the particle–particle collisions and the component associated with the presence of fluid pressure gradient.⁴² The solid pressure term corresponding to the particle–particle collisions has been modeled empirically by a number of researchers. In this study, we use one such empirical model, namely, the Gidaspow–Ettehadieh model, to account for the contribution of the particle–particle collision to the solid pressure.⁶⁹ 

Boundary conditions:

The boundary conditions in the presence of which the governing equations [Eqs. (1)–(3)] are solved are as follows:

Here, i represents the phase and it can take values of i = c and i = d for continuous and dispersed phases, respectively.

In Table I, we provide the parameters used in our study for the chosen material combination of plasma (continuous phase in our study) and polystyrene (dispersed phase in our study):

In this section, we investigate the Dean flow characteristics of a two-phase fluid–particle system (consisting of a dispersed phase and a continuous phase) flowing through a spiral microchannel. We provide the analysis of this Dean flow for three separate cases. First, we study the Dean flow for a single-phase fluid flow. Second, we investigate the effect of the particle volume fraction $ϕd,$ which determines the amounts of solid particles in the system, for the case of the two-phase fluid–particle Dean flow. Finally, we analyze the effect of the particle diameter (d), which determines the interaction forces, on the two-phase fluid–particle Dean flow.

Before probing the effect of the presence of the particles and the interaction forces, we first study the flow behavior of the single-phase Newtonian liquid (this liquid serves as the continuous medium for the case of the particle–liquid two-phase flow, discussed later) inside the spiral channel of rectangular cross section and a radius of curvature of 2 mm. For this purpose, we switch off the particle–liquid interaction terms and set the particle volume fraction to be zero. In Fig. 2, we plot the velocity magnitude (|U_c|), the z-velocity (u_cz) of the flow (with arrows providing the direction) in the cross section, and the vortices in the cross section. This cross section is shown in Fig. 1(b): henceforth, the locations x/W = 0 and x/W = 1 correspond to the locations of the outer (concave) and inner (convex) wall of the cross section, while z/H= 0 and z/H = 1 correspond to the locations of the bottom and top walls of the cross section [here W and H denote the width and the height of the microchannel, see Fig. 1(b)]. In a curved rectangular channel, a typical Poiseuille flow profile in the absence of any particles is disturbed due to the presence of centrifugal effects due to the curvature of the channel: this is referred to as the Dean flow effect. Such a centrifugal force causes the shift of the fluid mass toward the concave (outer) wall of the channel,^14,70–73 or effectively, there is a shift of the maximum velocity away from the center of the channel. In a typical Poiseuille's flow, without the effect of any curvature or centrifugal force, the maximum velocity is always at the channel center. The tendency of curvature-driven shift of the maximum fluid velocity is evident in Fig. 2. In Fig. 2, one can clearly see that with a progressive increase in the Dean number (De), the location of the maximum in |U_c| shifts toward the concave (or outer) wall of the channel (identified as x/W = 0). This shift is also confirmed by Fig. 3(a) in which the x coordinate of the location of the maximum velocity progressively decreases (smaller x implies a location closer to the outer wall) with an increase in De. Also, Fig. 3(b), which provides the profile of U_c (corresponding to a given z, i.e., z/H = 0.5) across the channel cross section for different De, confirms that with an increase in De, the location of the peak of the velocity profile shifts more toward the outer wall toward a smaller x value. In Fig. 3(a), we also plot the result for the variation of the x-location of the maximum velocity as reported in the simulation studies of Nivedita et al.¹⁴ The simulation results from Nivedita et al.¹⁴ correlate well with our numerical prediction, thereby providing confidence in our numerical model (with a maximum error of ∼5%).

To better understand the effect of the finite curvature on the overall fluid flow, it is useful to study the corresponding variation of the pressure change $[Δpx=px−pinner$⁠, where p_inner is the pressure at the inner wall or x/W = 1], or the dimensionless variation of $Δpxi.e.,Cp=Δp(x)12ρcUm2$⁠. For a given x (or a given radial location), the velocity peaks at the channel centerline (z/H = 0.5), but progressively decreases to zero due to no-slip condition as one approach the bottom wall (z/H = 0) or the top wall (z/H = 1). Accordingly, the centrifugal force, for a given x, is much smaller near z/H = 0 or z/H = 1, since the velocity is smaller at these locations. On the other hand, as the fluid velocity increases at the centerline (z/H = 0.5) for a given x, the centrifugal force increases at such locations. Therefore, for a given x, to balance this centrifugal force on the fluid mass, one would need a smaller $Δp$ (or a smaller C_p) at locations near the top and bottom walls (i.e., z/H = 1 or z/H = 0), as compared to the $Δp$ (or C_p) at the channel centerline (i.e., z/H = 0.5). Hence, Figs. 4(a) and 4(b) show that for a given x, C_p increases as z/H decreases from 0.9 to 0.5. We can similarly explain the variation of $Δp$ (or C_p) with x by relating the variation of $Δp$ with the corresponding variation of the centrifugal force. For a given z, the centrifugal force increases as one moves toward the outer wall (i.e., as one moves from x/W = 1 to x/W = 0). Therefore, to balance the effect of the centrifugal force acting on the fluid mass, $Δp$ must also increase as one moves from x/W = 1 to x/W = 0, for a given z. This is evident in Figs. 4(a)–4(c). Next, we look at the effect of the Dean number (or the corresponding Reynolds number) on the pressure difference across the channel [see Fig. 4(c)]. Here, we observe an increase in the pressure difference with an increase in the Dean number to balance the corresponding increase in the centrifugal force experienced by the fluid mass.

Next, in Figs. 5(a) and 5(b), we plot the radial velocity (u_cx or the x-component of the velocity) of the single-phase fluid flow. In a curved channel, due to the centrifugal force, the fluid mass tends to move toward the outer wall (x/W ∼ 0). As the flow is contained within an enclosed space, the inbound mass of fluid near the outer wall (x/W ∼ 0) will cause the relatively stagnant fluid that is present near the wall (x/W ∼ 0) to move away. This movement will occur from the center (z/H = 0.5) toward the bottom or top wall (z/H = 0 or z/H = 1), and not in the other direction. This direction of fluid motion is anticipated since it implies the motion of the liquid from a region of higher pressure (z/H = 0.5) to a region of lower pressure (z/H = 0 or z/H =1) [see Figs. 4(a) and 4(b)]. This mass of fluid that has been driven to locations near the bottom and top walls (z/H = 0 or z/H = 1), subsequently, gets directed toward the inner wall (x/W = 1) owing to the existing radial pressure difference between the outer wall (x/W = 0) as compared to the inner wall (x/W= 1). The fluid gets directed toward the outer wall (x/W = 0). This happens for two separate fluid masses between center (z/H = 0.5) and the top wall (z/H = 1) and between the centerline (z/H = 0.5) and the bottom wall (z/H = 0). This motion creates two counter-rotating vortices, also known as Dean vortices, in the two halves of the channel cross section, where one half is between z/H = 0.5 and z/H = 1 and the other half is between z/H = 0.5 and z/H = 0. Such counter-rotating vortices, which are also confirmed in Fig. 2, have been extensively reported in the context of the Dean flow.^74,75

Next, we study the evolution of Dean vortices with the Dean number (see Fig. 2). As explained above, the radial pressure difference across the cross section creates two symmetric counter-rotating vortices. With an increase in De, due to a corresponding increase in the centrifugal force, we observe the eye of the vortex moves toward the outer wall (x/W = 0) [see Figs. 2(a) and 2(b)]. However, at larger De, which gives rise to larger centrifugal force, an increased pressure differential is required to balance the effect of the centrifugal force. The large pressure difference makes it such that two counter-rotating vortices alone cannot balance the opposing forces. Under such circumstances, the formation of additional vortices near the outer (concave) wall occurs, and the fluid mass gets preferentially skewed more toward the outer wall. The presence of these additional vortices, which are observed in Figs. 2(c) and 2(d), is consistent with the study by Boutabaa et al.⁷⁴ and the experimental study by Bara et al.⁷⁵ It is also interesting to note that these additional vortices grow in size with De to compensate for the corresponding increased centrifugal force at larger De [see Fig. 2(d)].

In this subsection, we report the findings from the CFD simulations of the particle–liquid two-phase flow in spiral microchannels. In the  Appendix, we have provided a validation of our CFD model by comparing our findings with the experimental study of the transport of fluorescent microparticles in a U-shaped channel.⁴⁵ Through our CFD simulation results (obtained from the simulations conducted for the spiral microchannels), we compare and analyze the effect of the particle volume fraction on the overall flow by investigating the flow fields in different directions, the induced vortices, the interaction forces experienced by the continuous phase (from the dispersed phase), and the way these interactions are dictated by the velocity fields. We study the cases for different volume fractions of suspended particles $ϕd=0,0.05,0.1,0.15,0.25, and 0.4$ and for different values of the Dean number (De). In addition to the centrifugal force triggered by the curvature of the channel, in a fluid–particle two-phase flow, both the continuous phase (fluid) and the dispersed phase (particle) experience the interphase interaction forces, namely, the inter-phase momentum transfer force [modeled using the Schiller–Neumann drag model, see Eq. (7)] and the shear lift force [modeled using the Saffman lift force model, please see Eq. (7)].

The multiphase flow behavior is probed by quantifying the magnitude of the velocity, the location of the maximum velocity, the structure of the vortices constituting the flow field, the axial (u_cz) and the radial (u_cx) velocity profiles, and the dispersed-phase-continuous-phase interaction forces. In Figs. 6–9, we plot the velocity magnitude (|U_c|), the z-velocity or the axial velocity of the continuous phase (u_cz) (with arrows providing the directions), and streamlines of the flow in the cross-sectional plane to elucidate the vortices [identified in Fig. 1(b)]. In the absence of the interaction forces for the case of the single-phase flow, we observe (see Figs. 2 and 3) that the centrifugal force caused a distinct shift of the fluid mass toward the outer wall (x/W= 0), resulting in larger velocities near the outer wall. This result is also evident for the case of $ϕd=0$ in Figs. 6–9. However, with an increase in the volume fraction of the dispersed phase, we observe a shift in the fluid mass toward the center (i.e., away from x/W = 0). This is caused by the presence of the dispersed-phase-continuous-phase interaction force opposing the centrifugal force acting on the continuous phase. Consequently, the location of the maximum velocity gets shifted away from the outer wall (x/W = 0) and moves toward the center (x/W= 0.5) of the cross section. This is clearly seen in Fig. 10(a), where we plot the location of the maximum velocity magnitude as a function of De and $ϕd$⁠. This is also evident from the plots of the velocity magnitudes along the x direction at a height of z/H = 0.5 for different Dean numbers and volume fractions [please see Figs. 10(b)–10(e)]. For a given De, we clearly see a shift of the location of the peak of the velocity away from the outer wall (x/W = 0) for an increasing volume fraction of the dispersed phase (⁠$ϕd$⁠). This multiphase flow governed by the competitive interplay of the centrifugal and interphase interaction forces could be better understood by investigating the two main components of the interaction forces: the shear-lift force and the interphase momentum exchange (or drag) force.

It is worthwhile to note here that the interphase interaction forces and the velocity field are highly coupled. For instance, the x-component of the drag force is $FDx∝ϕducx$⁠, and the x-component of the lift force is $FLx∝ϕd∇xucy$⁠. Here, $∇x$ denotes the gradient in x direction and u_cy is the velocity of the continuous phase in y direction. Therefore, it is important to probe the components of the velocity and the gradients of velocity components to understand the interphase interaction forces. First, in Fig. 11, we plot the x-component of the velocity of the continuous phase (⁠$ucx$⁠) for different values of $ϕd$ and the Dean number. $ucx$ is found to be positive (i.e., directed along the positive x axis) near the top and the bottom walls (z/H ∼ 1 and z/H ∼ 0), and negative (i.e., directed along the negative x axis) near the center (z/H ∼ 0.5). This is seen for the case of single-phase flow as well (see Fig. 5). This flow profile can be attributed to the centrifugal force acting on the fluid, which causes the fluid mass to move toward the outer wall (x/W = 0). However, very much like the case where no particles are present, here too as the fluid is flowing in an enclosed space (in x direction), the fluid moving toward the outer wall causes the relatively stagnant fluid mass near the outer wall (x/W ∼ 0) to move away from the center (z/H ∼ 0.5) (such a movement implies a negative u_cx near z/H ∼ 0.5). This fluid mass moving away from the center (z/H ∼ 0.5) gets recirculated toward the inner wall (x/W ∼ 1) at locations near the top (z/H ∼ 1) and bottom walls (z/H ∼ 0) (such a movement implies a positive u_cx near z/H ∼ 0 and z/H ∼ 1). This explains the u_cx velocity profile reported in Fig. 11 for the different combinations of the particle volume fractions and the Dean number. We also note that the magnitude of the $ucx$ increases with an increase in the Dean number (or Reynolds number) for a given volume fraction of the dispersed phase. This is because the increase in Reynolds number leads to an increase in the centrifugal force, which in turn increases the magnitude of $ucx$ (which is also the radial velocity). We also observe that increasing the volume fraction of the dispersed phase (⁠$ϕd$⁠) reduces the magnitude of the radial velocity (⁠$ucx$⁠). This can be attributed to the net interphase interaction force, which opposes the centrifugal force that causes the shift of fluid mass toward the outer wall (x/W ∼ 0). The components of the interphase interaction forces and their quantitative dependence on $ϕd$ will be discussed later.

In Fig. 12, we plot the y-velocity component (⁠$ucy$⁠) of the continuous phase across the radial direction (across x) at three different heights z/H = 0.5, 0.7, and 0.9, in addition to plotting the x-gradient of this velocity component $∇xucy=∂ucy∂x$ at locations near the outer wall (x/W = 0.2) and the inner wall (x/W = 0.8). The distribution of y-velocity is very similar to that of the velocity magnitude plotted in Fig. 10, as y-velocity is the most dominant velocity component. It can be seen that $∂ucy∂x$ is generally negative at locations near the outer wall (x/W = 0) and positive at locations near the inner wall (x/W = 1). This is due to the fact that for a given cross section (xz plane) at a given height (or given z), the transverse velocity (⁠$ucy$⁠) is zero at the side walls (x/W = 0,1) showing a parabolic (or skewed parabolic) profile (of negative velocity) depending on the Dean number and the volume fraction. The skewness of the parabolic profile is due to the centrifugal force pushing the mass of fluid toward the outer wall (x/W = 0). It can also be seen that the magnitude of the velocity gradient $∂ucy∂x$ increases with an increase in the Dean number (or Reynolds number). This is clearly due to the increase in the overall velocity, which in turn increases the gradient (stemming from the fact that the velocity at the walls is always zero). Next, we observe a non-monotonic variation of the gradient with respect to the volume fraction of the dispersed phase. This can be attributed to the effect of the net interphase interaction force, which will be explained later. It is interesting to note that the gradient becomes positive at a few locations (z/H ∼ 0.5), especially at higher Dean number and smaller values of $ϕd$⁠. This is due to lower volume fractions and weaker opposing (i.e., opposing the centrifugal forces) interphase interactions. Also, at higher Dean numbers, the centrifugal force, shifting the fluid mass outwards, is much larger. The combined effect of these two factors results in a maximum velocity (magnitude) of the y-component of the velocity at locations much closer to the outer wall. For these cases, the negative peak of the velocity is at an x-location below x/W = 0.2, and hence, we see a positive gradient of y-velocity as the slope is increasing at the location of x/W = 0.2.

Next, in Figs. 13(a-i)–13(a-iv), we plot the x-component of the normalized interphase momentum exchange force [or normalized drag force, $FDx/ρcUm2/R$] for different Dean numbers as a function of the volume fraction of the dispersed phase (⁠$ϕd$⁠) at x/W = 0.2 across the height of the channel. The interphase drag force exerted on the particle by the fluid causes the particles to move in the direction of the flow. Therefore, the direction of the drag force exerted by the particle on the fluid is opposite to the direction of the flow at any location. This is evident from Fig. 13(a), where we see that the profile of the x-component of the drag force on the continuous phase is exactly opposite the profile of the typical radial flow velocity (u_cx) seen in Sec. III A as well as the profiles shown in Fig. 11. For smaller Dean numbers (De = 34, 68), we witness a monotonic increase in the drag force with an increase in the volume fraction of the dispersed phase (⁠$ϕd$⁠). This is consistent with the monotonic increase in the radial velocity with increasing volume fraction of the dispersed phase (⁠$ϕd$⁠). However, we witness a non-monotonic variation of the drag force magnitude (increasing volume fraction) at locations closer to the wall (x/W= 0.2), especially for larger Dean number (De = 206, 343). Such a non-monotonic variation of $FDx$ with increasing volume fraction is due to the combined effect of non-monotonic variation of velocity and the increase in the volume fraction. It is important to note that the direction of the x-component of the drag force depends only on the z-location. At all x-locations, the drag force on the continuous phase is positive at locations near the center (z/H ∼ 0.5) and negative at locations near the top and bottom walls (z/H = 1 and z/H = 0): this profile is exactly opposite to the profile of the x-component of the continuous-phase velocity (⁠$ucx$⁠).

Finally, in Figs. 13(b-i)–13(b-iv), we plot the variation of x-component of the normalized shear-lift force (experienced by the continuous phase) for different Dean numbers as a function of the volume fraction of the dispersed phase (⁠$ϕd$⁠) at a location of x/W = 0.2 across the height of the channel. The shear-lift force (magnitude) experienced by the particle is directly proportional to the volume fraction of the solid particles and the x-gradient of the y-component of the continuous phase velocity, that is, $FLx∝ϕd∇xucy$ [please see Eq. (7)]. It is interesting to note that the x-gradients of other velocity components also contribute to the lift force; however, we have only shown the dominant component, namely, $∇xucy,$ that contributes to the lift force. The shear-lift force on the particle is directed toward the location where the relative velocity between the particle and the continuous phase is maximum. For a particle moving with the fluid in a channel flow, the maximum relative velocity is near the wall (since at the wall the fluid velocity is zero due to the no slip and no-penetration conditions), and hence, the particle gets directed toward the wall.^32,55,76 This force on the particle (dispersed phase) comes from the continuous phase. Therefore, a force that is equal in magnitude but opposite in direction gets applied to the fluid (continuous phase) by the particles (dispersed phase). Hence, the continuous phase must experience a shear-lift force that is directed away from the wall. In other words, the direction of the shear-lift force is opposite to the direction of the x-gradient of y-velocity component (see Fig. 12). This is evident from Figs. 13(b-i)–13(b-iv), where we observe a shear lift force component (force in x-direction) experienced by the continuous phase that is always directed away from the wall (x/W = 0) and toward the center (x/W = 0.5). Therefore, there is a positive x-force at location x/W = 0.2 and a negative x-force at the location x/W = 0.8 [see Fig. 13(b-iii), inset]. Also, we observe a monotonic increase in the lift force with increasing volume fraction, since an increase in the number of particles results in an increase in the force exerted on the fluid by these particles (⁠$FL∝ϕd$⁠), for all the different chosen Dean number values. In Fig. 13(c), where we provide a schematic representation of the typical direction of the different interaction forces acting on the dispersed phase (particle) and the continuous phase (fluid), this shear-lift force and the drag force (discussed previously) have been depicted.

Similar to the single-phase flow, in the multiphase flow we also observe two symmetric counter-rotating vortices across the channel cross section (see the right column in Figs. 6–9). We observe that with an increase in the volume fraction, which causes an increase in the net interaction forces, there occurs a shift in the eye of the vortices from the outer wall toward the inner wall. It is also interesting to observe that the splitting of vortices to balance the larger centrifugal force, which is shown for the single-phase flow at higher Dean number (De = 206, 343), is no longer observed for the case of the two-phase particle–liquid transport for such Dean number values. This again can be ascribed to the presence of an increased interaction force opposing the centrifugal force.

In this subsection, we report the findings from the CFD simulations of the effect of the particle diameter on the flow profiles for the particle–liquid two-phase transport in spiral microchannels. Since many of the applications for spiral microchannel designs involve blood cell separation, we consider the cases of three different particle diameters $d=0.5, 1,2μm$ that are good representative sizes for platelets, red blood cells, and circulating tumor cells, respectively, for different Dean number values. While the Dean number affects the centrifugal force induced by the curvature of the channel, the change in the particle diameter impacts the interphase interaction forces, namely, the inter-phase momentum transfer force (modeled using Schiller–Neumann drag model) and the shear lift force (modeled using Saffman lift force model).

In Figs. 14–17, we plot the velocity magnitude |U_c|, axial velocity (z-component velocity) of the continuous phase (arrows indicating the direction of the flow), and the streamline (showing the counter-rotating vortices) in a cross-sectional plane shown in Fig. 1(b) elucidating the effect of the particle size for a given particle volume fraction (⁠$ϕd=0.15$⁠). It is clear from these plots (Figs. 14–17) that the location of the maximum velocity magnitude moves further away from the outer wall (x/W = 0) with an increase in the particle diameter. This means the fluid mass is directed further away from the outer wall with increasing particle diameter. Therefore, an increase in the particle diameter leads to an increase in the particle-induced interaction forces on the continuous phase that progressively oppose the curvature induced centrifugal force (on the continuous phase). Later, we shall see that it is primarily the shear-induced lift force that is augmented with the increase in the particle size, enforcing a greater nullification of the effect of the centrifugal force. It is also clear from Fig. 18(a), where we plot the location of the maximum velocity magnitude. Furthermore, the plots of the velocity magnitude at the location of z/H= 0.5 across the cross section also shows a distinctive shift in the location of the peak of the velocity magnitude away from the outer wall (x/W = 0) with an increasing particle diameter d [see Figs. 18(b)–18(e)]. This behavior could be understood by investigating the net effect of the centrifugal forces and the particle size-dependent interphase interaction forces.

It is essential to probe the flow velocity and its gradients to understand the components of the interaction forces. For instance, we see that the x-component of the drag force $FDx∝ucxd2$ and the x-component of the lift force $FLx∝∇xucyd$⁠. First, in Fig. 19, we probe the radial velocity component (⁠$ucx$⁠) of the continuous phase for various particle sizes and Dean numbers. It can be observed that $ucx$ for a given x-location, for smaller particle sizes, is directed toward the outer wall near the centerline (z/H = 0.5) and directed toward the inner wall at locations near the top and the bottom walls (z/H = 1 and z/H = 0). As elucidated in Sec. III A, this behavior can be attributed to two factors. First is the nature of the centrifugal force acting on the fluid due to the curvature of the channel. The second factor is that the flow in the radial direction is occurring in an enclosed space. Therefore, the incoming fluid mass (toward the outer wall), due to the centrifugal force, causes the relatively stagnant fluid mass to move away from the centerline and eventually move toward the inner wall. For all the chosen values of De, the net effect of the increase in the particle diameter is to push the fluid further away from the outer wall. The components of the net interaction forces causing such a shift would be discussed in detail later.

Next, we investigate the z-variation of the x-gradient of the y-component of velocity of the continuous phase, that is, $∂ucy∂x$⁠, for various Dean numbers and particle sizes (⁠$d$⁠) (see Fig. 20). In Fig. 20, in addition to plotting z-variation of $∂ucy∂x$ at two specific x locations (x/W= 0.2 and x/W = 0.8), we also plot the radial profile of the y-velocity component of the continuous phase (u_cy) at different heights (z/H= 0.5, 0.7, and 0.9). Here, we witness a negative gradient at locations near the outer wall (x/W = 0.2), and a positive gradient at locations near the inner wall (x/W = 0.8) for most of the chosen parameters. As explained in Sec. III B, the combined effect of the centrifugal force and the interphase interaction force results in a skewed parabolic profile [see the y-velocity profiles in Figs. 20(a)–20(d)] causing such positive and negative gradients at different x locations. It can also be seen that the profile of $ucy$ is similar to that of the velocity magnitude (see Fig. 18), as $ucy$ is the dominant velocity component. The effect of the variation of the Dean number on $∂ucy∂x$ is clearly seen in Fig. 20, where we witness an increase in the gradient with an increase in the Dean number. This is due to the increase in the velocity, which in turn increases the gradient for a given channel dimension (due to the fact that the velocity is zero at the walls). The effect of the particle diameter seems to have a non-monotonic effect especially at larger Dean number values. This can be associated with the combined effect of the centrifugal and the interphase interaction forces (see below for a more detailed discussion).

Next, we probe the radial (or x) components of the interaction forces. First, in Figs. 21(a-i)–21(a-iv) we plot the z-variation of radial component of the normalized interphase drag force modeled using the Schiller–Neumann model for different particle sizes and Dean number for a dispersed-phase volume fraction of $ϕd=0.15$ at x/W = 0.2. As explained in Sec. III B, the drag force imparted on the particle by the fluid is in the same direction as the flow, and hence, the drag force experienced by the liquid due to the particle is equal in magnitude and opposite in direction. This is clearly observed in Figs. 21(a-i)–21(a-iv), where the direction of the x-component of the drag force on the liquid is opposite to that of $ucx$ (see Fig. 19). It is clear from Figs. 21(a-i)–21(a-iv) that the x-component of the normalized interphase drag force monotonically decreases with an increase in the particle size. This is consistent with the monotonic decrease in the radial velocity (u_cx) with increasing particle size and the fact that $FDx∝ucxd2.$

In Figs. 21(b-i)–21(b-iv), we plot the variation of the x-component of the normalized shear-lift force for different Dean numbers for varying particle size (⁠$d$⁠) at a x-location of x/W = 0.2. The magnitude of the shear-lift is directly proportional to the x-gradient of velocity (⁠$∂ucy∂x$⁠) and inversely proportional to the particle size (d). In addition to this, the shear-lift force is also dependent on the relative x-velocity component (⁠$ucx−udx$⁠), that is, $FLx∝∇xucyducx−udx$⁠. First, it is worthwhile to note that the direction of the shear-lift acting on the particle is in the same direction as $∂ucy∂x$⁠. Hence, the direction of the shear-lift force acting on the fluid is in the direction opposite to that of $∂ucy∂x:$ this is clear from Figs. 20 and 21. Also, it can be seen from Fig. 21(b) that the magnitude of shear-lift force increases monotonically with the particle size. This is because of the following reasons. First, with increasing particle sizes, we witness a decrease in $∂ucy∂x$ (for lower De; see Fig. 20). Also, $FLx∝1d$⁠. These two factors combine to ensure that an increase in the particle diameter (d) will decrease $FLx$⁠. On the other hand, we find that that the relative x-velocity component, $ucx−udx,$ increases with the particle diameter (d), which indicates that larger particles have a larger slip velocity. This can be seen from the inset of Fig. 21(b-i) where we plot the product $∇xucyucx−udx$⁠. This competition between the decrease in $∇xucyd$ and the increase in $ucx−udx$ with the particle diameter results in a monotonic trend of the lift force with d, as witnessed in Fig. 21. Such a particle size dependence of $FLx$ supports our observation (for a wide range of De values) of a sharp shift in the location of the maximum fluid velocity from a location in the vicinity of the outer wall (x/W= 0) to a location significantly away from the outer wall for the case when the particle diameter (d) increases from $0.5$ to $1 μm$ (see Figs. 14–18), while such a shift is less distinct as the particle diameter (d) increases from $1$ to $2 μm.$

Next, we investigate the evolution of the vortices for different Dean numbers and particle diameters of the dispersed phase. We observe that with an increase in the particle diameter, the net interaction forces increase to cause a shift in the eye of the vortices from the outer wall toward the inner wall (see Figs. 14–17). It is also interesting to observe that the splitting of vortices (or the formation of additional or secondary vortices), typically witnessed at higher Dean numbers (De = 206, 343) to compensate for higher centrifugal force in a single-phase flow (see Fig. 2), is no longer observed for any given particle diameter for the chosen volume fraction (⁠$ϕd=0.15$⁠) in the multiphase flow (see Figs. 14–17). This, too, can be ascribed to the effect of the increased interaction force opposing the centrifugal force.

In this paper, a CFD strategy has been employed for computationally studying for the first time the two-phase liquid–particle transport in spiral microchannels giving rise to inherent centrifugal forces and Dean flow effect. Therefore, our study produces the first set of original fluid mechanics results in the context of the mechanistic understanding of the behavior of the two-phase liquid–particle transport in microfluidic Dean flow in spiral microchannels.

Curvature-induced forces have been known to cause a large asymmetry in the flow profile, skewing the flow toward the outer wall of the microchannel. This type of flow profile also contains two counter-rotating vortices across the microchannel cross section. Increasing Dean numbers increase the extent of asymmetry and also distort these vortices by generating secondary vortices. Our CFD platform, validated against the experimental findings of Nivedita et al.,¹⁴ confirms these flow physics aspects associated with the single-phase Dean flow in spiral microchannels.

Subsequently, our CFD framework is used to probe the particle–liquid two-phase transport in such spiral microchannels. The presence of the particles invariably reduces the asymmetry-inducing tendency of the Dean flow. These particles ensure that the liquid is now subjected to the drag and the lift forces from the particles in reaction to the drag and lift forces that the fluid flow exerts on the particles. These forces on the particles counter and diminish the effect of the curvature-driven centrifugal force, leading to a dampening of the asymmetry-inducing effect of the Dean flow. Furthermore, the generation of the secondary vortices, a characteristic of the profile of the single-phase microfluidic Dean flow at large Dean number values, is also no longer observed when the fluid is seeded with particles. We investigate the cases corresponding to different particle volume fractions and different spherical particle diameters. For each case, we quantify the drag and the lift forces exerted on the continuous phase and the way these forces could be related to the corresponding velocity profiles. This analysis allows us to effectively quantify the impact of particles in spiral microchannels and the coupled particle–liquid transport in two-phase Dean flow.

This work was supported by the U.S. National Science Foundation (Award No. 1935814).

Please note, the mention of commercial products, their sources, or their use in connection with material reported herein is not to be constructed as either an actual or implied endorsement of such products by the U.S. Department of Health and Human Services. The methodology discussed here also does not constitute any regulatory policy endorsed by the U.S. Food and Drug Administration. Nor is it intended for replacing any existing policies pertaining to submissions that are currently under review at the U.S. Food and Drug Administration.

The authors have no conflicts to disclose.

The data that support the findings of this study are available within the article.

In this section, we provide the results of the grid convergence study performed to estimate the appropriate the mesh size required. For this purpose, we took the case with large velocity gradients and secondary flow component, that is, the case with a dispersed phase volume fraction of $ϕd=0.05$ and Dean number De = 343. We compare the velocity magnitude of the continuous phase $Uc$⁠, radial velocity component (⁠$ucx$⁠) and x-component of the normalized interaction forces for three different hexagonal meshes of different grid sizes. From the Fig. 22, it can be seen that the grid convergence occurs with grid 3.

In this section, we attempt to validate our Eulerian multiphase model against the experimental findings of Ramachandraiah et al.⁴⁵ In that study, the authors quantified the lateral positioning of the fluorescent microparticles due to the inertial effects in the presence of very small Reynolds number flow in a curved channel. They studied the Dean flow in a U-shaped channel with different curvatures, channel width, and Dean numbers (De). They quantify the location of the particle by capturing the fluorescence intensity at each location. For the purpose of validation of our model, we chose the case of De = 36 inside a U-shaped channel with a width of $250 μm$ and curvature of 2 mm.

In Fig. 23, we plot and compare the normalized volume fraction denoting the concentration of the dispersed phase obtained from the simulation (conducted by employing our Euler–Euler model for the two-phase flow in a U-shaped channel of same dimensions as considered by Ref. 45) with the fluorescent intensity (a measure of concentration) of the particles (reported by Ref. 45). The experimental data are extracted from Ref. 43 using GRABIT digitizer function in MATLAB.⁷⁷ We do not have a direct measure of the volume fraction of the particles for this experimental study.⁴⁵ However, the fluorescent intensity distribution across the channel provides a good estimate of the concentration of the particles (assuming that this intensity is proportional to the particle concentration). To obtain a direct comparison between the fluorescent intensity and the volume fraction (obtained from the simulation), we rescale the volume fraction (obtained from the simulation) to a scale from 0 to 1. This is done using the following relation to obtain the normalized volume fraction, $ϕn=ϕd−minϕdmax(ϕd)−minϕd$⁠. $ϕn$ shows a good quantitative match with the experimentally observed fluorescent intensity of Ref. 45 (see Fig. 23). It is to be noted here that we validate our simulation model against the findings of an experimental study⁴⁵ that considers relatively dilute concentration of particles; this stems from the apparent lack of experimental studies of Dean's flow with significantly large particle concentration.

In Fig. 23, we notice some differences between the experimental and the simulation-predicted (based on equivalent concentration values) intensity values. These discrepancies in the values could have stemmed from the following reasons: (1) approximations in the parameters and properties that were not explicitly reported in the experimental work (e.g., in the experimental study, only the range of the particle volume fraction was reported instead of the actual value), (2) approximations involved in modeling such particle flows using empirical models describing particle behaviors, and (3) errors associated with the fluorescent light-based experiments (such as reflectance) that were used to generate the intensity values.