key: cord-1026919-g7divujb
authors: Huang, Lizhong; Du, Jiayou; Zhu, Zefei
title: Neutrally Buoyant Particle Migration in Poiseuille Flow Driven by Pulsatile Velocity
date: 2021-09-06
journal: Micromachines (Basel)
DOI: 10.3390/mi12091075
sha: 0aaeda5b311b27a0b1b46e86ef4c2bcded0860ef
doc_id: 1026919
cord_uid: g7divujb

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 ([Formula: see text]) and blockage ratio [Formula: see text] 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 [Formula: see text] , while [Formula: see text] has little impact. The quasi frequency of such oscillation has almost no business with [Formula: see text] and [Formula: see text]. Moreover, [Formula: see text] 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.

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 (Re = 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 Re increased and revealed additional inner annulus when Re was greater than 600. Moreover, if Re 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 The equilibrium distribution function can be calculated by [44] :

where is the weight factor with 4 9 ⁄ , ~ 1 9 ⁄ and ~ 1 36 ⁄ , is the fluid density and is fluid velocity which can be determined by:

, , 1 , .

For low Mach number, the Navier-Stokes equations can be derived from the lattice Boltzmann equation by utilizing Chapman-Enskog expansion [45] .

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 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. The equilibrium distribution function can be calculated by [44] :

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:

For low Mach number, the Navier-Stokes equations can be derived from the lattice Boltzmann equation by utilizing Chapman-Enskog expansion [45] .

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. 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 can be determined by the information of node S, D, A, B, C.

Take fluid node as an example, after the streaming step, three unknown distribution functions ( , and denoted by three blue arrows shown in Figure 2 ) need to be determined by applying improved bounce-back scheme. Therefore, the distribution function can be calculated by:

where is the nearest solid node along direction; and are the two nearest fluid nodes along direction; the blue star denotes the boundary location which is the intersection point of particle boundary and line , meanwhile, can be determined by | | | | ⁄ ; is the velocity of the boundary location ; √3 ⁄ is the speed of sound. The other unknown distribution functions of fluid nodes around particle boundary can be solved in the similar method.

Let us continue to take the boundary location 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:

where is the position of the particle. When the particle is moving in the lattice grid from time to ∆ , 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] :

Moreover, in order to avoid unphysical overlapping between the particle and the channel wall, the extra lubrication force model is needed [49] : 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.

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:

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 AS, meanwhile, q can be determined by q = |AD|/|AS|; 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.

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:

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] : Moreover, in order to avoid unphysical overlapping between the particle and the channel wall, the extra lubrication force model is needed [49] :

where ν = c 2 s (τ − 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:

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.

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 out = ρ f c 2 s [19] . No-slip boundary conditions are imposed at the top and bottom channel walls. , ℎ ℎ ,

where 0.5 Δ is the kinematic viscosity of the fluid; is the diameter of the particle;

is the particle velocity towards the wall; ℎ is the minimum gap between the particle and the wall; ℎ 1.5Δ 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:

where , , , , , , and are the translational acceleration, velocity, rotational acceleration, rotational velocity, angle, mass and moment inertia of the particle.

The configuration of one circular particle migrating in Poiseuille flow is shown in Figure 3 . A parabolic velocity profile with the maximum velocity is set at the left inlet boundary in the positive direction, and the velocity in direction is zero. For , there are two cases: non-pulsatile or pulsatile. If pulsatile, 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 , (systolic period) and zero at , (diastolic period) is utilized which can be seen from Figure 3 . In the systolic period, 0.3 s, and in the diastolic period, 0.5 s, which means one cardiac cycle lasts 0.8 s.

The unit conversion factor from lattice time to physical time is 10 . At the right boundary, the normal derivative of the velocity is zero and the pressure is set to be [19] . No-slip boundary conditions are imposed at the top and bottom channel walls. The simulation domain size is fixed as , while the particle is located at , initially. The diameter of the particle is , represents the maximum velocity. For pulsatile case, will change by time. 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, Re = 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 Xcoordinate 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.

To validate the accuracy of our LBM code, two benchmark cases are implemented below. In the first case, Re is set to be 50, and the terminal particle Reynolds number Re 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 Re 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] . If there are no additional statements, the channel length is 500 and the height is 100 . TEP is basically unchanged for different channel length 300 , 400 , 500 , 600 , 700 and 2000 . So 500 is adopted same as Wen et al. [17] . As shown in Figure 3 , the circular particle center is located at , initially. The diameter of the circular particle in this work is 15 Driven by the fluid velocity, the particle will always travel along 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 thecoordinate of the particle exceeds Δ , 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 axis in positive direction should add up one lattice spacing when performing this shift once.

To validate the accuracy of our LBM code, two benchmark cases are implemented below.

