Next Article in Journal
Microfluidic Device for Droplet Pairing by Combining Droplet Railing and Floating Trap Arrays
Next Article in Special Issue
Visualization Experimental Study on Silicon-Based Ultra-Thin Loop Heat Pipe Using Deionized Water as Working Fluid
Previous Article in Journal
Controllable Fast and Slow Light in Photonic-Molecule Optomechanics with Phonon Pump
Previous Article in Special Issue
Formation and Elimination of Satellite Droplets during Monodisperse Droplet Generation by Using Piezoelectric Method
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Neutrally Buoyant Particle Migration in Poiseuille Flow Driven by Pulsatile Velocity

School of Mechanical Engineering, Hangzhou Dianzi University, Hangzhou 310018, China
*
Author to whom correspondence should be addressed.
Micromachines 2021, 12(9), 1075; https://doi.org/10.3390/mi12091075
Submission received: 4 August 2021 / Revised: 1 September 2021 / Accepted: 3 September 2021 / Published: 6 September 2021
(This article belongs to the Special Issue Advances in Heat and Mass Transfer in Micro/Nano Systems)

Abstract

:
A neutrally buoyant circular particle migration in two-dimensional (2D) Poiseuille channel flow driven by pulsatile velocity is numerical studied by using immersed boundary-lattice Boltzmann method (IB-LBM). The effects of Reynolds number ( 25 R e 200 ) and blockage ratio ( 0.15 k 0.40 ) on particle migration driven by pulsatile and non-pulsatile velocity are all numerically investigated for comparison. The results show that, different from non-pulsatile cases, the particle will migrate back to channel centerline with underdamped oscillation during the time period with zero-velocity in pulsatile cases. The maximum lateral travel distance of the particle in one cycle of periodic motion will increase with increasing R e , while k has little impact. The quasi frequency of such oscillation has almost no business with R e and k . Moreover, R e plays an essential role in the damping ratio. Pulsatile flow field is ubiquitous in aorta and other arteries. This article is conducive to understanding nanoparticle migration in those arteries.

1. Introduction

Particle two-phase flow is a very complex problem, which ubiquitously exists in nature, industry, hemodynamics, such as the formation and movement of sand dunes, haze (PM2.5), ventilation dusting system, spread of virus (COVID-19), inertial microfluidics, drug delivery in blood, etc. Numerous researches have revealed the behavior of the particles in fluid flow depends on Reynolds number ( R e = U m H / ν , where U m is the maximum inlet velocity, H is the channel width, ν is the fluid kinematic viscosity) and blockage ratio ( k = D p / H , donates the ratio of particle diameter D p and the channel width), whether or not the particles are neutrally buoyant [1,2,3,4,5,6]. The density of the neutrally buoyant particles is the same as the suspension fluid, which means the particles will suspend in the fluid.
Segré and Silberberg [7] first discovered experimentally neutrally buoyant spherical particles would migrate to a radial equilibrium position in a pipe flow and form the Segré and Silberberg (SS) annulus, which is known as SS effect. This phenomenon prompted a lot of correlation research to reveal the underlying mechanism. Ho and Leal [8] theoretically studied particle migration in two-dimensional (2D) Poiseuille flow. Asmolov [9] interpreted that the particles migration was due to the effect of inertial lift by using the matched asymptotic expansions method. The results indicated the wall induced inertial lift became significant in the thin layers near the channel wall, and such lift could be neglected when the particles are far away from the wall. Matas et al. [10] found that the particles would move closer to the circular tube wall as R e increased and revealed additional inner annulus when R e was greater than 600. Moreover, if R e exceeded 700, the particles in the inner annulus accounted for the majority. Matas et al. [11] also utilized the matched asymptotic expansions method to calculate the lateral force in the pipe geometry (used to be in the plane geometry), but they did not find the second zero lateral force intersection point which indicates the inner annulus. Thus, they concluded the inner annulus was most likely due to finite-size effect. Hood et al. [12] calculated the lateral forces by a perturbation analysis. Morita et al. [13] predicted that all particles in the inner annulus would return to the SS annulus according to their experimental results when R e is less than 1000 and the tube is long enough. Due to their equipment limitations, the tube length was only 500 times the particle diameter. Then, Nakayama et al. [14] increased the length of the tube to 1000 times of the particle diameter, and the results showed that three regimes were eventually formed: only the SS annulus, only the inner annulus or those two annuli exist at the same time. The transition between these three regimes was determined by the critical R e which decreases with increasing of k .
Besides the aforementioned theoretical and experimental methods, computational fluid dynamics (CFD) simulation has become a powerful tool in analyzing particle-fluid interaction. Feng et al. [15] adopted finite-element method to investigate a circular particle migration in Poiseuille flow. Their simulations agreed qualitatively with the results of perturbation theories and pertinent experiments. By using LBM [16,17,18,19], the same problem was studied, and the SS effect was reconstructed. Shao et al. [20] found the inner annulus for elevated R e by using the fictitious domain method. Abbas et al. [21] mentioned that the equilibrium position (TEP) depends exclusively on R e and k . Recently, inertial microfluidics can precisely separate particles with or without extra external force field by realizing SS effect [22,23,24,25,26,27,28].
All the research works mentioned above are based on a non-pulsatile flow field, but in the artery, the blood pumped by the heart behaves as a pulsatile flow field. To the best of our knowledge, there are no articles related to particle migration in the pulsatile flow field.
Cancer is one of the leading causes of death in the world; nanoparticles are widely used for cancer therapy. Ideally, the therapeutic nanoparticles system should be able to deliver drugs just to the tumor and have no severe side effects on the body [29]. However, nanoparticle movement in non-pulsatile blood flow field is quite different from that in the pulsatile blood flow field.
In this study, we perform a series of CFD simulations to investigate the migration of one neutrally buoyant circular particle in 2D Poiseuille flow driven by pulsatile velocity. The influence of R e and k on the particle migration is analyzed in detail. The difference between particle migration driven in pulsatile and non-pulsatile velocity is illustrated. Engineering precision therapeutic nanoparticles system in artery [30] can be attainable by understanding the particle migration in pulsatile blood flow field. Furthermore, it can help to optimize therapeutic nanoparticles system for achieving precise medical against cancer near arteries.
The organization of this article is as follows. The numerical method, boundary conditions, and problem description are introduced in Section 2. In Section 3, we simulated TEPs of one particle at several specific parameters, and our computational code is validated by comparing with the published paper. Afterwards, the simulation results are presented and discussed in Section 4. Finally, conclusions are provided in Section 5.

