*Article* **Formation and Elimination of Satellite Droplets during Monodisperse Droplet Generation by Using Piezoelectric Method**

**Zejian Hu <sup>1</sup> , Shengji Li 2, \* , Fan Yang 1 , Xunjie Lin 1 , Sunqiang Pan 3 , Xuefeng Huang 1, \* and Jiangrong Xu 1**

> 1 Institute of Energy, Department of Physics, Hangzhou Dianzi University, Hangzhou 310018, China; zejianhu@hdu.edu.cn (Z.H.); yf@hdu.edu.cn (F.Y.); 17072215@hdu.edu.cn (X.L.); jrxu@hdu.edu.cn (J.X.)

<sup>2</sup> College of Materials and Environmental Engineering, Hangzhou Dianzi University, Hangzhou 310018, China <sup>3</sup> Division of Biological and Chemical Metrology, Zhejiang Institute of Metrology, Hangzhou 310018, China;

pansunqiang@hotmail.com **\*** Correspondence: shengjili@hdu.edu.cn (S.L.); xuefenghuang@hdu.edu.cn (X.H.); Tel.: +86-571-8687-8677 (X.H.)

**Abstract:** One of the key questions in the generation of monodisperse droplets is how to eliminate satellite droplets. This paper investigates the formation and elimination of satellite droplets during the generation of monodisperse deionized water droplets based on a piezoelectric method. We estimated the effects of two crucial parameters—the pulse frequency for driving the piezoelectric transducer (PZT) tube and the volume flow rate of the pumping liquid—on the generation of monodisperse droplets of the expected size. It was found that by adjusting the pulse frequency to harmonize with the volume flow rate, the satellite droplets can be eliminated through their coalescence with the subsequent mother droplets. An increase in the tuning pulse frequency led to a decrease in the size of the monodisperse droplets generated. Among three optimum conditions (OCs) (OC1: 20 mL/h, 20 kHz; OC2: 30 mL/h, 30 kHz; and OC3: 40 mL/h, 40 kHz), the sizes of the generated monodisperse deionized water droplets followed a bimodal distribution in OC1 and OC2, whereas they followed a Gaussian distribution in OC3. The average diameters were 87.8 µm (OC1), 85.9 µm (OC2), and 84.8 µm (OC3), which were 8.46%, 6.14%, and 4.69% greater than the theoretical one (81.0 µm), respectively. This monodisperse droplet generation technology is a promising step in the production of monodisperse aerosols for engineering applications.

**Keywords:** monodisperse droplet generation; satellite droplets; piezoelectric method; droplet coalescence

#### **1. Introduction**

In recent years, monodisperse droplet generation technology has been widely used in the fields of additive manufacturing [1,2], inkjet printing [3,4], electronic packaging [5,6], bioengineering [7–9], instrument calibration [10,11], etc. It has prompted many researchers to become engaged in developing droplet generation techniques to meet new requirements such as droplets that are highly uniform and monodisperse in terms of their size, shape, density, and surface characteristics, with a variety of solutes and solvents.

Many attempts have been made to generate monodisperse droplets, such as hot bubble [12], mechanical [13], pneumatic [14,15], piezoelectric [16], electromagnetic [17], and droplet-based microfluidic [18–20] technologies. Among these, the piezoelectric droplet generation method is one of the best choices to obtain monodisperse droplets. On the basis of the piezoelectric method, the generation of droplets depends on the control of the pulse waveform for the driving of the PZT tube. Li et al. [21] reported that monodisperse droplets can be obtained through adjusting the frequency and amplitude of a rectangular pulse waveform at a high operating pressure of 3.5 MPa. Fan et al. [22] reported that monodisperse droplets were generated by controlling the upper and lower limits of the pulse

**Citation:** Hu, Z.; Li, S.; Yang, F.; Lin, X.; Pan, S.; Huang, X.; Xu, J. Formation and Elimination of Satellite Droplets during Monodisperse Droplet Generation by Using Piezoelectric Method. *Micromachines* **2021**, *12*, 921. https:// doi.org/10.3390/mi12080921

Academic Editors: Junfeng Zhang and Ruijin Wang

Received: 10 July 2021 Accepted: 28 July 2021 Published: 31 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

amplitude. However, during droplet generation, many satellite droplets were produced, leading to a wide size distribution of the droplets. Shin et al. [23] also reported that, for low viscosity liquid, satellites were generated when a single-pulse waveform was applied, whereas when a double-pulse waveform was utilized, the satellites were eliminated. A simple change in the time separation can precisely control the droplet size. If the time separation is shortened, the droplets' size becomes smaller, since the second pulse reduces the mass and momentum of the ejected liquid. However, the double-pulse waveform can possibly result in the coalescence of two adjacent mother droplets. Lin et al. [24] further studied the influence of pulse frequency, positive voltage time, and voltage magnitude on droplet ejection velocity using a double-pulse voltage pattern, but they estimated it to have a scarce effect on droplet size.

Regardless of the control associated with the single-pulse waveform or the doublepulse waveform, the breakup of fluid filaments ejected from the nozzle leads to an array of uniformly spaced large droplets (called mother droplets) with smaller droplets (known as satellite droplets) in between them. Therefore, to generate uniform mother droplets based on piezoelectric method, the problem of eliminating a large number of satellite droplets has to be solved. Otherwise, this situation would result in a nonuniform size distribution of the breakup droplets. There is no doubt that one of the key steps involved in generating uniform and monodisperse droplets is eliminating the satellite droplets and avoiding the coalescence of two adjacent mother droplets.

In this work, we have aimed to investigate the formation and elimination of satellite droplets during the generation of monodisperse deionized water droplets based on the piezoelectric method. The emphasis of this study was on harmonizing the relationship between the frequency of the square-pulse waveform and the volume flow rate in order to obtain monodisperse droplets of the expected size. Therefore, this work involves: (1) observing the formation and elimination processes of satellite droplets during the generaton of monodisperse deionized water droplets through a high-resolution imaging system; (2) estimating the effect of pulse frequency and volume flow rate parameters on the droplet size and its distribution; (3) revealing the mechanisms involved in the generation and elimination of satellite droplets related to these two crucial parameters.

#### **2. Experimental Methods**

A schematic of the experimental setup is shown in Figure 1a. It was mainly composed of a micropump and syringe, filter, nozzle, controller, high-speed camera, and an LED lamp. The micropump was used to feed the deionized water stored in the syringe and to control the volume flow rate with an accuracy of ± 2%. The deionized water (the properties of which are listed in Table 1) was transported to the nozzle through the connecting pipe and filter. The droplets were then extruded by means of the PZT tube in the nozzle, controlled by the controller. The high-speed camera (Phantom M310, Vision Research Inc., Wayne, NJ, USA) with a lens (AT-X M100 AF PRO, Tokina, Japan) was utilized to record the pinchoff process of the liquid filaments, the formation of mother droplets, and the generation of satellite droplets. The recording frame was set to 5000 fps. The LED lamp with an adjustable luminance was used to illuminate the field of view for clear imaging using the high-speed camera.