In the first case, is set to be 50, and the terminal particle Reynolds number , / is about 9.63, where , is the terminal particle velocity in direction. It is very close to while the particle is driven by pressure difference [17] . Our simulating results with 0.25 and 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 ( ⁄ 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 is 50. The second case is adopted to validate over the entire range of 20, 40, 100, 200. The particle diameter is 22, the channel width is 200, and the channel length is 1000, which leads to the blockage ratio being 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 in 2D Poiseuille flow. The foregoing results are validated only when Re is 50. The second case is adopted to validate over the entire range of Re = 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 Re in 2D Poiseuille flow. 

After validation, we first simulate the particle migrating in non-pulsatile flow at 0.25 and 50. Figure 6a -c, shows the contour view of the dimensionless pressure ⁄ , the dimensionless fluid velocity in direction ⁄ and in direction ⁄ , 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] . 

Several dominant forces acting on the particle drive it to TEP in 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 

After validation, we first simulate the particle migrating in non-pulsatile flow at k = 0.25 and Re = 50. Figure 6a -c, shows the contour view of the dimensionless pressure p = (p − p out )/p out , 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] . 

After validation, we first simulate the particle migrating in non-pulsatile flow at 0.25 and 50. Figure 6a -c, shows the contour view of the dimensionless pressure ⁄ , the dimensionless fluid velocity in direction ⁄ and in direction ⁄ , 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] . 

Several dominant forces acting on the particle drive it to TEP in 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 

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 Re = 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. in non-pulsatile flow travels in 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 0.25 and 0.35 at 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. Overall, TEPs are all closer to the channel centerline when or 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 * is time in this cardiac cycle when the particle locates the lowest position. In general, the influence of on , ∆, * is greater than . As 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 * all increase monotonously with increasing . The effect of on , ∆, * is a little more complicated. As increases, the particle inertia increases, and the drag force also increases. Shear gradient force increases with increasing , 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 increases, the maximum of , ∆ occurs at 0.2. Moreover, * was mostly unchanged with increasing . (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 Re 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 Re on δ, ∆, t * is greater than k. As Re 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 Re. 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. Micromachines 2021, 12, x FOR PEER REVIEW 9 of 15 

Observed in Figure 9 , the particle always rotates clockwise while migrating in the channel, irrespective of , , 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 100, TEP of the particle is very similar, so they rotate almost at the same speed 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 . Furthermore, the moment inertia 

Observed in Figure 9 , the particle always rotates clockwise while migrating in the channel, irrespective of Re, 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 Re < 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 Re.

Micromachines 2021, 12, x FOR PEER REVIEW 10 of 15 of the particle is proportional to the square of the particle diameter, i.e., . Consequently, by comparing Figure 9a ,b, the influence of on is much larger than . 

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 direction is preferred, because the particle velocity curve is not center symmetry in direction. For better illustration, the particle velocity in direction is normalized by: ⁄ . The curves of in Figure 10 are much like a spring-mass system which is underdamped. Figure 10a shows, when is small, damps out rapidly after several quasi periods. If is small enough, like 1, over damped is expected. Obviously,

in Figure 10a -f, the damping effect becomes weaker with increasing . We fit the upside and downside envelope curve by using exponential function for purpose. For example, . .

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 . And the absolute values of those constants decrease with increasing . The influence of is also analyzed, but it can be seen from Figure 11a -f that has almost no effect on 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: Figure 10 are much like a spring-mass system which is underdamped. Figure 10a shows, when Re is small, U x damps out rapidly after several quasi periods. If Re is small enough, like Re ≈ 1, over damped is expected. Obviously, in Figure 10a -f, the damping effect becomes weaker with increasing Re. We fit the upside and downside envelope curve by using exponential function for purpose. For example, A up = e −15.99t +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 Re. And the absolute values of those constants decrease with increasing Re. 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 ( ) of oscillation at different and . It can be seen from the subplot that the quasi frequency increases as increasing , but the impact of on the quasi frequency is relatively small. However, the quasi frequency will not change when increases except for 25. As mentioned earlier, in Figure 10a , when 25, 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 25.

As shown in Figure 12b , the damping ratio 2 ⁄ , 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 oscillation basically does not Figure 12a gives the quasi frequency ( ) of oscillation at different and . It can be seen from the subplot that the quasi frequency increases as increasing , but the impact of on the quasi frequency is relatively small. However, the quasi frequency will not change when increases except for 25. As mentioned earlier, in Figure 10a , when 25, 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 25.

As shown in Figure 12b , the damping ratio 2 ⁄ , 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 oscillation basically does not Figure 12a gives the quasi frequency ( f q ) of U x oscillation at different Re 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 Re increases except for Re = 25. As mentioned earlier, in Figure 10a , when Re = 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 Re = 25. 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 is fitted as: 0.83 . . Figure 12 . The quasi frequency and damping ratio at different and .

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 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 varies only from 25 to 200. Smaller leads to bigger SRT, thus the simulation will not be accurate. In contrast, bigger results in smaller SRT which the simulation will diverge. Furthermore, range only from 0.15 to 0.40. Because the channel height is set to be 100 lattice units, if is smaller than 0.15, the simulation resolution will not be adequate. On the other hand, if 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 and . 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. 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 Re, i.e., Re 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 Re is fitted as: ζ = 0.83Re −0.64 .

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 Re 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 Re varies only from 25 to 200. Smaller Re leads to bigger SRT, thus the simulation will not be accurate. In contrast, bigger Re 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 Re 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. 

Influence of non-Newtonian power law rheology on inertial migration of particles in channel flow

Stability condition of self-organizing staggered particle trains in channel flow

Poiseuille flow-induced vibrations of a cylinder in subcritical conditions

Dynamics of particle migration in channel flow of viscoelastic fluids

Equilibrium positions of the elasto-inertial particle migration in rectangular channel flow of Oldroyd-B viscoelastic fluids

Fundamentals of Differential Particle Inertial Focusing in Symmetric Sinusoidal Microchannels

Behaviour of macroscopic rigid spheres in Poiseuille flow Part 2. Experimental results and interpretation

Migration of rigid spheres in a two-dimensional unidirectional shear flow of a second-order fluid

The inertial lift on a spherical particle in a plane Poiseuille flow at large channel Reynolds number

Inertial migration of rigid spherical particles in Poiseuille flow

Lateral force on a rigid sphere in large-inertia laminar pipe flow

Inertial migration of a rigid sphere in three-dimensional Poiseuille flow

Equilibrium radial positions of neutrally buoyant spherical particles over the circular cross-section in Poiseuille flow

Three regimes of inertial focusing for spherical particles suspended in circular tube flows

Direct simulation of initial value problems for the motion of solid bodies in a Newtonian fluid. Part 2. Couette and Poiseuille flows

Force evaluations in lattice Boltzmann simulations with moving boundaries in two dimensions

Lattice-type-dependent momentum-exchange method for moving boundaries

An investigation on momentum exchange methods and refilling algorithms for lattice Boltzmann simulation of particulate flows

Flow between parallel walls containing the lines of neutrally buoyant circular cylinders

Inertial migration of spherical particles in circular Poiseuille flow at moderately high Reynolds numbers

Migration of finite sized particles in a laminar square channel flow from low to high Reynolds numbers

Investigation on the thermophoretic sorting for submicroparticles in a sorter with expansion-contraction microchannel

Investigation on the Thermophoresis-Coupled Inertial Sorting of Submicrometer Particles in a Microchannel

Computational inertial microfluidics: A review

Oscillatory inertial focusing in infinite microchannels

Nanoparticle design strategies for enhanced anticancer therapy by exploiting the tumour microenvironment

Engineering precision nanoparticles for drug delivery

Lattice-Boltzmann Method for Complex Flows

Direct Numerical Simulation of Turbulent Flows Laden with Droplets or Bubbles

Lattice-Boltzmann Simulations of Particle-Fluid Suspensions

A comparative study of fluid-particle coupling methods for fully resolved lattice Boltzmann simulations

Fully-resolved prolate spheroids in turbulent channel flows: A lattice Boltzmann study

A coupled lattice Boltzmann method and discrete element method for discrete particle simulations of particulate flows

Immersed boundary lattice Boltzmann simulation of turbulent channel flows in the presence of spherical particles

Sedimentation of elliptical particles using Immersed Boundary -Lattice Boltzmann Method: A complementary repulsive force model

Lattice Boltzmann simulation of dense rigid spherical particle suspensions using immersed boundary method

Lattice-Boltzmann simulations of particle transport in a turbulent channel flow

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

A scalable interface-resolved simulation of particle-laden flow using the lattice Boltzmann method

Fully-resolved direct numerical simulations of particle motion in a turbulent channel flow with the lattice-Boltzmann method

Lattice BGK Models for Navier-Stokes Equation

Lattice Boltzmann Model for the Incompressible Navier-Stokes Equation

Lattice Boltzmann method for moving boundaries

Force evaluation in the lattice Boltzmann method involving curved geometry

Direct analysis of particulate suspensions with inertia using the discrete Boltzmann equation

Rheology of hydrodynamically interacting concentrated hard disks

On the choice of outlet boundary conditions for patient-specific analysis of aortic flow using computational fluid dynamics

Inertial migration of circular particles in Poiseuille flow of a power-law fluid

Simulation of sedimentation of two spheres with different densities in a square tube

The motion of a single and multiple neutrally buoyant elliptical cylinders in plane Poiseuille flow

Inertial migration of aerosol particles in three-dimensional microfluidic channels

The Magnus or Robins effect on rotating spheres

Chaotic and Hyperchaotic Complex Jerk Equations

Particle pairs and trains in inertial microfluidics

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

The authors report that they have no conflict of interest.