2. Method and Problem

LBM is widely used to simulate particle migration, turbulence, flow in porous media, multiphase flow, non-Newtonian rheology, and so on [31,32,33,34,35,36,37,38,39,40,41,42,43]. Because of its advantages in efficiency, easy to code, and parallel run, LBM has become a very popular CFD numerical tool.

2.1. Lattice Boltzmann Method

In this work, the single-relaxation time (SRT) lattice Bhatnagar–Gross–Krook (LBGK) Boltzmann method is adopted to solve particles migration in incompressible viscous flow [44]:
f i ( x + e i Δ t , t + Δ t ) = f i ( x , t ) 1 τ [ f i ( x , t ) f i ( e q ) ( x , t ) ] ,
where f i ( x , t ) is the distribution function at space coordinate x = ( X , Y ) and time t in the i th direction; f i ( e q ) ( x , t ) is the correspond equilibrium distribution function; τ is the SRT; Δ t is the time step; the discrete velocities e i of 2D nine-velocity (D2Q9) model are shown in Figure 1; c = Δ x / Δ t is the lattice velocity; Δ x is the lattice spacing. For coding conveniently, both the time step and the lattice spacing are set to be equal to 1, which result in Δ t = Δ x = c = 1 .
The equilibrium distribution function can be calculated by [44]:
f i ( e q ) ( x , t ) = ω i ρ f [ 1 + 3 e i · u c 2 + 4.5 ( e i · u ) 2 c 4 1.5 u 2 c 2 ] ,
where ω i is the weight factor with ω 0 = 4 / 9 , ω 1 ~ 4 = 1 / 9 and ω 5 ~ 8 = 1 / 36 , ρ f is the fluid density and u is fluid velocity which can be determined by:
ρ f =   f i ( x , t ) ,   u = 1 ρ f   e i f i ( x , t ) .
For low Mach number, the Navier–Stokes equations can be derived from the lattice Boltzmann equation by utilizing Chapman–Enskog expansion [45].

2.2. Improved Bounce-Back Scheme