The nozzle (MDG100) was purchased from TSI Co. Ltd., USA, and was based on a squeeze mode design. A PZT tube was wrapped outside a glass tubular reservoir, shown in Figure 1b. The orifice diameter at the bottom of the nozzle was 50 µm. The PZT tube was driven by the controller using a transistor–transistor logic (TTL) signal to modulate the width and frequency of periodic rectangle-wave pulse as shown in Figure 2. The period of the pulse waveform is *T*, and its voltage amplitude is *V*. The duty ratio is *Tp*/*T*, where *T<sup>p</sup>* is the applied time of the pulse. The single-pulse waveform is a rectangle wave, which is defined by a pulse width, a rising edge, and a falling edge. The applying time of the pulse can be divided into three parts—the rising time (*tr*), dwelling time (*t<sup>d</sup>* ), and falling time (*tf* ). When the rising (falling) edge of the waveform is applied to the PZT tube, a negative (positive) pressure wave is produced, expanding (contracting) the glass tubular reservoir, so that the fluid filaments will be ejected and broken up. These parameters of the driving pulse for the PZT tube are listed in Table 2. negative (positive) pressure wave is produced, expanding (contracting) the glass tubular reservoir, so that the fluid filaments will be ejected and broken up. These parameters of the driving pulse for the PZT tube are listed in Table 2. was driven by the controller using a transistor–transistor logic (TTL) signal to modulate the width and frequency of periodic rectangle-wave pulse as shown in Figure 2. The period of the pulse waveform is *T*, and its voltage amplitude is *V*. The duty ratio is *Tp*/*T*, where *T<sup>p</sup>* is the applied time of the pulse. The single-pulse waveform is a rectangle wave,

The nozzle (MDG100) was purchased from TSI Co. Ltd., USA, and was based on a squeeze mode design. A PZT tube was wrapped outside a glass tubular reservoir, shown in Figure 1b. The orifice diameter at the bottom of the nozzle was 50 μm. The PZT tube was driven by the controller using a transistor–transistor logic (TTL) signal to modulate the width and frequency of periodic rectangle-wave pulse as shown in Figure 2. The period of the pulse waveform is *T*, and its voltage amplitude is *V*. The duty ratio is *Tp*/*T*, where *T<sup>p</sup>* is the applied time of the pulse. The single-pulse waveform is a rectangle wave, which is defined by a pulse width, a rising edge, and a falling edge. The applying time of the pulse can be divided into three parts—the rising time (*tr*), dwelling time (*td*), and falling time (*tf*). When the rising (falling) edge of the waveform is applied to the PZT tube, a

The nozzle (MDG100) was purchased from TSI Co. Ltd., USA, and was based on a squeeze mode design. A PZT tube was wrapped outside a glass tubular reservoir, shown in Figure 1b. The orifice diameter at the bottom of the nozzle was 50 μm. The PZT tube

*Micromachines* **2021**, *12*, 921 3 of 12

*Micromachines* **2021**, *12*, 921 3 of 12

**Figure 1.** (**a**) Schematic of experimental setup; (**b**) schematic of nozzle. **Figure 1.** (**a**) Schematic of experimental setup; (**b**) schematic of nozzle.

**Figure 1.** (**a**) Schematic of experimental setup; (**b**) schematic of nozzle.

**Table 1.** Physical properties of deionized water.


**Figure 2.** Schematic of periodic rectangle-wave pulse waveform. **Figure 2.** Schematic of periodic rectangle-wave pulse waveform.


of monodisperse droplets. In this work, a digital image processing method was used. The

The measurement of droplet size is crucial to estimate performance in the generation of monodisperse droplets. In this work, a digital image processing method was used. The elaborate operation procedures have been presented in our previous work [25]. Briefly, (1) a magnified image of the calibrating ruler was acquired; (2) pictures of the breakup droplets were taken under the same imaging conditions; (3) a self-programming digital imaging treatment program was then used to compare the images of the breakup droplets with calibrated pixel pitches. If the droplets appeared ellipsoidal or non-spherical due to deformation, a characteristic dimension *d*<sup>21</sup> (called the equivalent diameter) was calculated by *d*<sup>21</sup> = 4*Sd*/*L<sup>d</sup>* , where *S<sup>d</sup>* and *L<sup>d</sup>* represent the area and the perimeter of each droplet, respectively. Finally, the statistical average diameter of multiple droplets (*d*) was determined by the equation *d* = ∑ *nid<sup>i</sup>* ∑ *n<sup>i</sup>* . elaborate operation procedures have been presented in our previous work [25]. Briefly, (1) a magnified image of the calibrating ruler was acquired; (2) pictures of the breakup droplets were taken under the same imaging conditions; (3) a self-programming digital imaging treatment program was then used to compare the images of the breakup droplets with calibrated pixel pitches. If the droplets appeared ellipsoidal or non-spherical due to deformation, a characteristic dimension *d*21 (called the equivalent diameter) was calculated by *d*21 = 4*Sd*/*Ld*, where *Sd* and *Ld* represent the area and the perimeter of each droplet, respectively. Finally, the statistical average diameter of multiple droplets (*d*) was determined by the equation *i i i n d d n* .

#### **3. Results and Discussion 3. Results and Discussion**

#### *3.1. Pinch-Off of Liquid Filament and Generation of Satellite Droplets 3.1. Pinch-Off of Liquid Filament and Generation of Satellite Droplets*

*Micromachines* **2021**, *12*, 921 4 of 12

**Table 2.** Parameters of pulse waveform for driving the PZT tube.

Figure 3 demonstrates a snapshot taken during the pinch-off process of the liquid filament at a pulse frequency of 10 kHz and a volume flow rate of 30 mL/h. It can be observed that surface waves cover the liquid filament, and a mother droplet has been produced; another droplet is being generating at the liquid neck (the narrowest position along the liquid filament). At the bottom of the picture, when a mother droplet is detached from the liquid filament, a satellite droplet is simultaneously generated. Figure 3 demonstrates a snapshot taken during the pinch-off process of the liquid filament at a pulse frequency of 10 kHz and a volume flow rate of 30 mL/h. It can be observed that surface waves cover the liquid filament, and a mother droplet has been produced; another droplet is being generating at the liquid neck (the narrowest position along the liquid filament). At the bottom of the picture, when a mother droplet is detached from the liquid filament, a satellite droplet is simultaneously generated.

**Figure 3.** Snapshot of the pinch-off of a deionized water filament and the generation of satellite **Figure 3.** Snapshot of the pinch-off of a deionized water filament and the generation of satellite droplets.

droplets. The Reynolds and Weber numbers can be expressed as *Re* = *ρvd/μ* and *We* = *ρv2d/σ*, in which *ρ*, *v*, *d*, *σ* and *μ* represent the density, velocity, orifice diameter, surface tension, and viscosity of the deionized water, respectively. In this case, the Reynolds number was ~186.2 and the Weber number was ~1.4. This suggests that the pinch-off of the deionized water filament generally exhibited a laminar regime, and the surface tension played a key The Reynolds and Weber numbers can be expressed as *Re* = *ρvd/µ* and *We* = *ρv <sup>2</sup>d/σ*, in which *ρ*, *v*, *d*, *σ* and *µ* represent the density, velocity, orifice diameter, surface tension, and viscosity of the deionized water, respectively. In this case, the Reynolds number was ~186.2 and the Weber number was ~1.4. This suggests that the pinch-off of the deionized water filament generally exhibited a laminar regime, and the surface tension played a key role in the detachment.

role in the detachment. The evolution of the liquid filament was considered as a function of the viscosity ratio (*p*) of the fluids and the initial wavenumber of the interface perturbation, and the satellite The evolution of the liquid filament was considered as a function of the viscosity ratio (*p*) of the fluids and the initial wavenumber of the interface perturbation, and the satellite droplets were generated around the neck region of a highly deformed filament. Tjahjadt et al. [26] numerically and experimentally explained that, in low-viscosity ratio systems, *p* < *O* (0.1), the breakup mechanism depends on self-repetition in multiple breakup

sequences in which every pinch-off is always associated with the formation of a neck; the neck undergoes pinch-off, and the process repeats. In this case, unlike multiple breakup sequences, at each period, we observed one breakup sequence and a generated satellite droplet. This means that, at each corresponding period of the controlled rectangle-wave pulse of PZT tube, the duty ratio, the frequency, and the sequence of time applied will be crucial in order to affect the formation of mother and satellite droplets.

#### *3.2. Elimination of Satellite Droplets*

Theoretically, if the liquid filament breaks up into monodisperse mother droplets without the generation of satellite droplets, the input volume per unit time should be equal to the total volume of the monodisperse droplets (mother droplets) generated per unit time:

$$q \times \frac{10^{-6}}{3600} = \frac{\pi}{6} d\_{\rm th}^3 \times 10^{-18} \times f \times 10^3 \tag{1}$$

In which *dth*, *q*, and *f* represent the theoretical average diameter of the generated monodisperse droplets (µm), the volume flow rate of the feeding liquid (mL/h), and the pulse frequency for the driving of the PZT tube (kHz), respectively.

Thus, the theoretical average droplet diameter of the generated monodisperse droplets in an ideal state can be calculated as follows:

$$d\_{th} = \sqrt[3]{\frac{q}{600\pi f}} \times 10^3 \tag{2}$$

The volume flow rate of deionized water was still set to 30 mL/h. As the pulse frequency was adjusted from 10 kHz to 30 kHz, we observed the coalescence of satellite droplets and mother droplets, as shown in Figure 4. During the pinch-off process of the liquid filament, a satellite droplet was generated (Figure 4a), but it immediately merged with the subsequent mother droplet (Figure 4b,c). Finally, monodisperse droplets were generated with a highly uniform size distribution (Figure 4d). The average diameter of the generated monodisperse droplets was 85.9 µm in this case, which is 6.14% larger than the theoretical one. This suggests that the adjustment of the pulse frequency could eliminate the satellite droplets by means of coalescence, validating the feasibility of this methodology. *Micromachines* **2021**, *12*, 921 6 of 12 This indicates that the slope of the voltage variation should be less than 15 V/μs. Additionally, the voltage amplitude should be high enough to ensure that the droplet is ejected, but not so high that the ejection becomes chaotic and it becomes difficult to obtain monodisperse droplets.

**Figure 4.** Snapshots of the pinch-offs of liquid filaments and the elimination of satellite droplets. (**a**) Formation of satellite droplets. (**b**,**c**) Satellite droplets merged by a subsequent mother droplet and eliminated; (**d**) a highly uniform size distribution of monodisperse droplets were obtained (volume flow rate of 30 mL/h; pulse frequencies of 30 kHz). **Figure 4.** Snapshots of the pinch-offs of liquid filaments and the elimination of satellite droplets. (**a**) Formation of satellite droplets. (**b**,**c**) Satellite droplets merged by a subsequent mother droplet and eliminated; (**d**) a highly uniform size distribution of monodisperse droplets were obtained (volume flow rate of 30 mL/h; pulse frequencies of 30 kHz).

**Figure 5.** Schematic of the theoretical and practical generation of monodisperse droplets.

To generate monodisperse droplets, the rising time (*tr*) and the falling time (*tf)* should also harmonize with the detaching time of each droplet. Under the same operating

It is worthy of notice that the number of monodisperse droplets is theoretically equal to the pulse frequency, i.e., only a single droplet is generated as a pulse is applied to the PZT tube. Thus, the pulse width and the duty ratio of each pulse need to be negotiated for the detachment of each droplet from the liquid filament, as shown in Figure 5. The generation time of monodisperse droplets (*tgen*) is equal to the theoretical pulse width, or the applied time of the practical pulse, i.e., *tgen*=*Tp*=*t<sup>r</sup>* + *t<sup>d</sup>* + *t<sup>f</sup>* , in which *t<sup>r</sup>* , *t<sup>d</sup>* , and *t<sup>f</sup>* are the rising time, dwelling time, and falling time, respectively. The flight time of the monodisperse droplets (*tfly*) is the same as the vacant time of the pulse, i.e., *tfly*=*T* − *Tp*, in which *T* and *T<sup>p</sup>* are the period of the pulse waveform and the applied time of the pulse, respectively. Lin et al. [24] suggested that a sufficient time of *t<sup>r</sup>* and *t<sup>f</sup>* is necessary to reach the desired voltage amplitude and to make the PZT tube expand and contract enough. The empirical expression of the voltage amplitude of the pulse waveform must meet the following requirements [24]: **Figure 4.** Snapshots of the pinch-offs of liquid filaments and the elimination of satellite droplets. (**a**) Formation of satellite droplets. (**b**,**c**) Satellite droplets merged by a subsequent mother droplet and eliminated; (**d**) a highly uniform size distribution of monodisperse droplets were obtained (volume

This indicates that the slope of the voltage variation should be less than 15 V/μs. Additionally, the voltage amplitude should be high enough to ensure that the droplet is ejected, but not so high that the ejection becomes chaotic and it becomes difficult to obtain

$$V/t\_r \le 15, V/t\_f \le 15\tag{3}$$

*Micromachines* **2021**, *12*, 921 6 of 12

monodisperse droplets.

**Figure 5.** Schematic of the theoretical and practical generation of monodisperse droplets. **Figure 5.** Schematic of the theoretical and practical generation of monodisperse droplets.

To generate monodisperse droplets, the rising time (*tr*) and the falling time (*tf)* should also harmonize with the detaching time of each droplet. Under the same operating This indicates that the slope of the voltage variation should be less than 15 V/µs. Additionally, the voltage amplitude should be high enough to ensure that the droplet is ejected, but not so high that the ejection becomes chaotic and it becomes difficult to obtain monodisperse droplets.

To generate monodisperse droplets, the rising time (*tr*) and the falling time (*t<sup>f</sup>* ) should also harmonize with the detaching time of each droplet. Under the same operating conditions, the detaching time depends mainly on the viscosity of the liquid. The larger the viscosity, the longer the detaching time becomes. The dwelling time (*t<sup>d</sup>* ) has been given by Bogy et al. [27], i.e., the optimum pulse width (*topt*), which can be calculated as follows:

$$t\_d = t\_{opt} = L/v\_{acc,liq} \tag{4}$$

in which *L* and *vaco,liq* are the length of nozzle and the velocity of acoustic wave propagation in the liquid, respectively. In this work, the theoretical value of the optimum dwelling time (*t<sup>d</sup>* ) was ~12 µs. In Figure 4d, the central spacing between two adjacent droplets was

twice as large as the droplet diameter, since the duty ratio of the rectangular pulse was 0.5. To avoid the coalescence of two spherical mother droplets without satellite droplets, the central spacing must be greater than the droplet diameter. However, since the droplets exhibited severe deformation after being detached from the liquid filament (Figure 4a–c), the maximum ratio of droplet length to width reached 1.8, and the central spacing was close to the droplet diameter. This suggests that the duty ratio should be below 0.8. in the liquid, respectively. In this work, the theoretical value of the optimum dwelling time (*td*) was ~12 μs. In Figure 4d, the central spacing between two adjacent droplets was twice as large as the droplet diameter, since the duty ratio of the rectangular pulse was 0.5. To avoid the coalescence of two spherical mother droplets without satellite droplets, the central spacing must be greater than the droplet diameter. However, since the droplets exhibited severe deformation after being detached from the liquid filament (Figure 4a–c), the maximum ratio of droplet length to width reached 1.8, and the central spacing was

*Micromachines* **2021**, *12*, 921 7 of 12

conditions, the detaching time depends mainly on the viscosity of the liquid. The larger the viscosity, the longer the detaching time becomes. The dwelling time (*td*) has been given by Bogy et al. [27], i.e., the optimum pulse width (*topt*), which can be calculated as follows:

*d opt aco liq* ,

in which *L* and *vaco,liq* are the length of nozzle and the velocity of acoustic wave propagation

*t t L v* (4)

#### *3.3. Effect of Pulse Frequency on the Elimination of Satellite Droplets* close to the droplet diameter. This suggests that the duty ratio should be below 0.8.

Figure 6 demonstrates the snapshots of droplet generation as the pulse frequency ranged from 5 to 45 kHz at the volume flow rate of 30 mL/h. At the pulse frequencies of 5 kHz, 10 kHz, and 15 kHz, it was observed that the satellite droplets were difficult to eliminate (Figure 6b3,c2), and the two adjacent mother droplets coalesced (Figure 6a5). When the pulse frequency was enhanced to 20 kHz or 25 kHz, satellite droplets with greater size were generated (Figure 6d4,e3). At 30 kHz and 35 kHz, the satellite droplets were almost eliminated (Figure 6f3,g3). However, as the pulse frequency exceeded 40 kHz, the generated droplets were chaotic, and they had a wide size distribution. *3.3. Effect of Pulse Frequency on the Elimination of Satellite Droplets*  Figure 6 demonstrates the snapshots of droplet generation as the pulse frequency ranged from 5 to 45 kHz at the volume flow rate of 30 mL/h. At the pulse frequencies of 5 kHz, 10 kHz, and 15 kHz, it was observed that the satellite droplets were difficult to eliminate (Figure 6b3,c2), and the two adjacent mother droplets coalesced (Figure 6a5). When the pulse frequency was enhanced to 20 kHz or 25 kHz, satellite droplets with greater size were generated (Figure 6d4,e3). At 30 kHz and 35 kHz, the satellite droplets were almost eliminated (Figure 6f3,g3). However, as the pulse frequency exceeded 40 kHz, the generated droplets were chaotic, and they had a wide size distribution.

**Figure 6.** Snapshots of droplet generation at the volume flow rate of 30 mL/h at different pulse frequencies ((**a1**–**a5**): 5 kHz; (**b1**–**b3**): 10 kHz; (**c1**–**c2**): 15 kHz; (**d1**–**d4**): 20 kHz; (**e1**–**e3**): 25 kHz; (**f1**–**f3**): 30 kHz; (**g1**–**g3**): 35 kHz; (**h1**–**h2**): 40 kHz; (**i1**– **i2**): 45 kHz). **Figure 6.** Snapshots of droplet generation at the volume flow rate of 30 mL/h at different pulse frequencies ((**a1**–**a5**): 5 kHz; (**b1**–**b3**): 10 kHz; (**c1**–**c2**): 15 kHz; (**d1**–**d4**): 20 kHz; (**e1**–**e3**): 25 kHz; (**f1**–**f3**): 30 kHz; (**g1**–**g3**): 35 kHz; (**h1**–**h2**): 40 kHz; (**i1**–**i2**): 45 kHz).

Experiments on the influence of pulse frequency on the formation and elimination of satellite droplets were also carried out at the volume flow rates of 20 mL/h and 40 mL/h, respectively. According to the results, it can be concluded that, at a certain volume flow rate, the pulse frequency must be adjusted to an optimum range to generate monodisperse droplets. At the volume flow rate of 20 mL/h, the optimum range of the pulse frequency was 18–22 kHz. For the volume flow rates of 30 mL/h and 40 mL/h, the optimum ranges of the pulse frequency were 25–35 kHz and 35–50 kHz, respectively. Experiments on the influence of pulse frequency on the formation and elimination of satellite droplets were also carried out at the volume flow rates of 20 mL/h and 40 mL/h, respectively. According to the results, it can be concluded that, at a certain volume flow rate, the pulse frequency must be adjusted to an optimum range to generate monodisperse droplets. At the volume flow rate of 20 mL/h, the optimum range of the pulse frequency was 18–22 kHz. For the volume flow rates of 30 mL/h and 40 mL/h, the optimum ranges of the pulse frequency were 25–35 kHz and 35–50 kHz, respectively.

#### *3.4. Effect of Pulse Frequency on the Average Diameter of Droplets*

According to Equation (2), the diameter of the generated droplets depends on the volume flow rate and the pulse frequency. As the volume flow is fixed, the desired diameter of the droplets can be obtained through adjusting the pulse frequency. Figure 7 shows the average diameters of the generated droplets at the volume flow rate of 30 mL/h, at different pulse frequencies of 25–40 kHz, with an interval of 5 kHz. The theoretical diameter of the droplets was calculated using Equation (2), and these are also shown in Figure 7. In general, the experimental result displayed a similar trend to the theoretical one. As the pulse frequency for driving the PZT tube was enhanced, the average diameter of the generated droplets decreased. The average droplet diameters based on the experiments were 5–11% greater than those based on the theoretical calculation. This difference can be attributed to the random error associated with dealing with the images and the systematic calibration error.

*Micromachines* **2021**, *12*, 921 8 of 12

*3.4. Effect of Pulse Frequency on the Average Diameter of Droplets* 

with the images and the systematic calibration error.

According to Equation (2), the diameter of the generated droplets depends on the volume flow rate and the pulse frequency. As the volume flow is fixed, the desired diameter of the droplets can be obtained through adjusting the pulse frequency. Figure 7 shows the average diameters of the generated droplets at the volume flow rate of 30 mL/h, at different pulse frequencies of 25–40 kHz, with an interval of 5 kHz. The theoretical diameter of the droplets was calculated using Equation (2), and these are also shown in Figure 7. In general, the experimental result displayed a similar trend to the theoretical one. As the pulse frequency for driving the PZT tube was enhanced, the average diameter of the generated droplets decreased. The average droplet diameters based on the experiments were 5%–11% greater than those based on the theoretical calculation. This difference can be attributed to the random error associated with dealing

**Figure 7.** Average droplet diameters at the volume flow rate of 30 mL/h at different pulse frequencies of 25 kHz–40 kHz. **Figure 7.** Average droplet diameters at the volume flow rate of 30 mL/h at different pulse frequencies of 25 kHz–40 kHz.

#### *3.5. Size Distribution of Monodisperse Droplets under Optimum Operating Conditions 3.5. Size Distribution of Monodisperse Droplets under Optimum Operating Conditions*

If the ratio of the volume flow rate and the optimum pulse frequency are kept constant, based on Equation (2), the diameter of the generated droplet would be theoretically the same. This was also supported by the experimental results. Figure 8 illustrates the snapshots of the generated monodisperse droplets under three optimum conditions (OCs) (OC1: 20 mL/h, 20 kHz; OC2: 30 mL/h, 30 kHz; and OC3: 40 mL/h, 40 kHz). In the optimum pulse frequency range, no satellite droplets were observed. If the ratio of the volume flow rate and the optimum pulse frequency are kept constant, based on Equation (2), the diameter of the generated droplet would be theoretically the same. This was also supported by the experimental results. Figure 8 illustrates the snapshots of the generated monodisperse droplets under three optimum conditions (OCs) (OC1: 20 mL/h, 20 kHz; OC2: 30 mL/h, 30 kHz; and OC3: 40 mL/h, 40 kHz). In the optimum pulse frequency range, no satellite droplets were observed. *Micromachines* **2021**, *12*, 921 9 of 12

distribution.

kHz).

**Figure 8.** Snapshots of monodisperse droplets generated under three optimum conditions (OC1: 20 mL/h, 20 kHz; OC2: 30 mL/h, 30 kHz; and OC3: 40 mL/h, 40 kHz). **Figure 8.** Snapshots of monodisperse droplets generated under three optimum conditions (OC1: 20 mL/h, 20 kHz; OC2: 30 mL/h, 30 kHz; and OC3: 40 mL/h, 40 kHz).

The size distributions of t**he** monodisperse droplets generated under the three optimum conditions were further extracted and analyzed in detail, as shown in Figures 9–11. Under OC1, the size of generated monodisperse droplets were mainly distributed in the range of 80–95 μm, with a small number of droplets of 95–110 μm (Figure 9). As for OC2, the size distribution of generated monodisperse droplets ranged from 80 to 93 μm The size distributions of the monodisperse droplets generated under the three optimum conditions were further extracted and analyzed in detail, as shown in Figures 9–11. Under OC1, the size of generated monodisperse droplets were mainly distributed in the range of 80–95 µm, with a small number of droplets of 95–110 µm (Figure 9). As for OC2, the size distribution of generated monodisperse droplets ranged from 80 to 93 µm

(Figure 10). In OC3, the droplet size was mainly distributed in the range of 75–95 μm

droplets had a better dispersion uniformity, and the droplet size followed a Gaussian

**Figure 9.** Size distribution of generated monodisperse droplets under OC1 (*q* = 20 mL/h, *f* = 20

(Figure 10). In OC3, the droplet size was mainly distributed in the range of 75–95 µm (Figure 11). We observed two-peak profiles for both OC1 and OC2. Compared with the bimodal size distribution in OC1 and OC2, under OC3, the generated monodisperse droplets had a better dispersion uniformity, and the droplet size followed a Gaussian distribution. (Figure 10). In OC3, the droplet size was mainly distributed in the range of 75–95 μm (Figure 11). We observed two-peak profiles for both OC1 and OC2. Compared with the bimodal size distribution in OC1 and OC2, under OC3, the generated monodisperse droplets had a better dispersion uniformity, and the droplet size followed a Gaussian distribution.

**Figure 8.** Snapshots of monodisperse droplets generated under three optimum conditions (OC1: 20

The size distributions of t**he** monodisperse droplets generated under the three optimum conditions were further extracted and analyzed in detail, as shown in Figures 9–11. Under OC1, the size of generated monodisperse droplets were mainly distributed in the range of 80–95 μm, with a small number of droplets of 95–110 μm (Figure 9). As for OC2, the size distribution of generated monodisperse droplets ranged from 80 to 93 μm

*Micromachines* **2021**, *12*, 921 9 of 12

mL/h, 20 kHz; OC2: 30 mL/h, 30 kHz; and OC3: 40 mL/h, 40 kHz).

**Figure 9.** Size distribution of generated monodisperse droplets under OC1 (*q* = 20 mL/h, *f* = 20 **Figure 9.** Size distribution of generated monodisperse droplets under OC1 (*q* = 20 mL/h, *f* = 20 kHz).

kHz).

kHz).

kHz).

respectively.

**Figure 10.** Size distribution of generated monodisperse droplets under OC2 (*q* = 30 mL/h, *f* = 30 **Figure 10.** Size distribution of generated monodisperse droplets under OC2 (*q* = 30 mL/h, *f* = 30 kHz).

**Figure 11.** Size distribution of generated monodisperse droplets under OC3 (*q* = 40 mL/h, *f* = 40

The droplets generated under the three optimum conditions were generally well dispersed and had a good uniformity of dispersion. The average diameters of the generated monodisperse droplets were 87.8 μm (OC1), 85.9 μm (OC2), and 84.8 μm (OC3), respectively, as shown in Figure 12. Compared to the theoretical diameter (81.0 μm), the experimental results were 8.46%, 6.14%, and 4.69% greater than the theoretical one, kHz).

kHz).

*Micromachines* **2021**, *12*, 921 10 of 12

**Figure 11.** Size distribution of generated monodisperse droplets under OC3 (*q* = 40 mL/h, *f* = 40 **Figure 11.** Size distribution of generated monodisperse droplets under OC3 (*q* = 40 mL/h, *f* = 40 kHz).

**Figure 10.** Size distribution of generated monodisperse droplets under OC2 (*q* = 30 mL/h, *f* = 30

The droplets generated under the three optimum conditions were generally well dispersed and had a good uniformity of dispersion. The average diameters of the generated monodisperse droplets were 87.8 μm (OC1), 85.9 μm (OC2), and 84.8 μm (OC3), respectively, as shown in Figure 12. Compared to the theoretical diameter (81.0 μm), the The droplets generated under the three optimum conditions were generally well dispersed and had a good uniformity of dispersion. The average diameters of the generated monodisperse droplets were 87.8 µm (OC1), 85.9 µm (OC2), and 84.8 µm (OC3), respectively, as shown in Figure 12. Compared to the theoretical diameter (81.0 µm), the experimental results were 8.46%, 6.14%, and 4.69% greater than the theoretical one, respectively. *Micromachines* **2021**, *12*, 921 11 of 12

experimental results were 8.46%, 6.14%, and 4.69% greater than the theoretical one,

**Figure 12.** Average diameters of theoretically and experimentally generated monodisperse droplets. **Figure 12.** Average diameters of theoretically and experimentally generated monodisperse droplets.

#### **4. Conclusions 4. Conclusions**

optimum operating conditions.

Market Regulation (20200103).

Experiments on the generation of monodisperse deionized water droplets based on the piezoelectric method were conducted. Our observations demonstrated that satellite droplets were generated because the mother droplets detached from the liquid filament with surface waveform perturbations. As the frequency of the square-wave pulse was tuned to match the volume flow rate of the pumping liquid, the satellite droplets merged with the subsequent mother droplets, and we were thus able to generate monodisperse Experiments on the generation of monodisperse deionized water droplets based on the piezoelectric method were conducted. Our observations demonstrated that satellite droplets were generated because the mother droplets detached from the liquid filament with surface waveform perturbations. As the frequency of the square-wave pulse was tuned to match the volume flow rate of the pumping liquid, the satellite droplets merged with the

droplets with a uniform size distribution. If the pulse frequency is relatively low, the

the generated droplets are chaotic. This results in a wide size distribution. Under different volume flow rates, the optimal ranges of pulse frequency required to generate monodisperse droplets are different. The size of the monodisperse droplets decreases with the increase in the pulse frequency, as the volume flow rate is given. The sizes of the monodisperse droplets generated follow bimodal and Gaussian distributions under

**Author Contributions:** Conceptualization, X.H. and S.L.; methodology, X.H. and Z.H.; experiments, Z.H., F.Y., and X.L.; writing, X.H. and S.L.; review and editing, S.P. and S.L.; supervision, J.X. All

**Funding:** This work was financially supported by the National Natural Science Foundation of China (51976050 and 41805108), the Fundamental Research Funds for the Provincial Universities of Zhejiang (GK209907299001-004), and the Scientific Research Projects of Zhejiang Administration for

authors have read and agreed to the published version of the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest

subsequent mother droplets, and we were thus able to generate monodisperse droplets with a uniform size distribution. If the pulse frequency is relatively low, the satellite droplets cannot merge with the subsequent mother droplets. It is also possible that the two adjacent mother droplets can coalesce if the pulse frequency is so high that the generated droplets are chaotic. This results in a wide size distribution. Under different volume flow rates, the optimal ranges of pulse frequency required to generate monodisperse droplets are different. The size of the monodisperse droplets decreases with the increase in the pulse frequency, as the volume flow rate is given. The sizes of the monodisperse droplets generated follow bimodal and Gaussian distributions under optimum operating conditions.

**Author Contributions:** Conceptualization, X.H. and S.L.; methodology, X.H. and Z.H.; experiments, Z.H., F.Y., and X.L.; writing, X.H. and S.L.; review and editing, S.P. and S.L.; supervision, J.X. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was financially supported by the National Natural Science Foundation of China (51976050 and 41805108), the Fundamental Research Funds for the Provincial Universities of Zhejiang (GK209907299001-004), and the Scientific Research Projects of Zhejiang Administration for Market Regulation (20200103).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Lizhong Huang , Jiayou Du and Zefei Zhu \***

School of Mechanical Engineering, Hangzhou Dianzi University, Hangzhou 310018, China; ricky@hdu.edu.cn (L.H.); abc@hdu.edu.cn (J.D.)

**\*** Correspondence: zzf\_3691@163.com; Tel.: +86-571-8691-9007

**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 ≤ *Re* ≤ 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 *Re*, while *k* has little impact. The quasi frequency of such oscillation has almost no business with *Re* and *k*. Moreover, *Re* 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.

**Keywords:** lattice Boltzmann method; inertial migration; Poiseuille flow; pulsatile velocity

**Citation:** 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

Academic Editor: Myung-Suk Chun

Received: 4 August 2021 Accepted: 3 September 2021 Published: 6 September 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **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 (*Re* = *UmH*/*ν*, where *U<sup>m</sup>* is the maximum inlet velocity, *H* is the channel width, *ν* is the fluid kinematic viscosity) and blockage ratio (*k* = *Dp*/*H*, donates the ratio of particle diameter *D<sup>p</sup>* and the channel width), whether or not the particles are neutrally buoyant [1–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

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 *Re* 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 *Re* 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–19], the same problem was studied, and the SS effect was reconstructed. Shao et al. [20] found the inner annulus for elevated *Re* by using the fictitious domain method. Abbas et al. [21] mentioned that the equilibrium position (TEP) depends exclusively on *Re* and *k*. Recently, inertial microfluidics can precisely separate particles with or without extra external force field by realizing SS effect [22–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 *Re* 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–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(\mathbf{x} + \mathbf{e}\_i \Delta t, t + \Delta t) = f\_i(\mathbf{x}, t) - \frac{1}{\tau} \left[ f\_i(\mathbf{x}, t) - f\_i^{(\neq q)}(\mathbf{x}, t) \right],\tag{1}$$

$$\begin{array}{ccccccccc}\hline f & & & f & & f & & f\\ \hline f & & & f & & f & & f & \\ \hline \end{array}$$

where *fi*(*x*, *t*) is the distribution function at space coordinate *x* = (*X*,*Y*) and time *t* in the *i*th direction; *f* (*eq*) *i* (*x*, *t*) is the correspond equilibrium distribution function; *τ* is the SRT; ∆*t* is the time step; the discrete velocities **e***<sup>i</sup>* 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. ሺ, ሻ ൌ ሺ, ሻ ሺሻ ሺ, ሻ Δ ൌ Δ Δ ⁄ Δ Δ ൌ Δ ൌ ൌ 1

– –

ሺ, ሻቃ,

**Figure 1.** D2Q9 Cartesian lattice and discrete velocities.

The equilibrium distribution function can be calculated by [44]:

$$f\_l^{(\mathbf{\hat{e}}\boldsymbol{\eta})}(\mathbf{x},t) = \omega\_l \rho\_f \left[ 1 + \frac{3\mathbf{e}\_l \cdot \boldsymbol{u}}{c^2} + \frac{4.5\left(\mathbf{e}\_l \cdot \boldsymbol{u}\right)^2}{c^4} - \frac{1.5\boldsymbol{u}^2}{c^2} \right],\tag{2}$$

 ൌ4 9⁄ ଵ~ସ ൌ1 9⁄ ହ~଼ ൌ 1 36 ⁄ where *ω<sup>i</sup>* is the weight factor with *ω*<sup>0</sup> = 4/9, *ω*1∼<sup>4</sup> = 1/9 and *ω*5∼<sup>8</sup> = 1/36, *ρ<sup>f</sup>* is the fluid density and *u* is fluid velocity which can be determined by:

$$\rho\_f = \sum f\_i(\mathbf{x}, t), \ \mathbf{u} = \frac{1}{\rho\_f} \sum \mathbf{e}\_i f\_i(\mathbf{x}, t). \tag{3}$$

– 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.

ሺሻ ൌ

⎩ ⎨

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

 <sup>ଷ</sup> <sup>ସ</sup> Take fluid node *A* as an example, after the streaming step, three unknown distribution functions (*f*3, *f*<sup>7</sup> and *f*<sup>4</sup> denoted by three blue arrows shown in Figure 2) need to be determined by applying improved bounce-back scheme. Therefore, the distribution function *f*<sup>7</sup> can be calculated by:

$$f\_{\mathcal{T}}(\boldsymbol{A}) = \begin{cases} \begin{array}{c} f\_{\mathcal{T}}(\boldsymbol{A}) = \begin{cases} \frac{q(1+2\eta)f\_{5}(\boldsymbol{S}) + \left(1-4\eta^{2}\right)f\_{5}(\boldsymbol{A}) - q(1-2\eta)f\_{5}(\boldsymbol{B}) - 2\omega\varepsilon\_{\mathcal{T}}\rho\_{f}\frac{\mathbf{e}\_{3}\cdot\mathbf{u}\_{D}}{\varepsilon\_{s}^{2}}, & q < \frac{1}{2}, \\\ \frac{1}{q(2\eta^{2}\boldsymbol{I})}f\_{5}(\boldsymbol{\Theta}) + \frac{2q-1}{q}f\_{7}(\boldsymbol{B})\boldsymbol{I} - \frac{2q-1}{2q+1}f\_{7}(\boldsymbol{C}) - \frac{2\omega\varepsilon\_{\mathcal{T}}\rho\_{f}}{q(2q+1)}\frac{\mathbf{e}\_{3}\cdot\mathbf{u}\_{D}}{\varepsilon\_{s}^{2}}, & q > \frac{1}{2}, \end{cases} \end{cases} \tag{4}$$

1 ሺ2 1ሻ ହሺሻ 2 െ 1 ሺሻ െ 2 െ 1 2 1 ሺሻ െ 2ହ ሺ2 1ሻ ହ ⋅ ௦ ଶ , ⩾ 1 2 , <sup>ହ</sup> ൌ || | ⁄ | <sup>௦</sup> ൌ ⁄√3 where *S* is the nearest solid node along **e**<sup>5</sup> direction; *B* and *C* are the two nearest fluid nodes along **e**<sup>7</sup> 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<sup>D</sup>* is the velocity of the boundary location *D*; *c<sup>s</sup>* = *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:

$$\begin{array}{l} F^{(h)}(\mathbf{x} + q\mathbf{e}\_{5}, t) = \mathbf{e}\_{5} [f\_{5}(\mathbf{x} + \mathbf{e}\_{5}, t) + f\_{7}(\mathbf{x}, t)], \\ T^{(h)}(\mathbf{x} + q\mathbf{e}\_{5}, t) = \begin{pmatrix} \mathbf{x} + q\mathbf{e}\_{5} - \mathbf{x}\_{p} \end{pmatrix} \times F^{(h)}(\mathbf{x}, t), \end{array} \tag{5}$$

 ሺሻሺ <sup>ହ</sup> where *x<sup>p</sup>* is the position of the particle.

, ሻ ൌ ൫ <sup>ହ</sup> െ ൯ൈሺሻሺ, ሻ, ∆ 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]:

ሺ௨ሻሺ, ሻ ൌ ൫ െ ൯ൈሺ௨ሻሺ, ሻ.

$$\begin{aligned} \mathbf{F}^{(c)}(\mathbf{x},t) &= \rho\_f(\mathbf{x},t)\boldsymbol{\mu}(\mathbf{x},t), \\ \mathbf{T}^{(c)}(\mathbf{x},t) &= (\mathbf{x} - \mathbf{x}\_p) \times \mathbf{F}^{(c)}(\mathbf{x},t), \\ \mathbf{F}^{(u)}(\mathbf{x},t) &= -\rho\_f(\mathbf{x},t)\boldsymbol{\mu}(\mathbf{x},t), \\ \mathbf{T}^{(u)}(\mathbf{x},t) &= \left(\mathbf{x} - \mathbf{x}\_p\right) \times \mathbf{F}^{(u)}(\mathbf{x},t). \end{aligned} \tag{6}$$

Moreover, in order to avoid unphysical overlapping between the particle and the channel wall, the extra lubrication force model is needed [49]:

$$\mathbf{F}^{(l)} = \begin{cases} 0, & h \ge h\_{\mathcal{C}} \\ -1.5\pi \rho\_f \nu \left[ D\_p \left( \frac{1}{h} - \frac{1}{h\_{\mathcal{C}}} \right) \right]^{1.5} \mathbf{U}\_{p\nu} & h < h\_{\mathcal{C}} \end{cases} \tag{7}$$
 
$$\mathbf{U}^{(l)} := \begin{array}{cccccccccc} & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \end{array} \tag{8}$$

where *ν* = *c* 2 *s* (*<sup>τ</sup>* − 0.5)∆*<sup>t</sup>* is the kinematic viscosity of the fluid; *<sup>D</sup><sup>p</sup>* is the diameter of the particle; *U<sup>p</sup>* is the particle velocity towards the wall; *h* is the minimum gap between the particle and the wall; *h<sup>c</sup>* = 1.5∆*x* is the cutoff distance whether to consider the lubrication force or not. െ1.5 ൬ 1 ℎ െ 1 ℎ ൰൨ଵ.ହ , ℎ ൏ ℎ , ൌ<sup>௦</sup> ଶ ሺ െ 0.5ሻΔ ℎ

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: ℎ ൌ 1.5Δ

$$\begin{aligned} \mathbf{a}\_p^{t+\Delta t} &= \left(\sum \mathbf{F}^{(h)} + \sum \mathbf{F}^{(c)} + \sum \mathbf{F}^{(u)} + \mathbf{F}^{(l)}\right) / m\_p, \\ \mathbf{a}\_p^{t+\Delta t} &= \mathbf{a}\_p^t + 0.5 \left(\mathbf{a}\_p^{t+\Delta t} + \mathbf{a}\_p^t\right) \Delta t, \\ \mathbf{x}\_p^{t+\Delta t} &= \mathbf{x}\_p^t + 0.5 \left(\mathbf{a}\_p^{t+\Delta t} + \mathbf{a}\_p^t\right) \Delta t, \\ \mathbf{a}\_p^{t+\Delta t} &= \left(\sum \mathbf{T}^{(h)} + \sum \mathbf{T}^{(c)} + \sum \mathbf{T}^{(u)}\right) / I\_{p\prime} \\ \mathbf{w}\_p^{t+\Delta t} &= \mathbf{w}\_p^t + 0.5 \left(\mathbf{a}\_p^{t+\Delta t} + \mathbf{a}\_p^t\right) \Delta t, \\ \mathbf{\theta}\_p^{t+\Delta t} &= \mathbf{\theta}\_p^t + 0.5 \left(\mathbf{w}\_p^{t+\Delta t} + \mathbf{w}\_p^t\right) \Delta t, \end{aligned} \tag{8}$$

where *ap*, *up*, *αp*, *wp*, *θp*, *mp*, and *I<sup>p</sup>* 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 *Um* is set at the left inlet boundary in the positive *X* direction, and the velocity in *Y* direction is zero. For *Um*, there are two cases: non-pulsatile or pulsatile. If pulsatile, *Um* 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*<sup>1</sup> − *t*<sup>0</sup> = 0.3 s, and in the diastolic period, *t*<sup>2</sup> − *t*<sup>1</sup> = 0.5 s, which means one cardiac cycle lasts 0.8 s. The unit conversion factor from lattice time to physical time is 10−<sup>5</sup> . At the right boundary, the normal derivative of the velocity is zero and the pressure is set to be *pout* = *ρ<sup>f</sup> c* 2 *s* [19]. No-slip boundary conditions are imposed at the top and bottom channel walls. ሾ , <sup>ଵ</sup> ሿ ሾ<sup>ଵ</sup> , <sup>ଶ</sup> ሿ <sup>ଵ</sup> െ ൌ 0.3 s <sup>ଶ</sup> െ <sup>ଵ</sup> ൌ 0.5 s 0.8 s 10ିହ ௨௧ ൌ <sup>௦</sup> ଶ

 ൈ ሾ<sup>௦</sup> , <sup>௦</sup> ሿ **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 [*Xs*, *Ys*] initially. The diameter of the particle is *Dp*, *Um* represents the maximum velocity. For pulsatile case, *Um* 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<sup>s</sup>* , *Ys*] 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* = *Dp*/*H* = 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, respectively. Meanwhile, *Re* = *UmH*/*ν* = 25, 50, 75, 100, 150, 200 is studied by changing the kinematic viscosity of the fluid. ሾ<sup>௦</sup> , <sup>௦</sup> ሿ 15 20 25 30 35 40 ൌ⁄ ൌ 0.15, 0.20, 0.25, 0.30, 0.35, 0.40 ൌ/ ൌ 25, 50, 75, 100, 150, 200 <sup>௦</sup> Δ

100 ሺ300 400

500 600 700 2000ሻ ൌ 500

500

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<sup>s</sup>* + ∆*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, *Re* is set to be 50, and the terminal particle Reynolds number *Re<sup>p</sup>* = *Ux*,*pDp*/*ν* is about 9.63, where *Ux*,*<sup>p</sup>* is the terminal particle velocity in *X* direction. It is very close to *Re<sup>p</sup>* 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 (*Ys*/*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]. ௫,/ ௫, ൌ 0.25 ൌ 0.35 ௦⁄ ൌ 0.20, 0.25, 0.35, 0.40, 0.45

ሺሻ ൌ 0.25 ሺሻ ൌ 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.

 ൌ 20, 40, 100, 200 ൌ 22 ൌ 200 ൌ 1000 ൌ 0.11 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<sup>p</sup>* = 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.

 **Figure 5.** Comparison of TEPs of the particle migrating in Poiseuille flow at different *Re*.

#### **4. Results and Discussions**

#### *4.1. Fluid and Particle Interaction*

 ൌ 0.25 ൌ 50 <sup>ᇱ</sup> ൌ ሺെ௨௧ሻ⁄௨௧ <sup>௫</sup> <sup>ᇱ</sup> ൌ ௫⁄ <sup>௬</sup> <sup>ᇱ</sup> ൌ ௬⁄ 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* − *pout*)/*pout*, the dimensionless fluid velocity in *X* direction *U*′ *<sup>x</sup>* = *Ux*/*U<sup>m</sup>* and in *Y* direction *U*′ *<sup>y</sup>* = *Uy*/*Um*, 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]. ൌ 0.25 ൌ 50 <sup>ᇱ</sup> ൌ ሺെ௨௧ሻ⁄௨௧ <sup>௫</sup> <sup>ᇱ</sup> ൌ ௫⁄ <sup>௬</sup> <sup>ᇱ</sup> ൌ ௬⁄

ሺሻ ᇱ ሺሻ <sup>௫</sup> ᇱ ሺሻ <sup>௬</sup> ᇱ ሺሻ ᇱ ሺሻ <sup>௫</sup> ᇱ ሺሻ <sup>௬</sup> ᇱ ൌ 0.25 ൌ 50 **Figure 6.** The contour view of (**a**) the dimensionless pressure *p* ′ , (**b**) the dimensionless fluid velocity in *X* direction *U*′ *<sup>x</sup>* and (**c**) in *Y* direction *U*′ *y* in stable state of particle migrating in non-pulsatile flow at *k* = 0.25 and *Re* = 50.

#### ൌ 0.25 ൌ 50 *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 *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. ൌ 0.25 ൌ 0.35 ൌ 50

 ൌ 0.25 ൌ 0.35 ൌ 50 **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 *Re* = 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.

ሺ <sup>ᇱ</sup> ൌ 10<sup>ହ</sup> ⁄ ሻ ൌ 0.25 ሺ, ሻ ൌ 50 ሺ, ሻ ሺ, ሻ ሺ <sup>ᇱ</sup> ൌ ሾ8.8, 9.6ሿሻ ∆ ∗ , ∆, <sup>∗</sup> , ∆ ∗ , ∆, <sup>∗</sup> , ∆ , ∆ ൌ 0.2 Figure 8a,c gives the particle center trajectory versus time *t* ′ = *t*/10<sup>5</sup> variation with *Re* at *k* = 0.25, and (*b*, *d*) donates different *k* at *Re* = 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 *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*.

<sup>∗</sup>

 ሺ, ሻ ൌ 0.25 ሺ, ሻ ൌ 50 , ∆, <sup>∗</sup> ሺሻ ሺሻ **Figure 8.** Time history of the particle center trajectory at different *Re* (**a**,**c**) at *k*= 0.25, and different *k* (**b**,**d**) at *Re* = 50. The variation of *δ*, ∆, *t* ∗ with *Re* (**e**) and *k* (**f**).

#### *4.3. Orientation*

 ൏ 100 ሺሻ 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*. 

 ሺ, ሻ ሺ, ሻ ሺ, ሻ ሺ, ሻ **Figure 9.** The particle orientation at different *Re* (**a**,**c**) and *k* (**b**,**d**). Meanwhile, (**a**,**b**) is in non-pulsatile and (**c**,**d**) in pulsatile flow.

#### *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 (*Ux*) is preferred, because the particle velocity curve is not center symmetry in *Y* direction.

 <sup>௫</sup> <sup>ᇱ</sup> ൌ ௫⁄ <sup>௫</sup> ᇱ <sup>௫</sup> ᇱ ൎ 1 ௨ ൌ ିଵହ.ଽଽ௧ᇲାଵସଶ.ଷଽ ൌ 15.99 For better illustration, the particle velocity in *X* direction is normalized by: *U*′ *<sup>x</sup>* = *Ux*/*Um*. The curves of *U*′ *x* in Figure 10 are much like a spring-mass system which is underdamped. Figure 10a shows, when *Re* is small, *U*′ *<sup>x</sup>* 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, *Aup* = *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 *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.

௫ <sup>ᇱ</sup> ൌ 0.25 ൌ ሺሻ 25, ሺሻ 50, ሺሻ 75, ሺሻ 100, ሺሻ 150, ሺሻ 200 **Figure 10.** Damping of *U*′ *x* (solid line in sky-blue color) of the particle at *k* = 0.25 and different *Re* = (**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.

௫ <sup>ᇱ</sup> ൌ 50 ൌ ሺሻ 0.15, ሺሻ 0.20, ሺሻ 0.25, ሺሻ 0.30, ሺሻ 0.35, ሺሻ 0.40 **Figure 11.** Damping of *U*′ *<sup>x</sup>* at *Re* = 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.

<sup>ᇱ</sup> ൌ 50 ൌ ሺሻ 0.15, ሺሻ 0.20, ሺሻ 0.25, ሺሻ 0.30, ሺሻ 0.35, ሺሻ 0.40

<sup>௫</sup>

<sup>௫</sup>

௫

ሺሻ 25, ሺሻ 50, ሺሻ 75, ሺሻ 100, ሺሻ 150, ሺሻ 200

ൌ 25 <sup>௫</sup>

ൌ 25

ൌ 25

ᇱ

 <sup>௫</sup> <sup>ᇱ</sup> ൌ 25 <sup>௫</sup> <sup>ᇱ</sup> ൌ 25 ൌ 25 <sup>௫</sup> ᇱ Figure 12a gives the quasi frequency (*fq*) of *U*′ *<sup>x</sup>* 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*′ *<sup>x</sup>* 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.

ൎ ൫2 ⁄ ൯

ᇱ

ൎ ൫2 ⁄ ൯

ᇱ

'

ൌ 0.83ି.ସ

ሺሻ ሺሻ **Figure 12.** The quasi frequency (**a**) and damping ratio (**b**) at different *Re* and *k*.

 As shown in Figure 12b, the damping ratio *ζ* ≈ *λ*/ 2*π fq* , 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*′ *<sup>x</sup>* 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.83*Re*−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 *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.

**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**

