Numerical investigation of elasto-inertial particle focusing patterns in viscoelastic microfluidic devices

TL;DR

This paper uses 3D computational simulations to explore how particles behave in viscoelastic microchannels, revealing new particle dynamics and providing a phase diagram for device design.

Contribution

It offers the first detailed 3D analysis of particle focusing in viscoelastic flows considering inertial and elastic effects, including novel particle behaviors.

Findings

Identified equilibrium positions for particles across various flow conditions

Provided a phase diagram predicting particle dynamics

Discovered new particle behaviors not previously observed

Abstract

Viscoelastic microfluidic devices are promising for various microscale procedures such as particle sorting and separation. In this study, we perform three dimensional computational investigation of a particle in a viscoelastic channel flow by considering combined inertial and elastic effects. We calculate equilibrium positions of the particle in the microchannel for a wide range of Reynolds and Weissenberg numbers, which are important in the design of microfluidic devices. The results provide an insight into the motion of cells and particles and explain the findings of previous experiments. Furthermore, we suggest new particle behaviors that have not been found before. Ultimately, a phase diagram is provided to predict the particle dynamics under a wide range of inertial and elastic effects.

Equations

$$\nabla .\text{u}=0,$$

$$\frac{\partial(\rho\textbf{u})}{\partial t}+\nabla.(\rho\textbf{uu})=-\nabla p+\nabla.\mbox{\boldmath$\tau$}+\textbf{F},$$

$$\mbox{\boldmath$\tau$}=\mbox{\boldmath$\tau_{s}$}+\mbox{\boldmath$\tau_{p}$}.$$

$$\lambda\overset{\bigtriangledown}{\mbox{\boldmath$\tau_{p}$}}+\mbox{\boldmath$\tau_{p}$}+\frac{\alpha\lambda}{\mu_{p}}\mbox{\boldmath$\tau_{p}$}.\mbox{\boldmath$\tau_{p}$}=\mu_{p}(\nabla\textbf{u}+\nabla\textbf{u}^{T}),$$

$$\overset{\bigtriangledown}{\mbox{\boldmath$\tau_{p}$}}\equiv\frac{\partial\mbox{\boldmath$\tau_{p}$}}{\partial t}+\textbf{u}.\nabla\mbox{\boldmath$\tau_{p}$}-\nabla\textbf{u}\mbox{\boldmath$\tau_{p}$}-\mbox{\boldmath$\tau_{p}$}\nabla\textbf{u}^{T}.$$

$${\text{F}}_{total}=-{\oint }_V\text{F}dv,$$

$$\textbf{F}_{el}=\int\mbox{\boldmath$\tau_{p}$}.\textbf{n}\,ds,$$

$${\text{F}}_{in}={\text{F}}_{total}-{\text{F}}_{el},$$

Full text

\alsoaffiliation

School of Mechanical Engineering, Purdue University, West Lafayette, IN 47907, USA

1 Introduction

Microfluidic devices have been widely used in biomedical devices. These devices have significantly enhanced therapeutic, diagnostic and many industrial procedures by increasing accuracy and accelerating the processes1, 2, 3, 4, 5, 6, 7. Isolation and separation of rare cells from a heterogeneous population of cells is a critical process in early diagnosis of fatal diseases such as cancer8 and malaria9. Furthermore, the isolated cells such as rare blood components can be used for therapeutic purposes. For instance, platelet-rich plasma are used in transfusion10, 1 or stem cells can be found in blood samples1. Enriching the cell population provides a platform to biologists to conduct physical and chemical analysis on cells11, 12, 13, 14, 15, 16.

In order to accomplish the desired tasks in aforementioned applications the precise control of particles is required17. Hence, developing new techniques for the manipulation of particle trajectory has been the subject of many studies in the past decade. Some of these techniques are designed based on the use of externally applied forces generated by electric18, magnetic19 and acoustic20 fields. These methods offer high sample processing rates1, however, there are many factors that prevent them from widely being used in clinical applications. Mainly, these methods work based on biochemical labeling of the cells that may affect the cell function and properties2, 1. Furthermore, the cost and complexity of the process can also be considered as an important downside for these methods1. Hence, there is a growing interest in developing label-free techniques that introduce advantages such as accurate analysis, low sample use and low cost operations2. In this regard, inertial microfluidic devices are used to manipulate particle trajectories in microchannels.

Segre and Silberberg21 first reported the transverse migration of particles due to inertial forces in a straight tube filled by a Newtonian fluid21. The hydrodynamic interaction between the flowing particles and the ambient fluid can be used to manipulate the trajectory of particles. This effect was studied theoretically by 22 in which they considered the dynamics of a rigid sphere in a 2D Poiseuille flow and a simple shear flow and predicted the equilibrium locations of the particle in these flows. Furthermore, they showed that the final location of the particle is independent of releasing point in the microchannel. 23 also theoretically investigated the lateral migration of a solid particle in a circular tube. They studied the motion of a buoyant and neutrally buoyant particle in a Poiseuille flow and showed that their findings agree well with experiments. Many microfluidic setups are designed based on inertial migration of particles to control the location of targeted cells in microchannels, particularly for sorting and separation of the cells24, 25, 26, 27, 28, 29. The efficiency of this method reduces in cases where the targeted cells are small or the flow rate is low30, 31, 17. Viscoelastic microfluics are developed to address this issue by replacing the Newtonian ambient fluid by a dilute polymer solution30. The solute polymer chains get deformed in the induced flow field and exert an additional elastic force on the particle. The resulting force affects the particle migration along with the inertial force17, 25. Previous studies show that the direction and the magnitude of generated elastic force depends on the rheology of the polymer solution and the volumetric flow rate of the suspending fluid30, 32, 33.

Transverse migration of a particle suspended in a viscoelastic fluid is caused by the lift force generated due to the interaction between ambient flow and the particle34, 35, 36. This lift force comprises an inertial lift force ($F_{in}$) and an elastic lift force ($F_{el}$). The inertial lift force can be decomposed into two forces: (i) shear-gradient lift force ($F_{s\_in}$) that arises from the non-uniform velocity profile across the channel and drives the particle away from the channel center and (ii) wall-induced lift force ($F_{w\_in}$) that is caused by the uneven distribution of the vorticity around the particle that leads to higher pressure in the gap between the wall and the particle and pushes the particle away from the wall37, 22, 38, 39. These forces have been investigated extensively in the literature and there are experimental40, 41, numerical42, 43, 44, 45 and analytical46, 22, 35, 36 works proposing scaling relations for the total inertial force in square and rectangular microchannels. On the other hand, the elastic force acts on the particle due to non-uniform distribution of normal stress difference across the channel 34, 47. In this phenomenon, the first normal stress difference ($N_{1}$) generates a stream-wise tension and the second normal stress difference ($N_{2}$) gives rise to a secondary flow in the cross section48. There are numerical49, 48, 50, 51, 52 and experimental17, 53 studies investigating the particle migration in a viscoelastic fluid. Here, in this work we focus on illustrating the distribution of elastic and inertial forces acting on a particle in a viscoelastic channel flow and calculating the equilibrium positions of the particle in the channel cross section.

In order to understand the mechanism of particle migration in a viscoelastic fluid, fully resolved 3D numerical simulations are conducted. In this work, we show the distribution of lift force acting on the particle in a viscoelastic fluid and investigate the influence of combined elastic and inertial forces on the particle behavior in a microchannel. Furthermore, the location of equilibrium points and their corresponding stability are determined for $1<Re<30$ and $0<Wi<3$ which is important for designing the microfluidic devices relying on viscoelastic effects. Our results explain various focusing pattern of particles observed in previous experimental works by scrutinizing the stability of equilibrium points. We also predict new behaviors that have not been discovered in the past studies.

2 Mathematical modeling

In this work, we study the dynamics of a solid spherical particle suspended in a straight, square microchannel. The particle is neutrally buoyant and the ambient fluid is viscoelastic. The radius of the sphere is represented by $a$ and the channel side and its length are $2w$ and $L$, respectively. In this problem, the origin of the reference frame is located at the channel center and $x$, $y$ and $z$ directions are aligned with the streamwise, horizontal and vertical directions, respectively, as illustrated in Fig. 1. The particle is initially at rest and a constant pressure gradient is applied in the $x$ direction, driving the flow in the microchannel. In order to simulate the hydrodynamic interaction between the particle and the surrounding fluid, an incompressible Navier-Stokes equation is numerically solved as follows:

[TABLE]

[TABLE]

where $\rho$ is the fluid density, u is the velocity vector, $t$ denotes the time and $p$ and $\tau$ represent the pressure and stress tensor, respectively. The particle density is assumed to be the same as the fluid density. A Distributed Lagrangian Multiplier (DLM) method is used to simulate the solid particle motion in the fluid. A forcing term F is added in eq. 2 to enforce the rigid body motion of the particle. The details of DLM method can be found in 54 54 The viscoelastic properties of the fluid can be modeled by splitting the stress tensor into two parts: (i) the contribution from solvent $\tau_{s}$ and (ii) that of polymer $\tau_{p}$. Hence, the total stress tensor can be written as:

[TABLE]

The Newtonian viscous stress is described as \mbox{\boldmath\tau_{s}}=\mu_{s}(\nabla\textbf{u}+\nabla\textbf{u}^{T}), where $\mu_{s}$ represents the solvent viscosity. The Giesekus constitutive equation 55 is used to model the viscoelastic behavior of the fluid. This model captures the shear-thinning behavior and constrained elongation of polymer chains in the fluid56. According to this model, the polymeric stress tensor is governed by:

[TABLE]

[TABLE]

Here $\lambda$ represents the fluid relaxation time and $\mu_{p}$ and $\alpha$ denote the polymeric viscosity and the mobility factor. The following equations are used to calculate the total, elastic and inertial forces acting on the particle.

[TABLE]

[TABLE]

[TABLE]

where n represents the unit vector normal to the particle surface and $ds$ and $dv$ denote differential elements of particle surface and volume, respectively. In this problem, the no-slip boundary condition is applied at the walls of the channel. These boundaries are normal to y and z directions. The periodic boundary condition is applied at the inlet and outlet of the computational domain. These boundaries are normal to the x direction. The microchannel length is set to $L=20a$ to ensure the channel is long enough and the particle does not interact with its periodic image. The details of the numerical methods used in this work and their validations are reported in our previous works57, 50, 51, 56.

We use $w$ as the characteristic length scale and $U_{0}$ as the characteristic velocity scale (where $U_{0}$ is the undisturbed flow velocity at the channel center filled with a Newtonian fluid). Accordingly, the governing dimensionless numbers can be defined as: (i) $Re=\frac{\rho U_{0}2w}{\mu}$, representing the ratio of inertial to the viscous forces in which the total viscosity is defined as $\mu=\mu_{s}+\mu_{p}$, (ii) the Weissenberg number $Wi=\frac{\lambda U_{0}}{w}$, representing the ratio of elastic to the viscous forces, (iii) $\beta=\frac{\mu_{p}}{\mu}$, representing the ratio of polymer viscosity to the total viscosity and (iv) blockage ratio $\frac{a}{w}$, describing the finite size of the flowing particle. In this work, $\beta$, $\frac{a}{w}$ and $\alpha$ are set to $0.5$, $0.3$ and $0.2$, respectively, unless otherwise stated. The particle has the same density as that of the fluid. It should be noted that the initial velocity of the fluid is set to the undisturbed flow velocity in the absence of the particle. In order to find the lift force distribution experienced by the particle across the channel cross section, the lateral position of the particle is fixed at the location where the lift force should be calculated. Consequently, the particle only travels along a line parallel to the $x$ direction and it rotates freely around all directions. The particle is consequently not force-free and the elastic and inertial forces balance with the transverse force required to fix the transverse location of the particle. Hence, the dynamics of this particle is not the same as a force-free particle dynamics due to non-linearity of the problem, but the equilibrium positions and their instability are the same as the ones for a force-free particle. The particle is released with a zero initial velocity in the microchannel and the simulation continues until the ambient flow field reaches a steady state and then the lift force is calculated. This method is used to determine the particle’s equilibrium locations and their corresponding stability across the microchannel in a Newtonian fluid58, 45, 59, 60, 61, 44. The domain and grid size independency tests and the validation of the employed method are shown in the supplementary material.

3 Results and discussion

The migration of particles in Newtonian fluids has been extensively investigated by many researchers. According to the previous numerical and experimental studies21, 5, the particles released in microchannel with a circular cross-section migrate toward an annulus ring with a radius of $\sim 0.6R$, where $R$ represents the radius of microchannel. However, the particle dynamics changes in a square microchannel due to the reduced level of symmetry in the flow field56, 5, 53. In this case, the particles focus at four discrete points along the main axes of the microchannel. Considering the flow structure, there are nine equilibrium points, where the lateral force becomes zero, among which only four are stable. The rest exhibit unstable behavior, meaning the particles migrate away upon any disturbance in the flow58, 45, 5. This phenomenon is further explained in details in the following section. In this work, we investigate the stability and location of equilibrium points for the particles suspended in a viscoelastic fluid and we show the focusing patterns of particles for a wide range of $Re$ and $Wi$ numbers. The variety observed in the focusing patterns of particles indicates the promising effect of addition of polymers in microfluidic devices and we show that this method can be used for different applications without applying any changes to the geometry of the microchannel.

3.1 Migration in a low inertial regime

In order to find the location and stability of the equilibrium points in a viscoelastic fluid, the force field experienced by the particle is calculated. The particle size is fixed at $\frac{a}{w}=0.3$ and the Reynolds number is set to $Re=5$ in this section. Due to the symmetry of the problem, the force field is illustrated only for one quarter of the channel cross section. Figure 2(a) shows the lateral force profile acting on the particle at $Re=5$ and $Wi=0$.

In this figure, the location of the equilibrium points are indicated by red circles. The observed profile explains the particle configuration observed in 26, 62 and 45 As shown in this figure, the radially directed forces drive the particle away from the center and the wall and push the particle toward an annulus ring across the channel marked with the dashed line. Hence, the particle primarily moves in the radial direction to reach the annulus ring. This motion is followed by the migration along the ring to focus at its equilibrium position along the main axes58. Considering the force-map illustrated in Fig. 2(a), there are four stable equilibrium points along the main axes and four unstable equilibrium points along the diagonal of the channel. The results indicate that the lateral force is also zero at the channel center, however, the force field around the channel center implies that this equilibrium point is unstable. The focusing pattern of particles for this force map that can be observed in experiments is presented in Fig. 2(b). In order to further investigate the location and stability of the equilibrium points, the lateral force profile along the y-direction (main axis of the microchannel, i.e., $z=0$ and $y=0$ lines in the cross section) is shown in Fig. 3 for $Re=1$ and $5$ and a wide range of $Wi$ number.

According to Fig. 3(b), the force profile for $Re=5$ and $Wi=0$ crosses the dashed horizontal line at two points (one at the center and the other at $\sim 0.52w$). These points are identified as equilibrium points on the main axis and correspond to the locations where the lateral force is zero in Fig. 2(a). Stability of equilibrium points depends on the slope of the force profile at those locations. Hence, the center of channel is an unstable equilibrium point due to its positive slope, while the off-center equilibrium point is stable due to the negative slope of the force profile in that region. Increasing the $Wi$ number changes the force field and consequently alters the focusing pattern of the particle. Figure 2(c) shows the lateral force acting on the particle for $Re=5$ and $Wi=0.1$. In the region constrained between two dashed curves the particle moves toward the center, while in the outer regions the particle moves toward the walls. A significant change in the force direction can be observed in the near wall region, where it is directed toward the wall in a viscoelastic fluid as opposed to a Newtonian fluid, where the particle is pushed away from the wall. Therefore, the corner becomes a basin of attraction for the particle at $Wi=0.1$. According to the observed force-map, the expected focusing pattern is illustrated in Fig. 2(d). This result also rationalizes the particle behavior found in 53, where the particles aggregate at the corner and on the diagonal of the channel. The induced lateral force results in the appearance of two equilibrium points along the diagonal line among which the one closer to the channel center is stable and the other one demonstrates an unstable behavior. Furthermore, there are three unstable equilibrium points on the main axis whose locations are shown in the force profile presented in Fig. 3(b). According to this profile, the channel center and the equilibrium point near the wall are unstable due to the positive slope of the force profile. The middle equilibrium point is also identified as an unstable point despite the negative slope of the curve. This behavior is attributed to the positive lateral force along the z-direction, which pushes the particle away from the main axis and turns this point to a sub-stable equilibrium point.

Figure 2(e) illustrates the force field for $Re=5$ and $Wi=0.5$. In this case the channel can be divided into two regions by a separatix. The region closer to the channel center (indicated by the region inside the dashed curve) attracts the particle toward the centerline, whereas the outer region pushes the particle toward the corner of the microchannel. Hence, only the channel center is stable and other equilibrium points along the diagonal of the channel and the main axes are unstable. It should be noted that the corner is also a basin of attraction for the particles. The corresponding focusing pattern is plotted in Fig. 2(f). The calculated force-map indicates the reason for the particle focusing pattern observed in previous experimental and numerical studies for the cases where the inertial effect is small and the elastic force is dominant53, 17, 62, 48. Figure 3(b) also shows a negative lateral force near the central region and a negative slope for the force profile at the center that leads to stability of the channel center. Another point that should be noted is the change in the force profile near the wall region as the lateral force decreases with increasing the elasticity effect for $Wi$ above $0.5$ in Fig. 3(b). As shown in Fig. 3, the change in $Wi$ number alters the force profile significantly, leading to various focusing patterns for the particles. This effect is more visible for smaller inertial effects as illustrated in Fig. 3(a) ($Re=1$). Increasing the $Wi$ number changes the convexity of the force profile along the main axis. In the Newtonian case ($Wi=0$), the convexity of the force is negative along the entire axis, while it becomes positive for higher $Wi$ numbers. Furthermore, the channel center which is unstable for $Wi=0$, becomes stable for any $Wi$ number above $0.5$ and the corner becomes the basin of attraction. This change is observed for both $Re=1$ and $5$. According to Fig. 3(a) and (b), the location at which the force profile crosses the horizontal line in the range of $Wi>0.5$ shifts toward the channel center with increasing elasticity. Hence, the size of the separatix shrinks with $Wi$. This behavior can explain the results found in previous studies17, 48, 53 as larger fraction of particles get attracted to the corner with increasing the elastic effects ($Wi$).

In order to investigate the effect of elasticity, particularly in the case where inertia is small, we split the lateral force plotted in Fig. 2(e) into two components: (i) elastic force ($F_{el}$) and (ii) inertial force ($F_{in}$). Figure 4 shows the distribution of elastic and inertial forces for $Re=5$ and $Wi=0.5$. The elastic force (represented in Fig. 4(a)) drives the particle toward the center in the entire channel unless the particle is positioned near the wall region, where the elastic force direction reverses and the particle is pushed toward the wall. Contrary to the elastic force profile, the inertial force shown in Fig. 4(b) repels the particle from the center across the entire channel, however, the magnitude of the inertial force becomes negligible compared to that of elastic force for in the region close to the channel center. Hence, we observed a trapping region formed near the center of the microchannel and the particles that fall outside this region migrate to the corner. The comparison between elastic and inertial forces shown in Fig. 4 can explain the observed experimental results in previous studies. It should be noted that the magnitude of the elastic and inertial forces are significantly different for some cases and we use different force scale bars, which are shown in the figures.

Figure 3 shows that increasing $Wi$ number increases the magnitude of the lateral force near the channel center (for $Wi>0.1$). This effect accelerates the transverse migration of the particle across the microchannel and leads to a shorter critical length for the microchannel required for complete focusing of the suspending particles25, 48. The second normal stress is non-zero in a Giesekus fluid modeled in this paper and it causes a secondary flow across the microchannel which affects the particle migration in the channel56, 31, 48. In order to show the effect of the secondary flow, the velocity profile is illustrated in Figure 5 at the channel inlet (far from the particle) and the location of particle center at $Re=5$ and $Wi=3$.

3.2 Migration in an intermediate inertial regime

In this section, we investigate the particle focusing pattern at $Re=10$ for various Weissenberg numbers. The larger inertial effect leads to a significant change in the generated flow field compared to that of the previous section. This change can be reflected in the induced force-map shown in Fig. 6. Figure 6(a) illustrates the force field at $Wi=0.1$. The noticeable difference between this case and the Newtonian fluid is the induced lateral force near the channel face center that attracts the particle toward the wall leading to the existence of a basin of attraction at the channel face.

Furthermore, the generated lateral force turns the equilibrium point on the diagonal of the channel stable and those on the main axes and channel center unstable. This phenomenon can be observed in Fig. 7, where the force profile along the main axis is plotted, showing the locations where the lateral force becomes zero and their corresponding stability status. Hence, the expected focal pattern looks like the one illustrated in Fig. 6(b) which has not been discovered in previous studies. Increasing the $Wi$ number changes the force field significantly as displayed in Fig. 6(c). The radially directed lateral force drives the particle toward the channel corner across the entire channel. Hence, the corner becomes the only basin of attraction resulting in the focal pattern presented in Fig. 6(d). Correspondingly, the positive value of the force along the entire main axis for $Wi=0.5$ shown in Fig. 7 indicates that the particle is pushed away from the center regardless of its location in the microchannel. The calculated force field shown in Fig. 6(c) explains the particle configuration reported in previous experimental studies17, 53. For $Wi=3$, the lateral force has a distribution similar to that of a Newtonian fluid. In this case (shown in Fig. 6(e)) the radially directed force pushes the particle away from the wall and channel center and creates an annulus ring which is similar to the curve marked in Fig. 2(a). Consequently, the predicted focal pattern of the particle shown in Fig. 6(f) is similar to that of a Newtonian fluid.

According to Fig. 7, the fluid elasticity significantly affects the magnitude of the lateral force such that the force profile near the wall region markedly changes and its direction reverses for high $Wi$ numbers, leading to a stable equilibrium point on the main axes. On the other hand, the slope of the force profile is positive at the center for the entire range of Weissenberg number explored in this study, implying that the channel center is not a stable position for the particle in this range of parameters. In order to investigate the effect of elasticity on the total lateral force the force-map of inertial and elastic forces are plotted for $Wi=0.1$ to $3$ in Fig. 8.

The results show that the inertial force dominates the elastic force. The direction of total force matches that of the inertial force, while the elastic force has a similar profile across the entire channel for various values of $Wi$ number. This behavior indicates that changing the elasticity of the fluid changes the velocity field in the channel. Consequently, the inertial force is affected and the resulting focal pattern alters, while the elastic force itself remains relatively unchanged. In other word, the fluid elasticity affects the particle dynamic indirectly by changing the flow field and not directly by means of elastic force.

3.3 Migration in a high inertial regime

With increasing the Reynolds number, the particle dynamics changes significantly from those observed in previous sections. Figure. 9(a) shows the force profile for $Re=30$. The force has a similar trend along the channel main axis for the entire range of $Wi$ number. The force has a positive value near the center and it becomes negative in the near wall region, i.e., the particle is pushed away from the center and the wall, which leads to the existence of an off-center equilibrium point along the main axis. As opposed to previous sections, where the elasticity changes the force profile significantly, the force distribution across the microchannel remains unchanged for the entire range of $Wi$ number studied in this work.

The slope of the force profile at the center and the off-center equilibrium points indicates that the particle is unstable and stable at these points, respectively. According to the force-map shown in Fig. 9(b), there is an unstable equilibrium point along the diagonal of the channel at $Re=30$ and $Wi=0.5$ that is similar to that of a Newtonian fluid ($Wi=0$). Hence, the behavior discovered in the numerical and experimental studies by 53 and 56 can be rationalized using our computational results. As shown in Fig. 9(a), the force gets stronger along most of the main axis as the $Wi$ number increases (for $Wi>0.1$). Hence, the particles is pushed away from the center and the wall faster, increasing the transverse migration toward the annulus ring and leading to a smaller critical length of microchannel required for particle focusing. However, this trend is not observed for $Wi=0.1$ as the force magnitude is smaller than that of a Newtonian fluid.

Figure 10 illustrates the distance of the off-center equilibrium points from the channel center and their stability for the entire range of $Re$ and $Wi$ numbers studied in this work along the main axis and the diagonal of the channel.

As shown in Fig. 10(a), the equilibrium points for $Wi=0$ are all stable on the main axis and their location does not change with $Re$, which is in agreement with previous results found in Li et al. 56. However, the results show that the equilibrium points for $Wi\neq 0$ are mostly unstable in the low inertial regime ($Re=1,5$). The distance of these points are smaller than that of a Newtonian fluid and they approach the channel center with increasing $Wi$ number (for $Wi>0.5$), indicating reduction of trapping area with increasing the elastic effects. Oppositely, the equilibrium points for $Wi\neq 0$ are all stable on the main axis in the high inertial regime ($Re=20,30,50$). In this range of $Re$ number, the equilibrium points at $Wi=0.1$ have relatively the same distance from the channel center as that of a Newtonian fluid which is due to the small effect of elasticity compared to inertial effect, while for a higher $Wi$ number this distance is larger than that in a Newtonian fluid. The results also show that the difference between the location of equilibrium points becomes smaller as the $Re$ number increases. This can be attributed to the dominant effect of the flow inertia compared to the elastic effects. In contrary to the stability of the equilibrium points on the main axis, Fig. 10(b) indicates that most of equilibrium points on the diagonal of the channel are unstable. The results show that the equilibrium points approach the channel center with increasing $Wi$ number for low inertial effects, while this behavior changes for larger inertia such that the distance between equilibrium points and channel center initially increases and subsequently the equilibrium point approaches the center with increasing the $Wi$ number.

The location of equilibrium points along the main axis and the diagonal of the channel does not change with the $Re$ number in a Newtonian fluid. In order to explain this phenomenon, the force profile along the channel main axis is plotted for $Wi=0$ in Fig. 11(a). According to this figure, the change in the $Re$ number changes the magnitude of the lateral force acting on the particle. As the inertial effect becomes stronger the lateral force magnitude increases significantly, however, the location at which the lateral force becomes zero is the same for the entire range of $Re$ number studied in this work. This explains the results reported by 56 In order to investigate the combined effects of inertia and elasticity on the stability of the channel centerline equilibrium point, the stability phase diagram is illustrated in Fig. 11(b). As shown in this figure, the channel center is a stable equilibrium point when the elastic effect dominates the inertial effect, while with increasing the $Re$ number the focal pattern looks similar to that of a Newtonian fluid, where the particle is not expected to travel toward the channel center.

4 Conclusions

In this work, we conduct 3D numerical simulations to find the distribution of lift force acting on the particle in viscoelastic fluids. As a result, we predict the location of equilibrium points and their corresponding stability for a wide range of parameter space which is important for designing microfluidic devices relying on viscoelastic effects. The results for low inertial regime show that the force field acting on the particle changes significantly with increasing the elasticity of the fluid. In the case of a Newtonian fluid the particle is pushed away from the wall and the channel center to reach an annulus ring across the channel. Subsequently, the particle moves toward the equilibrium points on the main axis. Our results show that only off-center equilibrium points along the main axes are stable and other equilibrium points are unstable in a Newtonian fluid. Increasing the elasticity changes the induced force field such that the main axis has no stable equilibrium point and the particle focuses on the corner and off-center points along the diagonal of the channel. Further increase in the Weissenberg number leads to dominance of the elastic force over the inertial lift force, shifting the off-center equilibrium points on the diagonal of the channel toward the channel center.

In the intermediate inertial regime, we also observe various focal patterns by changing the fluid elasticity. In the low $Wi$ number range, the particle focuses at the wall face center and off-center points along the diagonal of the channel. This configuration changes for larger $Wi$ numbers, where the particle is driven radially toward the corner in the entire channel cross-section, leading the corners to be the only basin of attraction in the microchannel. Further increase in the elasticity results in a configuration similar to that of a Newtonian fluid, in which the particle aggregates only at an off-center point on the main axis and equilibrium points on the center and diagonal of the channel are unstable. By splitting the total lift force experienced by the particle into inertial and elastic lift forces we conclude that the elastic effect modifies the velocity field in the microchannel. As a result, the inertial lift force changes accordingly, leading to various particle configurations. However, the direction and the magnitude of the elastic force remain relatively unchanged in the intermediate inertial regime implying that the fluid elasticity affects the particle dynamics indirectly by changing the velocity field. Our results for the high inertial flow indicate that the force profiles are similar for the entire range of $Wi$ number studied in this work. Accordingly the particle configuration is also similar to that of a Newtonian fluid, in which the particles aggregate only at the off-center point along the main axis.

5 Acknowledgement

This research was partially supported by a Grant from National Science Foundation [CBET-1705371]. This work used the Extreme Science and Engineering Discovery Environment (XSEDE) 63, which is supported by the National Science Foundation grant number ACI-1548562 through allocation TG-CTS180066.