In the LBM simulations, improved bounce-back scheme is of particular importance, and it allows to implement no-slip boundary condition on the surface of moving particle [46], which will be explained briefly below.
As shown in Figure 2, the sky-blue circles represent the fluid nodes, and the red thin diamonds donate the solid nodes. The orange triangles indicate uncovered fluid nodes means that those nodes are located inside the particle at time t and will locate outside the particle when the particle travels after one lattice time step. Similarly, the purple squares donate covered fluid nodes represent that those nodes will change from fluid nodes to solid nodes after the particle moves.
Take fluid node A as an example, after the streaming step, three unknown distribution functions ( f 3 , f 7 and f 4 denoted by three blue arrows shown in Figure 2) need to be determined by applying improved bounce-back scheme. Therefore, the distribution function f 7 can be calculated by:
f 7 ( A ) = { q ( 1 + 2 q ) f 5 ( S ) + ( 1 4 q 2 ) f 5 ( A ) q ( 1 2 q ) f 5 ( B ) 2 ω 5 ρ f e 5 · u D c s 2 ,   q < 1 2 , 1 q ( 2 q + 1 ) f 5 ( S ) + 2 q 1 q f 7 ( B ) 2 q 1 2 q + 1 f 7 ( C ) 2 ω 5 ρ f q ( 2 q + 1 ) e 5 · u D c s 2 ,   q 1 2 ,
where S is the nearest solid node along e 5 direction; B and C are the two nearest fluid nodes along e 7 direction; the blue star D denotes the boundary location which is the intersection point of particle boundary and line A S , meanwhile, q can be determined by q = | A D | / | A S | ; u D is the velocity of the boundary location D ; c s = c / 3 is the speed of sound. The other unknown distribution functions of fluid nodes around particle boundary can be solved in the similar method.

2.3. Force, Torque, and Particle Motion

Let us continue to take the boundary location D as an example. As shown in Figure 2, the hydrodynamic force and torque acting on the solid particle migrating in fluid can be integrated by adopting momentum exchange algorithm [46,47] as follows:
F ( h ) ( x + q e 5 , t ) = e 5 [ f 5 ( x + e 5 , t ) + f 7 ( x , t ) ] , T ( h ) ( x + q e 5 , t ) = ( x + q e 5 x p ) × F ( h ) ( x , t ) ,
where x p is the position of the particle.
When the particle is moving in the lattice grid from time t to t + t , as shown in Figure 2, the additional force and torque due to uncovered fluid node and covered fluid node exerted on the particle can be computed by [48]:
F ( c ) ( x , t ) = ρ f ( x , t ) u ( x , t ) , T ( c ) ( x , t ) = ( x x p ) × F ( c ) ( x , t ) , F ( u ) ( x , t ) = ρ f ( x , t ) u ( x , t ) , T ( u ) ( x , t ) = ( x x p ) × F ( u ) ( x , t ) .
Moreover, in order to avoid unphysical overlapping between the particle and the channel wall, the extra lubrication force model is needed [49]:
F ( l ) = { 0 , h h c , 1.5 π ρ f ν [ D p ( 1 h 1 h c ) ] 1.5 U p , h < h c ,
where ν = c s 2 ( τ 0.5 ) Δ t is the kinematic viscosity of the fluid; D p is the diameter of the particle; U p is the particle velocity towards the wall; h is the minimum gap between the particle and the wall; h c = 1.5 Δ x is the cutoff distance whether to consider the lubrication force or not.
By summation of the forces and torques acting on the particle, the movement of the particle can be solved explicitly using Newton’s second law:
a p t + Δ t = (   F ( h ) +   F ( c ) +   F ( u ) + F ( l ) ) / m p , u p t + Δ t = u p t + 0.5 ( a p t + Δ t + a p t ) Δ t , x p t + Δ t = x p t + 0.5 ( u p t + Δ t + u p t ) Δ t , α p t + Δ t = (   T ( h ) +   T ( c ) +   T ( u ) ) / I p , w p t + Δ t = w p t + 0.5 ( α p t + Δ t + α p t ) Δ t , θ p t + Δ t = θ p t + 0.5 ( w p t + Δ t + w p t ) Δ t ,
where a p , u p , α p , w p , θ p , m p , and I p are the translational acceleration, velocity, rotational acceleration, rotational velocity, angle, mass and moment inertia of the particle.

2.4. Problem

The configuration of one circular particle migrating in Poiseuille flow is shown in Figure 3. A parabolic velocity profile with the maximum velocity U m is set at the left inlet boundary in the positive X direction, and the velocity in Y direction is zero. For U m , there are two cases: non-pulsatile or pulsatile. If pulsatile, U m will change over time, otherwise, it will be a constant. For considering reproducibility of this work, patient-specific velocity profile [50] will not be adopted to impose at the inlet boundary. The half-period of the sinusoidal function at time interval [ t 0 ,   t 1 ] (systolic period) and zero at [ t 1 ,   t 2 ] (diastolic period) is utilized which can be seen from Figure 3. In the systolic period, t 1 t 0 = 0.3   s , and in the diastolic period, t 2 t 1 = 0.5   s , which means one cardiac cycle lasts 0.8   s . The unit conversion factor from lattice time to physical time is 10 5 . At the right boundary, the normal derivative of the velocity is zero and the pressure is set to be p o u t = ρ f c s 2   [19]. No-slip boundary conditions are imposed at the top and bottom channel walls.
If there are no additional statements, the channel length L is 500 Δ x and the height H is 100 Δ x . TEP is basically unchanged for different channel length ( 300 Δ x , 400 Δ x , 500 Δ x , 600 Δ x , 700 Δ x and 2000 Δ x ) . So L = 500 Δ x is adopted same as Wen et al. [17]. As shown in Figure 3, the circular particle center is located at [ X s ,   Y s ] initially. The diameter of the circular particle in this work is 15 Δ x , 20 Δ x , 25 Δ x , 30 Δ x , 35 Δ x , and 40 Δ x , which means k = D p / H = 0.15 ,   0.20 ,   0.25 ,   0.30 ,   0.35 ,   0.40 , respectively. Meanwhile, R e = U m H / ν = 25 ,   50 ,   75 ,   100 ,   150 ,   200 is studied by changing the kinematic viscosity of the fluid.
Driven by the fluid velocity, the particle will always travel along X axis in positive direction, so an infinite channel is needed to avoid the particle moving out of the simulation domain which can be achieved by moving domain technique [6,51,52]. When the X -coordinate of the particle exceeds X s + Δ x , the fluid field and the particle need to shift on lattice spacing left, which ensures that the particle will never travel too far from its original position. Meanwhile, it should be noted that the actual moving distance along X axis in positive direction should add up one lattice spacing when performing this shift once.

3. Validation

To validate the accuracy of our LBM code, two benchmark cases are implemented below.
In the first case, R e is set to be 50, and the terminal particle Reynolds number R e p = U x , p D p / ν is about 9.63, where U x , p is the terminal particle velocity in X direction. It is very close to R e p while the particle is driven by pressure difference [17]. Our simulating results with k = 0.25 and k = 0.35 are plotted in Figure 4a,b, respectively. The SS effect is quite obviously found and TEP of the particle has nothing to do with the initial horizontal position ( Y s / H = 0.20 ,   0.25 ,   0.35 ,   0.40 ,   0.45 ) which only changes the trajectory from initial position to TEP. All the results, even the curve shape shown in Figure 4a,b are consistent with those of Wen et al. [17].
The foregoing results are validated only when R e is 50. The second case is adopted to validate over the entire range of R e = 20 ,   40 ,   100 ,   200 . The particle diameter is D p = 22 , the channel width is H = 200 , and the channel length is L = 1000 , which leads to the blockage ratio being k = 0.11 . Figure 5 shows the comparison of the present results with previous ones simulated by Di Chen et al. [53] The comparison shows TEPs are in good agreement and the particle will be closer to channel centerline with increasing R e in 2D Poiseuille flow.

4. Results and Discussions

4.1. Fluid and Particle Interaction

After validation, we first simulate the particle migrating in non-pulsatile flow at k = 0.25 and R e = 50 . Figure 6a–c, shows the contour view of the dimensionless pressure p = ( p p o u t ) / p o u t , the dimensionless fluid velocity in X direction U x = U x / U m and in Y direction U y = U y / U m , respectively. Due to the moving domain method adopted here, the particle will always locate around the middle of the simulation domain which is significantly disturbed by the present of the particle. Obviously, the pressure strip can be observed upstream and downstream of the particle. The pressure at the upper left and down right side of the particle is higher than the upper right and down left corner, and hence generates a particle rotation in clockwise direction illustrated by the black arrow in Figure 6a. Due to non-slip boundary condition at the boundary location of particle, clockwise rotation of the particle will induce fluid flowing upward at left and downward at right as shown in Figure 6c. Moreover, it can be seen from Figure 6b that the particle follows with fluid movement quite well [54].

4.2. Trajectory

Several dominant forces acting on the particle drive it to TEP in Y direction. The first force is shear gradient lift force due to the parabolic velocity profile, which points to the wall. The second force is wall induced lift force due to the interaction between the particle and the channel wall calculated by using added lubrication force model introduced before which directs towards the channel centerline. The third force is called Magnus force [55] due to the rotation of the particle migrating in fluid. As demonstrated above, the particle in non-pulsatile flow travels in X direction and at the same time it rolls in clockwise direction. Thus, the Magnus force is directed to the channel centerline. Certainly, when the particle moves in fluid, it will be affected by the drag force. In summary, there are at least four forces governing particle migration in fluid suspension.
Figure 7 shows the particle center trajectory along the channel for k = 0.25 and k = 0.35 at R e = 50 . Obviously, several cardiac cycles later, the particle driven by pulsatile velocity migrates around TEP in non-pulsatile flow periodically. In the systolic period, the particle laterally migrates towards wall mainly affected by shear gradient lift force. While in the diastolic period, the shear gradient lift force disappears, and the particle still rotates because of inertia, so it will migrate back to the channel centerline under the Magnus force. The moving direction is illustrated by black arrows which are shown in the subplots of Figure 7. The particle oscillates in spiral-shaped structures and will travel towards the wall again in another systolic period.
Figure 8a,c gives the particle center trajectory versus time ( t = t / 10 5 ) variation with R e at k = 0.25 , and ( b , d ) donates different k at R e = 50 . Moreover, ( c , d ) is the enlarged view of ( a , b ) in one stable cardiac cycle ( t = [ 8.8 ,   9.6 ] ) illustrated by black dash line box, respectively. Overall, TEPs are all closer to the channel centerline when R e or k increases, whether or not in pulsatile flow. Consequently, the particle will take longer (or more cardiac cycles) to reach TEP. As shown in Figure 8c, we define as the dimensionless distance between the highest and lowest position in this stable cardiac cycle, δ is the signed dimensionless distance between TEP in non-pulsatile flow and the lowest position (negative value means the particle did not exceed TEP in non-pulsatile flow), and t is time in this cardiac cycle when the particle locates the lowest position. In general, the influence of R e on δ , , t is greater than k . As R e increases, the viscosity of the fluid decreases, thus the drag force acts on the particle decreases. The particle, therefore, can laterally migrate farther in systolic period, even exceeds TEP in non-pulsatile flow. Meanwhile, the particle will take longer to migrate from the highest position to the lowest position for its longer travel distance. As a result, Figure 8e shows that δ ,   and t all increase monotonously with increasing R e . The effect of k on δ , , t is a little more complicated. As k increases, the particle inertia increases, and the drag force also increases. Shear gradient force increases with increasing k , but decreases when the particle is closer to channel center. Consequently, it can be seen from Figure 8f, δ ,   increases at first and then decreases when k increases, the maximum of δ ,   occurs at k = 0.2 . Moreover, t was mostly unchanged with increasing k .

4.3. Orientation

Observed in Figure 9, the particle always rotates clockwise while migrating in the channel, irrespective of R e , k , and whether or not in pulsatile flow. In the diastolic period, the particle rotates much slower but still in the clockwise direction. If the particle is closer to the channel wall, it will experience larger gradient of fluid velocity. As a result, the smaller Re or k is, the faster the particle rotates. Figure 8 shows that, when R e < 100 , TEP of the particle is very similar, so they rotate almost at the same speed ( w ) which is shown in the subplot of Figure 9a. In addition, it can be seen in the subplot of Figure 9b, the particle rotate speed is almost linearly coherent with k . Furthermore, the moment inertia of the particle is proportional to the square of the particle diameter, i.e., k . Consequently, by comparing Figure 9a,b, the influence of k on w is much larger than R e .

4.4. Damping

As mentioned before, the particle oscillates in the diastolic period which is shown in Figure 7. We choose the particle velocity which indicates the interaction between the particle and fluid [56] for analysis. Moreover, the particle velocity in X direction ( U x ) is preferred, because the particle velocity curve is not center symmetry in Y direction.
For better illustration, the particle velocity in X direction is normalized by: U x = U x / U m . The curves of U x in Figure 10 are much like a spring-mass system which is underdamped. Figure 10a shows, when R e is small, U x damps out rapidly after several quasi periods. If R e is small enough, like R e 1 , over damped is expected. Obviously, in Figure 10a–f, the damping effect becomes weaker with increasing R e . We fit the upside and downside envelope curve by using exponential function for purpose. For example, A u p = e 15.99 t + 142.39 in Figure 10a, and the decay rate λ = 15.99 . The constants of fitting exponential function of the upside and downside are almost the same for each R e . And the absolute values of those constants decrease with increasing R e . The influence of k is also analyzed, but it can be seen from Figure 11a–f that k has almost no effect on damping.
Figure 12a gives the quasi frequency ( f q ) of U x oscillation at different R e and k . It can be seen from the subplot that the quasi frequency increases as increasing k , but the impact of k on the quasi frequency is relatively small. However, the quasi frequency will not change when R e increases except for R e = 25 . As mentioned earlier, in Figure 10a, when R e = 25 , U x damps out rapidly after several quasi periods. The first few quasi periods are longer than the last ones, which will result in smaller quasi frequency when R e = 25 .
As shown in Figure 12b, the damping ratio ζ λ / ( 2 π f q ) , in our cases, is always smaller than 1, which determines this is an underdamped system, and this system will die out slower when ζ decreases. Additionally, ζ of U x oscillation basically does not change with variation of k , but it decreases with increasing R e , i.e., R e is the decisive parameter on ζ [57]. A similar conclusion is made by C.A. Coulomb in 1784 by using several material plates hanging by a metal wire with an initial torsion angle and then released to start timing until the plate oscillation until dying out. ζ is only related to the viscosity of fluid and not the material of plate. The reason for damping is the friction inside the fluid which is called Newton’s law of friction, not the interaction friction between the plate and the fluid. Finally, in Figure 12b, the relation between ζ and R e is fitted as: ζ = 0.83 R e 0.64 .

5. Conclusions

IB-LBM was used to simulate one neutrally buoyant circular particle migration in 2D Poiseuille channel flow driven by pulsatile and non-pulsatile velocity. The moving domain technique is adopted to achieve that the particle can migrate in infinite channel. The results show that the particle moving in pulsatile flow is slightly different from that in non-pulsatile one. It will laterally migrate back to the channel centerline with small oscillations in the spiral-shaped structure during the diastole and move back toward TEP during the systole. The effect of R e on ζ is decisive. This research may shed some light on understanding the particle behavior in Poiseuille flow driven by pulsatile velocity for optimizing therapeutic nanoparticles system in arteries.
The limitation of this work is that R e varies only from 25 to 200. Smaller R e leads to bigger SRT, thus the simulation will not be accurate. In contrast, bigger R e results in smaller SRT which the simulation will diverge. Furthermore, k range only from 0.15 to 0.40. Because the channel height is set to be 100 lattice units, if k is smaller than 0.15, the simulation resolution will not be adequate. On the other hand, if k is bigger than 0.4, half channel width will be blocked by the particle. Parameters need to choose carefully to get a more varied range of R e and k . This will be a future direction of this work. Another future direction will be considering more particles or simulating particle migration in three-dimensional pipe pulsatile flow.

Author Contributions

Conceptualization, L.H. and Z.Z.; methodology, L.H.; software, L.H.; validation, L.H. and J.D.; writing—original draft preparation, L.H.; writing—review and editing, J.D. and Z.Z.; funding acquisition, J.D. and Z.Z. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Fundamental Research Funds for the Provincial Universities of Zhejiang (GK199900299012-026) and National Natural Science Foundation of China (Grant Nr. 11572107).

Data Availability Statement

The data that support the findings of this study are available from the first author upon reasonable request.

Conflicts of Interest

The authors report that they have no conflict of interest.

References

  1. Hu, X.; Lin, J.Z.; Chen, D.M.; Ku, X.K. Influence of non-Newtonian power law rheology on inertial migration of particles in channel flow. Biomicrofluidics 2020, 14, 14105. [Google Scholar] [CrossRef]
  2. Hu, X.; Lin, J.Z.; Chen, D.M.; Ku, X.K. Stability condition of self-organizing staggered particle trains in channel flow. Microfluid Nanofluidics 2020, 24, 25. [Google Scholar] [CrossRef]
  3. Jiang, R.J.; Lin, J.Z.; Tu, C.X. Poiseuille flow-induced vibrations of a cylinder in subcritical conditions. J. Fluids Struct. 2018, 82, 272–286. [Google Scholar] [CrossRef]
  4. Li, G.J.; McKinley, G.H.; Ardekani, A.M. Dynamics of particle migration in channel flow of viscoelastic fluids. J. Fluid Mech. 2015, 785, 486–505. [Google Scholar] [CrossRef] [Green Version]
  5. Yu, Z.S.; Wang, P.; Lin, J.Z.; Hu, H.H. Equilibrium positions of the elasto-inertial particle migration in rectangular channel flow of Oldroyd-B viscoelastic fluids. J. Fluid Mech. 2019, 868, 316–340. [Google Scholar] [CrossRef]
  6. Zhang, J.; Yuan, D.; Zhao, Q.B.; Teo, A.J.T.; Yan, S.; Ooi, C.H.; Li, W.H.; Nguyen, N.T. Fundamentals of Differential Particle Inertial Focusing in Symmetric Sinusoidal Microchannels. Anal. Chem. 2019, 91, 4077–4084. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  7. Segré, G.; Silberberg, A. Behaviour of macroscopic rigid spheres in Poiseuille flow Part 2. Experimental results and interpretation. J. Fluid Mech. 1962, 14, 136–157. [Google Scholar] [CrossRef]
  8. Ho, B.P.; Leal, L.G. Migration of rigid spheres in a two-dimensional unidirectional shear flow of a second-order fluid. J. Fluid Mech. 1976, 76, 783–799. [Google Scholar] [CrossRef] [Green Version]
  9. Asmolov, E.S. The inertial lift on a spherical particle in a plane Poiseuille flow at large channel Reynolds number. J. Fluid Mech. 1999, 381, 63–87. [Google Scholar] [CrossRef]
  10. Matas, J.P.; Morris, J.F.; Guazzelli, É. Inertial migration of rigid spherical particles in Poiseuille flow. J. Fluid Mech. 2004, 515, 171–195. [Google Scholar] [CrossRef] [Green Version]
  11. Matas, J.P.; Morris, J.F.; Guazzelli, É. Lateral force on a rigid sphere in large-inertia laminar pipe flow. J. Fluid Mech. 2009, 621, 59–67. [Google Scholar] [CrossRef]
  12. Hood, K.; Lee, S.; Roper, M. Inertial migration of a rigid sphere in three-dimensional Poiseuille flow. J. Fluid Mech. 2015, 765, 452–479. [Google Scholar] [CrossRef] [Green Version]
  13. Morita, Y.; Itano, T.; Sugihara Seki, M. Equilibrium radial positions of neutrally buoyant spherical particles over the circular cross-section in Poiseuille flow. J. Fluid Mech. 2017, 813, 750–767. [Google Scholar] [CrossRef]
  14. Nakayama, S.; Yamashita, H.; Yabu, T.; Itano, T.; Sugihara Seki, M. Three regimes of inertial focusing for spherical particles suspended in circular tube flows. J. Fluid Mech. 2019, 871, 952–969. [Google Scholar] [CrossRef]
  15. Feng, J.; Hu, H.H.; Joseph, D.D. Direct simulation of initial value problems for the motion of solid bodies in a Newtonian fluid. Part 2. Couette and Poiseuille flows. J. Fluid Mech. 1994, 277, 271–301. [Google Scholar] [CrossRef] [Green Version]
  16. Li, H.B.; Lu, X.Y.; Fang, H.P.; Qian, Y.H. Force evaluations in lattice Boltzmann simulations with moving boundaries in two dimensions. Phys. Rev. E 2004, 70, 26701. [Google Scholar] [CrossRef]
  17. Wen, B.H.; Li, H.B.; Zhang, C.Y.; Fang, H.P. Lattice-type-dependent momentum-exchange method for moving boundaries. Phys. Rev. E 2012, 85, 16704. [Google Scholar] [CrossRef]
  18. Tao, S.; Hu, J.J.; Guo, Z.L. An investigation on momentum exchange methods and refilling algorithms for lattice Boltzmann simulation of particulate flows. Comput. Fluids 2016, 133, 1–14. [Google Scholar] [CrossRef]
  19. Inamuro, T.; Maeba, K.; Ogino, F. Flow between parallel walls containing the lines of neutrally buoyant circular cylinders. Int. J. Multiph. Flow 2000, 26, 1981–2004. [Google Scholar] [CrossRef]
  20. Shao, X.M.; Yu, Z.S.; Sun, B. Inertial migration of spherical particles in circular Poiseuille flow at moderately high Reynolds numbers. Phys. Fluids 2008, 20, 103307. [Google Scholar] [CrossRef] [Green Version]
  21. Abbas, M.; Magaud, P.; Gao, Y.; Geoffroy, S. Migration of finite sized particles in a laminar square channel flow from low to high Reynolds numbers. Phys. Fluids 2014, 26, 123301. [Google Scholar] [CrossRef] [Green Version]
  22. Wang, R.J.; Sun, S.S.; Wang, W.; Zhu, Z.F. Investigation on the thermophoretic sorting for submicroparticles in a sorter with expansion-contraction microchannel. Int. J. Heat Mass Transf. 2019, 133, 912–919. [Google Scholar] [CrossRef]
  23. Wang, R.J.; Du, J.Y.; Guo, W.C.; Zhu, Z.F. Investigation on the Thermophoresis-Coupled Inertial Sorting of Submicrometer Particles in a Microchannel. Nanoscale Microscale Thermophys. Eng. 2016, 20, 51–65. [Google Scholar] [CrossRef]
  24. Du, J.Y.; Li, L.; Zhuo, Q.Y.; Wang, R.J.; Zhu, Z.F. Investigation on Inertial Sorter Coupled with Magnetophoretic Effect for Nonmagnetic Microparticles. Micromachines 2020, 11, 566. [Google Scholar] [CrossRef] [PubMed]
  25. Zhao, Q.B.; Yuan, D.; Zhang, J.; Li, W.H. A Review of Secondary Flow in Inertial Microfluidics. Micromachines 2020, 11, 461. [Google Scholar] [CrossRef] [PubMed]
  26. Razavi Bazaz, S.; Mashhadian, A.; Ehsani, A.; Saha, S.C.; Krüger, T.; Ebrahimi Warkiani, M. Computational inertial microfluidics: A review. Lab Chip 2020, 20, 1023–1048. [Google Scholar] [CrossRef] [PubMed]
  27. Mutlu, B.R.; Edd, J.F.; Toner, M. Oscillatory inertial focusing in infinite microchannels. Proc. Natl. Acad. Sci. USA 2018, 115, 7682–7687. [Google Scholar] [CrossRef] [Green Version]
  28. Di Carlo, D. Inertial microfluidics. Lab Chip 2009, 9, 3038–3046. [Google Scholar] [CrossRef]
  29. Dai, Y.L.; Xu, C.; Sun, X.L.; Chen, X.Y. Nanoparticle design strategies for enhanced anticancer therapy by exploiting the tumour microenvironment. Chem. Soc. Rev. 2017, 46, 3830–3852. [Google Scholar] [CrossRef]
  30. Mitchell, M.J.; Billingsley, M.M.; Haley, R.M.; Wechsler, M.E.; Peppas, N.A.; Langer, R. Engineering precision nanoparticles for drug delivery. Nat. Rev. Drug Discov. 2021, 20, 101–124. [Google Scholar] [CrossRef] [PubMed]
  31. Aidun, C.K.; Clausen, J.R. Lattice-Boltzmann Method for Complex Flows. Annu. Rev. Fluid Mech. 2010, 42, 439–472. [Google Scholar] [CrossRef]
  32. Elghobashi, S. Direct Numerical Simulation of Turbulent Flows Laden with Droplets or Bubbles. Annu. Rev. Fluid Mech. 2019, 51, 217–244. [Google Scholar] [CrossRef] [Green Version]
  33. Ladd, A.J.C.; Verberg, R. Lattice-Boltzmann Simulations of Particle-Fluid Suspensions. J. Stat. Phys. 2001, 104, 1191–1251. [Google Scholar] [CrossRef]
  34. Rettinger, C.; Rüde, U. A comparative study of fluid-particle coupling methods for fully resolved lattice Boltzmann simulations. Comput. Fluids 2017, 154, 74–89. [Google Scholar] [CrossRef] [Green Version]
  35. Eshghinejadfard, A.; Hosseini, S.A.; Thévenin, D. Fully-resolved prolate spheroids in turbulent channel flows: A lattice Boltzmann study. AIP Adv. 2017, 7, 95007. [Google Scholar] [CrossRef] [Green Version]
  36. Rettinger, C.; Rüde, U. A coupled lattice Boltzmann method and discrete element method for discrete particle simulations of particulate flows. Comput. Fluids 2018, 172, 706–719. [Google Scholar] [CrossRef] [Green Version]
  37. Eshghinejadfard, A.; Abdelsamie, A.; Hosseini, S.A.; Thévenin, D. Immersed boundary lattice Boltzmann simulation of turbulent channel flows in the presence of spherical particles. Int. J. Multiph. Flow 2017, 96, 161–172. [Google Scholar] [CrossRef]
  38. Karimnejad, S.; Amiri Delouei, A.; Nazari, M.; Shahmardan, M.M.; Mohamad, A.A. Sedimentation of elliptical particles using Immersed Boundary – Lattice Boltzmann Method: A complementary repulsive force model. J. Mol. Liq. 2018, 262, 180–193. [Google Scholar] [CrossRef]
  39. Thorimbert, Y.; Marson, F.; Parmigiani, A.; Chopard, B.; Lätt, J. Lattice Boltzmann simulation of dense rigid spherical particle suspensions using immersed boundary method. Comput. Fluids 2018, 166, 286–294. [Google Scholar] [CrossRef]
  40. Jebakumar, A.S.; Magi, V.; Abraham, J. Lattice-Boltzmann simulations of particle transport in a turbulent channel flow. Int. J. Heat Mass Transf. 2018, 127, 339–348. [Google Scholar] [CrossRef]
  41. Peng, C.; Ayala, O.M.; Wang, L.P. A comparative study of immersed boundary method and interpolated bounce-back scheme for no-slip boundary treatment in the lattice Boltzmann method: Part I, laminar flows. Comput. Fluids 2019, 192, 104233. [Google Scholar] [CrossRef] [Green Version]
  42. Geneva, N.; Peng, C.; Li, X.M.; Wang, L.P. A scalable interface-resolved simulation of particle-laden flow using the lattice Boltzmann method. Parallel Comput. 2017, 67, 20–37. [Google Scholar] [CrossRef]
  43. Jebakumar, A.S.; Premnath, K.N.; Magi, V.; Abraham, J. Fully-resolved direct numerical simulations of particle motion in a turbulent channel flow with the lattice-Boltzmann method. Comput. Fluids 2019, 179, 238–247. [Google Scholar] [CrossRef]
  44. Qian, Y.H.; D’Humières, D.; Lallemand, P. Lattice BGK Models for Navier-Stokes Equation. EPL 1992, 17, 479–484. [Google Scholar] [CrossRef]
  45. He, X.Y.; Luo, L.S. Lattice Boltzmann Model for the Incompressible Navier–Stokes Equation. J. Stat. Phys. 1997, 88, 927–944. [Google Scholar] [CrossRef]
  46. Lallemand, P.; Luo, L.S. Lattice Boltzmann method for moving boundaries. J. Comput. Phys. 2003, 184, 406–421. [Google Scholar] [CrossRef]
  47. Mei, R.W.; Da Yu, Z.; Shyy, W.; Luo, L.S. Force evaluation in the lattice Boltzmann method involving curved geometry. Phys. Rev. E 2002, 65, 41203. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  48. Aidun, C.K.; Lu, Y.N.; Ding, E.J. Direct analysis of particulate suspensions with inertia using the discrete Boltzmann equation. J. Fluid Mech. 1998, 373, 287–311. [Google Scholar] [CrossRef]
  49. Yuan, X.F.; Ball, R.C. Rheology of hydrodynamically interacting concentrated hard disks. J. Chem. Phys. 1994, 101, 9016–9021. [Google Scholar] [CrossRef]
  50. Pirola, S.; Cheng, Z.; Jarral, O.A.; O’Regan, D.P.; Pepper, J.R.; Athanasiou, T.; Xu, X.Y. On the choice of outlet boundary conditions for patient-specific analysis of aortic flow using computational fluid dynamics. J. Biomech. 2017, 60, 15–21. [Google Scholar] [CrossRef] [PubMed]
  51. Hu, X.; Lin, J.Z.; Ku, X.K. Inertial migration of circular particles in Poiseuille flow of a power-law fluid. Phys. Fluids 2019, 31, 73306. [Google Scholar]
  52. Nie, D.M.; Lin, J.Z. Simulation of sedimentation of two spheres with different densities in a square tube. J. Fluid Mech. 2020, 896, A12. [Google Scholar] [CrossRef]
  53. Di Chen, S.; Pan, T.W.; Chang, C.C. The motion of a single and multiple neutrally buoyant elliptical cylinders in plane Poiseuille flow. Phys. Fluids 2012, 24, 103302. [Google Scholar] [CrossRef]
  54. Qian, S.Z.; Jiang, M.Q.; Liu, Z.H. Inertial migration of aerosol particles in three-dimensional microfluidic channels. Particuology 2021, 55, 23–34. [Google Scholar] [CrossRef]
  55. Barkla, H.M.; Auchterlonie, L.J. The Magnus or Robins effect on rotating spheres. J. Fluid Mech. 1971, 47, 437–447. [Google Scholar] [CrossRef]
  56. Mahmoud, G.M.; Ahmed, M.E. Chaotic and Hyperchaotic Complex Jerk Equations. Int. J. Mod. Nonlinear Theory Appl. 2012, 01, 6–13. [Google Scholar] [CrossRef] [Green Version]
  57. Schaaf, C.; Stark, H. Particle pairs and trains in inertial microfluidics. Eur. Phys. J. E 2020, 43, 50. [Google Scholar] [CrossRef]
Figure 1. D2Q9 Cartesian lattice and discrete velocities.
Figure 1. D2Q9 Cartesian lattice and discrete velocities.
Micromachines 12 01075 g001
Figure 2. Improved bounce-back scheme boundary conditions. Lighter gray circle represents the location of particle at time t and darker gray circle at time t + t , which result in two orange triangles uncover to fluid node, three purple squares cover to solid node and four red thin diamonds remain solid node. The other sky-blue circles denote the fluid nodes. The distribution function f 7 can be determined by the information of node S ,   D ,   A ,   B ,   C .
Figure 2. Improved bounce-back scheme boundary conditions. Lighter gray circle represents the location of particle at time t and darker gray circle at time t + t , which result in two orange triangles uncover to fluid node, three purple squares cover to solid node and four red thin diamonds remain solid node. The other sky-blue circles denote the fluid nodes. The distribution function f 7 can be determined by the information of node S ,   D ,   A ,   B ,   C .
Micromachines 12 01075 g002
Figure 3. Configuration of a particle migrating in Poiseuille channel flow. The origin coordinate locates at left down corner, X in horizontal direction and Y in vertical direction. The simulation domain size is fixed as L × H , while the particle is located at [ X s ,   Y s ] initially. The diameter of the particle is D p , U m represents the maximum velocity. For pulsatile case, U m will change by time.
Figure 3. Configuration of a particle migrating in Poiseuille channel flow. The origin coordinate locates at left down corner, X in horizontal direction and Y in vertical direction. The simulation domain size is fixed as L × H , while the particle is located at [ X s ,   Y s ] initially. The diameter of the particle is D p , U m represents the maximum velocity. For pulsatile case, U m will change by time.
Micromachines 12 01075 g003
Figure 4. Lateral migration of the particle released from different initial positions in Poiseuille flow with ( a )   k = 0.25 , ( b )   k = 0.35 .
Figure 4. Lateral migration of the particle released from different initial positions in Poiseuille flow with ( a )   k = 0.25 , ( b )   k = 0.35 .
Micromachines 12 01075 g004
Figure 5. Comparison of TEPs of the particle migrating in Poiseuille flow at different R e .
Figure 5. Comparison of TEPs of the particle migrating in Poiseuille flow at different R e .
Micromachines 12 01075 g005
Figure 6. The contour view of ( a ) the dimensionless pressure p , ( b ) the dimensionless fluid velocity in X direction U x and ( c ) in Y direction U y in stable state of particle migrating in non-pulsatile flow at k = 0.25 and R e = 50 .
Figure 6. The contour view of ( a ) the dimensionless pressure p , ( b ) the dimensionless fluid velocity in X direction U x and ( c ) in Y direction U y in stable state of particle migrating in non-pulsatile flow at k = 0.25 and R e = 50 .
Micromachines 12 01075 g006
Figure 7. The particle center trajectory for k = 0.25 (two sky-blue curves with circle marker), k = 0.35 (two purple curves with square marker) at R e = 50 . Solid curve represents the case driven by pulsatile velocity, which is abbreviated as P, while NP is the acronym for non-pulsatile flow donated by dash curve.
Figure 7. The particle center trajectory for k = 0.25 (two sky-blue curves with circle marker), k = 0.35 (two purple curves with square marker) at R e = 50 . Solid curve represents the case driven by pulsatile velocity, which is abbreviated as P, while NP is the acronym for non-pulsatile flow donated by dash curve.
Micromachines 12 01075 g007
Figure 8. Time history of the particle center trajectory at different R e   ( a , c ) at k = 0.25 , and different k   ( b , d ) at R e = 50 . The variation of δ , , t with R e   ( e ) and k   ( f ) .
Figure 8. Time history of the particle center trajectory at different R e   ( a , c ) at k = 0.25 , and different k   ( b , d ) at R e = 50 . The variation of δ , , t with R e   ( e ) and k   ( f ) .
Micromachines 12 01075 g008
Figure 9. The particle orientation at different R e   ( a , c )   and k   ( b , d ) . Meanwhile, ( a , b ) is in non-pulsatile and ( c , d ) in pulsatile flow.
Figure 9. The particle orientation at different R e   ( a , c )   and k   ( b , d ) . Meanwhile, ( a , b ) is in non-pulsatile and ( c , d ) in pulsatile flow.
Micromachines 12 01075 g009
Figure 10. Damping of U x (solid line in sky-blue color) of the particle at k = 0.25 and different R e = ( a )   25 ,   ( b )   50 ,   ( c )   75 ,   ( d )   100 ,   ( e )   150 ,   ( f )   200 . The upside (dashed line in purple color) and downside (dash dotted line in orange color) envelope curves are fitted by exponential function.
Figure 10. Damping of U x (solid line in sky-blue color) of the particle at k = 0.25 and different R e = ( a )   25 ,   ( b )   50 ,   ( c )   75 ,   ( d )   100 ,   ( e )   150 ,   ( f )   200 . The upside (dashed line in purple color) and downside (dash dotted line in orange color) envelope curves are fitted by exponential function.
Micromachines 12 01075 g010
Figure 11. Damping of U x at R e = 50 and different = ( a )   0.15 ,   ( b )   0.20 ,   ( c )   0.25 ,   ( d )   0.30 ,   ( e )   0.35 ,   ( f )   0.40 . Envelope curves are fitted in the same way.
Figure 11. Damping of U x at R e = 50 and different = ( a )   0.15 ,   ( b )   0.20 ,   ( c )   0.25 ,   ( d )   0.30 ,   ( e )   0.35 ,   ( f )   0.40 . Envelope curves are fitted in the same way.
Micromachines 12 01075 g011
Figure 12. The quasi frequency ( a ) and damping ratio ( b ) at different R e and k .
Figure 12. The quasi frequency ( a ) and damping ratio ( b ) at different R e and k .
Micromachines 12 01075 g012
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Huang, L.; Du, J.; Zhu, Z. Neutrally Buoyant Particle Migration in Poiseuille Flow Driven by Pulsatile Velocity. Micromachines 2021, 12, 1075. https://doi.org/10.3390/mi12091075

AMA Style

Huang L, Du J, Zhu Z. Neutrally Buoyant Particle Migration in Poiseuille Flow Driven by Pulsatile Velocity. Micromachines. 2021; 12(9):1075. https://doi.org/10.3390/mi12091075

Chicago/Turabian Style

Huang, Lizhong, Jiayou Du, and Zefei Zhu. 2021. "Neutrally Buoyant Particle Migration in Poiseuille Flow Driven by Pulsatile Velocity" Micromachines 12, no. 9: 1075. https://doi.org/10.3390/mi12091075

APA Style

Huang, L., Du, J., & Zhu, Z. (2021). Neutrally Buoyant Particle Migration in Poiseuille Flow Driven by Pulsatile Velocity. Micromachines, 12(9), 1075. https://doi.org/10.3390/mi12091075

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop