The dynamics of a droplet in oscillatory and pulsating flows of a Newtonian fluid in a microchannel has been studied numerically. The effects of oscillation frequency, surface tension, and channel flow rate have been explored by simulating the drop within a microchannel. These types of flows introduce new equilibrium positions for the drop compared to steady flows with similar conditions. The simulation results are very sensitive to the grid resolution due to the unsteady behavior of the base flow. Therefore, a set of fine grids have been used in this study to capture the physics of this problem more accurately. However, these fine grids make the computations significantly expensive. Therefore, a multifidelity Gaussian processes method with two levels of fidelity has been used to predict the results of the remaining fine-grid simulations along with their uncertainties based on their correlations with those of the coarse-grid cases over a wide range of input parameters.

Transient dynamics of droplets and their lateral equilibrium positions in microchannels depend on several parameters including the rheology of the carrier fluid, the geometry of the channel, the channel flow rate, and the drop characteristics such as shape, size, interfacial tension, etc, which can be deployed for numerous biological and biomedical applications.1–6 

Most of the effective parameters underlying the physics of particle motion in microchannels and label-free approaches are inherent to the system, particle properties, and are out of human's complete control. Therefore, adding an extra parameter as a tool for more direct control over the dynamics of particles and better manipulating them can be very beneficial; especially, coming up with equilibrium positions between the channel center and the wall can be very useful for separation purposes.7 The other restriction in inertial microfluidics is the limitation of working with particle sizes of a few micrometers or larger only.8 This is because of the inverse relationship between the particle size and its traveling distance until focusing, which means that the separation of smaller particles requires longer channels.9 

Exploiting an oscillatory flow in the channel has been recommended as a solution to address the aforementioned issues.8,10 The oscillation frequency as the new extra parameter does not have the disadvantage of adversely affecting the biological properties of cells like most of the active methods.11 Besides, this type of flow as an alternative to traditional microfluidics has attracted researchers in recent years. Chaudhury et al.12 have observed a complicated trajectory of the droplet in the lateral direction suspended in an oscillatory flow as opposed to the smooth cross-stream migration under steady flow conditions. Pawłowska et al.13 have reported the use of hydrogel nanofilaments as an appropriate substitute for the long and deformable macromolecules and studied their dynamics in an oscillatory microchannel flow. They have shown that the final position of these particles fluctuates around the flow axis.14 Zhu et al.15 have illustrated that in a wall-bounded oscillating shear flow, the migration velocity of capsule decreases monotonically by increasing the frequency of applied shear and the distance of capsule from wall. Sarkar and Schowalter16 have investigated the dynamics of a viscoelastic drop in steady and oscillating extensional flows and have shown that although their deformation behaviors are naturally dissimilar, their maximum values are close in the long-time limit. Zhang et al.17 have found that the effect of capsule membrane viscosity is more significant in the oscillatory shear flow compared to that of the steady one, and the deformation of the capsule is influenced by both viscosity and elasticity and exhibits two modes. Lu et al.18 have performed a numerical and experimental study on an ultrasonic oscillatory air–water two-phase flow in a microchannel. They have observed highly unsteady behavior as the water and air interact with each other during the vibration cycles, which is drastically different from the steady flow in such microchannels.

The presented study investigates different aspects of the droplet migration in oscillatory and pulsating microchannel flows numerically and compares them with those in steady ones. The main challenge of performing this numerical work is using a sufficient grid resolution to resolve the underlying physics more accurately. Since the background flow is unsteady, the mesh has to be fine enough to capture the velocity gradients and lift forces acting on the drop correctly. Consequently, the transient dynamics of the drop and its final equilibrium position, which is the main goal of the presented study, are all affected by this choice. Therefore, after carrying out the mesh dependence study, a very fine grid has been selected to enable getting the realistic results we are interested in. However, running simulations with these grids can take a few months in some cases. Hence, it is impossible to get all the data in the wide range of input parameters we are looking for by just running all these simulations. An alternative approach is to produce some new data using a coarser grid resolution and implement the multifidelity Gaussian processes (MFGP) algorithm to predict the outputs of the finer mesh. Perdikaris et al.19 have proposed a probabilistic framework and a recursive version of MFGP. This method has been tested in several benchmark problems involving both synthetic and real multifidelity datasets such as the one employed by Babaee et al.20 The implementation of the recursive MFGP and the obtained results in this study are elaborated in Secs. II B and III C.

A single Newtonian droplet has been placed in a laminar flow of an incompressible Newtonian fluid in a rectangular microchannel with a square cross section. A schematic of the configuration is illustrated in Fig. 1(a). The density and viscosity ratios (η and λ, respectively) are set to 1 for most of the simulations. The front-tracking method21–28 is used to update the interface position. In this method, the main governing equations are solved in the fixed Eulerian grid, and the obtained information is used to update the properties across the droplet surface containing thousands of moving Lagrangian elements. The governing equations to be solved in the entire computational domain are the following:

·u=0,
(1)
(ρu)t+·(ρuu)=p+·τ+γκδ(xxi)ndA,
(2)

where ρ is the density of the fluid, p represents the pressure, u is the velocity vector, t is the time, τ=μ(u+uT) is the stress tensor in which μ is the fluid viscosity, κ is the curvature at the interface, γ is the interfacial tension, δ is the Dirac delta function, x is an arbitrary location in the whole computational domain, xi is such position on the drop interface, and n is the unit normal vector to a point on the interface. The given delta function is defined as

δ(x)=D̃(x)D̃(y)D̃(z),
(3)
D̃(x)=14Δ(1+cos(π2Δ(x))),|x|2Δ,
(4)

where Δ is the constant Eulerian grid size.

FIG. 1.

(a) Schematic of the problem setup and (b) an example of the droplet discretization.

FIG. 1.

(a) Schematic of the problem setup and (b) an example of the droplet discretization.

Close modal

To generate the oscillatory flow, a cosine wave of pressure gradient with a constant amplitude [in the form of P0cos(ωt)] is applied along the channel (x direction) to change the direction of the flow symmetrically. For the pulsating case, this pressure gradient has the shape of P0(a+bcos(ωt)), where a and b are the weights of the steady and oscillatory components, respectively, and a+b=1. The periodic boundary condition is applied in the x direction, and the no-slip condition is applied on the walls in the y and z directions. W and U0 (channel centerline velocity corresponding to the steady case) are used as the characteristic length and velocity, respectively. In other words, x*=xW,u*=uU0,t*=tWU0,P*=PμU0W,T*=TWU0 (where T is the period), and ω*=2πT*. Three dimensionless parameters describe and affect the motion of the drop: (i) Reynolds number Re=ρU02Wμ, expressing the ratio between inertial forces to viscous ones; (ii) Capillary number Ca=μU0γ, which denotes the ratio of viscous stress to the interfacial tension, where high Ca corresponds to a highly deformable drop; and (iii) the dimensionless oscillation frequency (ω*). The effect of the frequency can also be embedded in the Womersley number Wo=(ωW2ν)12, where ν is the kinematic viscosity of the fluid. This number compares the transient inertial effects to viscous forces.29 The blockage ratio of the drop (aW) is constant and equals to 0.3 for most of the cases studied here. The drop is assumed to have a spherical initial shape and is released at yW=0.55 and zW=1.07. The initial location of the drop is arbitrary since it does not alter its equilibrium position.9,12,30,31 The axes of symmetry have been avoided. A fine grid of 196×114×114 (high fidelity level) and a coarse grid of 128×76×76 (low fidelity level) in the x, y, and z directions, and with 29 578 and 13 038 triangular elements for the discretization of the droplet [as shown in Fig. 1(b)], respectively, is used for the simulations having Re = 10. The simulations with higher Re numbers require even a finer grid resolution (Eulerian grid of 256×152×152 with 48 050 Lagrangian elements). Since this adds another level of fidelity to the problem, we have omitted the effect of Re in the second part of the paper for more simplicity. Therefore, the input parameters that affect the output of interest, which is the distance of the equilibrium position of the drop from the channel center, are reduced to the Capillary number (Ca) and frequency (ω*).

Gaussian processes (GP) is an example of a continuous stochastic process in which, a multidimensional Gaussian distribution is assigned to the function values at the input points as the prior knowledge. In other words, the value of the function that we are aiming to predict follows a Gaussian probability density at each point:

f1:n|x1:nN(m(x1:n),K(x1:n,x1:n)),m(x1:n)=(m(x1),,m(xn)),K(x1:n,x1:n)=(k(x1,x1)k(x1,xn)k(xn,x1)k(xn,xn)),

where m is the mean vector, and K is the covariance matrix whose elements denote the correlations between the input points. These correlations are quantified by a covariance (kernel) function that encodes our prior belief and knowledge about the targeted function like its continuity, differentiability, etc. The choice of this function is very important for a good prediction; otherwise, we need to have sufficient amount of data to compensate for this lack of knowledge. One of the most typical examples for this function is the squared exponential (SE) kernel given below:

k(x,x)=vexp{12i=1d(xixi)2i2}.

In which, v is known as the signal strength, and i is the length scale of the ith input dimension of the GP. In this study, we also use the SE function as we believe the function of interest is continuous and infinitely differentiable.

If we have n observations (training set), consisting of x1:n=(x1,,xn) and y1:n=(y1,,yn), and n* test points (x1:n**) that we would like to predict the function values at (f1:n**), we use the following joint probability, according to the definition of GP, as our prior,

p(f1:n,f1:n**|x1:n,x1:n**)=N((f1:nf1:n**)|(m(x1:n)m(x1:n**)),(K(x1:n,x1:n)K(x1:n,x1:n**)K(x1:n**,x1:n)K(x1:n**,x1:n**))).

We then model our likelihood as p(y1:n|f1:n)N(f1:n,σ2In), which is also a Gaussian distribution centered at the known function values and has a variance of σ associated with the measurements noise. Since our data are from simulations, we can assume the observations are not noisy (σ = 0 and y1:n=f1:n). Thus, after performing the Bayes' rule, where D=(x1:n,y1:n) is the observed data, we get the posterior distribution for f* as below,

p(f1:n**|x1:n**,D)=N(f1:n**|mn(x1:n**),Kn(x1:n**,x1:n**)),

where the posterior mean function is

mn(x)=m(x)+k(x,x1:n)(K(x1:n,x1:n))1(y1:nm(x1:n)),

and the posterior covariance function is

kn(x,x)=k(x,x)k(x,x1:n)(K(x1:n,x1:n))1kT(x,x1:n).

With k(x,x1:n)=(k(x,x1),,k(x,xn)) being the cross-covariance vector.

Multi-fidelity modeling enables accurate inference of quantities of interest by synergistically combining realizations of low-cost/low-fidelity models with a small set of high-fidelity observations. This is particularly effective when the low and high-fidelity models exhibit strong correlations and can lead to significant computational gains over approaches that solely rely on high-fidelity models.19 The MFGP method introduced by Perdikaris et al.19 is fundamentally very similar to that of the Kennedy and O'Hagan,32 but it is also capable of capturing the more complex nonlinear cross-correlations between the data. Assuming that we have s fidelity levels (which is 2 in our case), based on this method, the GP has to be done s times recursively. In other words, the GP is done in its standard, usual way on the first level, and for the next levels, the outputs of the previous level t – 1 are the inputs for the level t. This is done in the way described below:19 

ft(x)=gt(x,f*t1(x)),

where f*t1(x) is the posterior distribution of the previous level t – 1 evaluated at the x, being the inputs of the current level t. The unknown function gt follows the description of GP:

gtGP(ft|0,kt((x,f*t1(x)),(x,f*t1(x));θt)).

According to Le Gratiet and Garnier,33 this scheme has the same posterior distribution predicted by the fully coupled scheme of Kennedy and O'Hagan.32 This procedure implies two conditions: (i) The training sets need to have a nested structure (i.e., xtxt1) (ii) The Markov property

cov{ft(x),ft1(x)|ft1(x)}=0,xx,

which translates into assuming that given the nearest point ft1(x), we can learn nothing more about ft(x) from any other model output ft1(x), for any xx.

Since f*t1 and x belong to inherently two different spaces, a more structured kernel function, coupling the elements from the same space together, is the following:19 

ktg=ktρ(x,x;θtρ)·ktf(f*t1(x),f*t1(x);θtf)+ktδ(x,x;θtδ),

where ktρ,ktf, and ktδ are valid covariance functions, and θtρ,θtf, and θtδ are their hyperparameters, respectively. For our application, we have chosen the SE kernel function

k(x,x;θt)=σt2exp{12i=1dWi,t(xixi)2},

where σt2 is a variance parameter, and {Wi,t}i=1d are the Automatic Relevance Determination (ARD) weights corresponding to the fidelity level t.

As previously discussed, the posterior distribution of the first level is obtained by the normal GP. Therefore, it is Gaussian. However, this is not generally the case for the subsequent levels, except for the case of t = 2. Therefore, for all the cases with t2, we have to perform predictions given uncertain inputs, where the uncertainty is propagated along each recursive step. Thus, the posterior distribution is given by19 

p(f*t(x*)):=p(ft(x*,f*t1(x*))|f*t1,x*,yt,xt)=p(ft(x*,f*t1(x*))|yt,xt,x*)p(f*t1(x*))dx*.

The predictive mean and variance of all posteriors p(f*t(x*)),t2 are calculated using Monte Carlo integration of this equation. We can then sample from this distribution.

Dynamics of the droplet in an oscillatory flow in the microchannel has been studied, and the effects of Re, Ca, and ω* have been investigated. Re ranges between 10 and 100, Ca ranges between 0.09 and 10, and ω* values are chosen such that for a channel with a square cross section of 100μm and water as the working fluid at room temperature, the frequency values range between 2 and 1600 Hz (or a corresponding Wo number between 0.2 and 6.3). Although a maximum frequency of 200 Hz is mostly reported in the literature,29 recent works have claimed of generating frequencies of around 1 KHz.34 The equilibrium position of the drop is a result of the competition between the lateral lift forces acting on it, including the wall effect and the deformation-induced lift forces, both acting toward the channel center, and the shear gradient force acting toward the wall.35–38 Magnus and Saffman lift forces are often very small compared to the other mentioned components and can be neglected.36,38–40 The boundary wall causes the particle to have rotational and translational velocities different from those of the adjacent fluid, which is caused by an uneven distribution of vorticities around the particle.41,42 This induces a higher pressure in the gap between the particle and the wall, which repels the particle away from the wall.36 The existing curvature of the fluid velocity profile makes the magnitude of the velocity of the fluid on the wall side much higher than the channel center side from the particle frame of reference. This inequality causes a low pressure on the wall side leading to a shear gradient lift force that pushes the particle toward the wall.36 Following the analytical results of Chan and Leal,43 the deformability-induced lift force for droplets or bubbles that have a distance higher than their diameter from the wall, which is the case in our simulations, is given by44 

FL,deformation=CapμVavga(aW)2dWf(λ),
(5)
f(λ)=128π(λ+1)3{11λ+10140(3λ2λ+8)3(19λ+16)14(3λ+2)(2λ2+λ1)},
(6)

where Cap=μU0γaW is the drop Capillary number, Vavg is the average velocity of the carrier fluid across the channel, d is the distance of the drop from the channel center, and λ is the viscosity ratio between the inner and outer fluids.

The code has been validated by comparing the drop deformations at Ca = 0.2 and different Deborah numbers with those of the Aggarwal and Sarkar.45 The results are in good agreement with a maximum error of 0.72%. To validate the inertial effects, the focal points of the drop at Re = 8.25, Ca = 0.18, aW=0.2 and Re = 21, Ca = 0.14, aW=0.3 have been compared with those presented by Marson et al.7 Our obtained focal points lie within their corresponding uncertainty bands. Furthermore, we have shown that the numerical results are independent of the distance between 2 consecutive drops in an infinite domain in the flow direction. This has been done by comparing the drop trajectory at Re = 10, Ca = 1, and ω*=0.1 for three different channel lengths of 4, 6, and 8 W in our simulation setup. The maximum difference between the drop trajectories for L=4 and L=6 is 0.0003W, and the one between those of L=4 and L=8 is 0.0005W. The results have also been shown to be grid independent, by comparing the equilibrium positions for the case of Re = 10, Ca = 1, and ω*=0.1 with two different grids of 196×114×114 and 256×152×152. The difference between their focal distances from the center is 0.0009W.

Figure 2 shows the distance of the droplet focal point from the channel center (d*) at different values of Wo, Ca, and Re. The drop focal point in the steady flow moves toward the center by increasing the Ca due to the increase in the deformation force9,30,46 and shifts toward the wall as Re increases because of the improvement in the strength of the shear gradient force.46–48 Since W and ν are constant in this work, ω* and Wo are directly related to each other. Hence, we can use them interchangeably. The droplet travels at locations far from the wall; so the wall lift can be neglected in our study.36,41 Both deformation and shear gradient lift forces, as the remaining active forces, depend on Vavg and d* that vary as the simulations proceed. The dependence of deformation force on these 2 parameters is apparent from Eq. (5). The parameter d* and the flow velocity at the drop location affect the magnitude of the difference between the velocities on the wall and center sides from the drop frame of reference. Hence, both d* and Vavg determine the magnitude of the shear gradient force. For non-steady flows, including oscillatory and pulsating ones, Vavg is time-dependent. The average of this Vavg in each half of a corresponding periodic cycle decreases as the ω* (or Wo) increases. Consequently, the averages of both forces in one periodic cycle change by changing the Wo value keeping other parameters fixed, leading to different equilibrium positions as we can see in Fig. 2. A complete explanation of the relationship between the focal point and parameters like ω*, Ca, and Re can be found in our previous work.10 According to this figure, the focal point is closest to the channel center at the highest Wo except for Re = 10 and Ca = 1 and Re = 10 and Ca = 1.67. This can be explained based on the shape of the flow velocity profile elaborated below.

FIG. 2.

Distance of the equilibrium position of the droplet from the channel center as a function of Wo for different cases.

FIG. 2.

Distance of the equilibrium position of the droplet from the channel center as a function of Wo for different cases.

Close modal

Stokes number in the oscillatory and pulsating flows is defined as St=Wων. O'Brien49 has solved these types of flows in rectangular channels analytically and quantified the shape of the velocity profile by calculating the ratio between the velocity at the center and its average across the cross section. Based on the values in Table I of that paper, the profile maintains its parabolic shape up to St = 3 for a square channel. Above this approximate value, the profile starts to become more like a flat, plug-like profile.49Figure 3 illustrates the shapes of the dimensionless averaged velocity along the flow direction near the drop focal points at St = 3 (Wo = 3.15) and St = 6 (Wo = 6.3). These shapes are consistent with the findings of Refs. 49 and 50. At Wo = 3.15 (ω*=2) the velocity shape is still parabolic and is similar to those of the lower Wo values. However, this shape changes to plug like at Wo = 6.3 (ω*=8). The average of Vavg at this Wo is the lowest because it has the highest frequency among others. Therefore, the value of deformation force on average is very small according to Eq. (5). However, due to the shape of the velocity profile, the relative flow velocity from the drop frame of reference is very small near the center and very large near the wall. Thus, there is a strong shear gradient force although the average of Vavg is small. Consequently, the focal point at Wo = 6.3 does not follow the trend observed for the lower Wo numbers and is pushed toward the wall.

FIG. 3.

Normalized flow direction averaged velocity profiles on the xy and xz planes at Re = 10, Ca = 1.67, and (a) ω*=2(Wo=3.15,St=3), (b) ω*=8(Wo=6.3,St=6), and (c) their shape comparison on the xy plane with the same constant z.

FIG. 3.

Normalized flow direction averaged velocity profiles on the xy and xz planes at Re = 10, Ca = 1.67, and (a) ω*=2(Wo=3.15,St=3), (b) ω*=8(Wo=6.3,St=6), and (c) their shape comparison on the xy plane with the same constant z.

Close modal

Furthermore, by taking a closer look at Fig. 3(c), we can see that the velocity has an opposite sign near the wall. Due to the present hysteresis in this type of flow, there is a lag in the response of fluid to the change of flow direction.12,51,52

The Taylor deformation parameter of the drop is defined as

D=LBL+B.
(7)

In which, L is the principal major axis, and B is the principal minor axis of an equivalent ellipsoidal particle. The parameter D is oscillatory for all the flows in this study except for the steady one. This parameter is proportional to the shear rate, which depends on the flow velocity. As a result, the average of D is lower at higher frequencies.53 This is reflected in Fig. 4 by visualizing the average of D in the corresponding periodic cycle as a function of time. This trend also implies that the amount of oscillations in the deformation is lower at higher frequencies since the minimum deformation in each cycle is zero. Moreover, it is apparent that the case with ω*=8 has a deformation of close to 0, which confirms that it has a very low Vavg. In fact, the droplet in this case travels about only 0.02W along the flow direction and remains almost spherical though the Ca = 1 corresponds to a very deformable drop. Nevertheless, it migrates around 0.08W from the initial location to its focal point.

FIG. 4.

Averaged deformation parameter vs averaged dimensionless time at Re = 10, Ca = 1, and different frequencies.

FIG. 4.

Averaged deformation parameter vs averaged dimensionless time at Re = 10, Ca = 1, and different frequencies.

Close modal

Previous works have reported that drops in the flow regimes of high Ca are elongated significantly leading to their breakup.7,9,30 This is the reason for the absence of any data for Ca = 10 and Re = 10 for steady and lower frequency flows in Fig. 2. The drop undergoes a very large deformation in these cases, which is not what it experiences at higher frequencies. A similar argument holds for Re = 37.8 and Ca = 1.67 since the drop is able to deform more easily at higher Re.54 

As it was discussed previously, the velocity in the presented type of flows is oscillatory. Consequently, all the active forces also experience fluctuations leading to oscillations in the trajectory of the drop except for the steady flows.12 This can be seen in Fig. 5 where the whole lateral migration patterns at Re = 10, Ca = 1, and all frequencies, ω*=1, and ω*=0.01, showing the oscillations in the trajectories are illustrated in Figs. 5(a)–5(c), respectively, and the amplitudes of oscillations after focusing (A*) for different cases are depicted in Fig. 5(d). In these figures, d* is the dimensionless, time-dependent distance of the drop from the channel center, and the insets of Figs. 5(b) and 5(c) show the trajectory at the last 2 periodic cycles. The plots depicted in the latter 2 figures qualitatively agree with those of a previous study and become more like a helical, spiral pathway as the frequency increases.12 The minimum velocity in each periodic cycle is zero, occurring when the flow direction changes. The higher average velocity at lower frequencies implies a higher maximum velocity in the corresponding cycle. The higher the difference between the maximum and minimum velocities, the higher is the oscillations amplitude in the forces and in the trajectory. Therefore, similar to the aforementioned discussion of deformation oscillations, the parameter A* decreases as ω* or Wo increases [Fig. 5(d)]. This is also the case for the oscillations amplitude along the flow direction [the comparison between Figs. 5(b) and 5(c)]. The existence of these oscillations at all Wo values is noteworthy.12 It is vital to mention that all the values reported in Fig. 2 are the averages of d* in the last periodic cycle.

FIG. 5.

(a) Migration patterns at Re = 10, Ca = 1, and various frequencies, (b) drop trajectory vs the flow direction at Re = 10, Ca = 1, and ω*=1, and (c) ω*=0.01, and (d) amplitude of oscillations around the equilibrium point after focusing for different cases.

FIG. 5.

(a) Migration patterns at Re = 10, Ca = 1, and various frequencies, (b) drop trajectory vs the flow direction at Re = 10, Ca = 1, and ω*=1, and (c) ω*=0.01, and (d) amplitude of oscillations around the equilibrium point after focusing for different cases.

Close modal

Focusing time can be considered as an important factor in the design and performance of the microfluidic system. Nevertheless, it strongly depends on where the particle is initially released. Therefore, the average migration velocity is a better parameter for a more meaningful comparison under different circumstances. The average migration velocity (v*) for different cases is shown in Fig. 6. This velocity is computed by calculating the distance between the initial and equilibrium positions and dividing it by the focusing time. The corresponding focusing time is determined when the trajectory reaches within 0.015W of the focal point. This figure expresses that the average migration velocity decreases as Wo increase. This pattern is also observed in the average velocity along the flow direction. Besides, the average migration velocity in the steady flows increases by increasing the Ca, which is in agreement with the findings of Alghalibi et al.55 

FIG. 6.

The average migration velocity vs Wo.

FIG. 6.

The average migration velocity vs Wo.

Close modal

Viscosity ratio in the range of 1λ<13.2 has a weak effect on the drop migration.7,30 The effect of density ratio (η) on the drop focal point is even less.47 Therefore, we have limited our study on the effect of these two parameters on the drop migration only to one case as shown in Fig. 7. The drop with a λ higher than 1 focuses closer to the wall for the Ca we are studying here.7,30,47 We can also see that the change in the η and λ does not change the distance between the focal points in the steady and oscillatory flows significantly.

FIG. 7.

The effect of density and viscosity ratios on the distance of the drop focal point from the channel center at Re = 10 and Ca = 1.67.

FIG. 7.

The effect of density and viscosity ratios on the distance of the drop focal point from the channel center at Re = 10 and Ca = 1.67.

Close modal

Changing the steady flow to oscillatory type helps achieve different focal points, which can have potential applications in cell sorting and separation as a great advantage. However, the oscillatory flow has a zero net flow rate leading to a lower throughput compared to that of traditional steady flows,8 which can be counted as a disadvantage considering many microfluidic applications depending on high-throughput systems. One solution to fix this issue is to combine both steady and oscillatory parts in the pressure gradient and make the flow regime to be pulsating.29 A pulsating flow has the advantages of having a non-zero net flow rate as well as introducing a new equilibrium point for the drop. The latter occurs because the drop dynamics in the pulsating flow is very similar to the one in a pure oscillatory flow with an equivalent frequency between 0 and that of the oscillatory portion of the pulsating pressure gradient. This claim is valid since the weights of steady and oscillatory portions add up to 1. We can see this equivalent frequency in Fig. 8(a) and the inset of Fig. 8(c) where the trajectories in the pulsating flows have frequencies of almost half of those of the oscillatory ones with the same frequencies. This feature also enables the existence of cases with high Ca or Re and low ω* that are absent in the figures of Sec. III A. This is because the equivalent frequency of the pulsating flow makes the average deformation lower compared to the pure steady flow or oscillatory flow with a lower frequency. It is crucial to mention that the highest weights for the steady portions of the pulsating cases depicted in Fig. 8 denote flows in which the drop experiences the highest feasible deformation without breaking up or being significantly elongated. This has been done to show the highest possible changes in the focal points.

FIG. 8.

Evolution of the oscillatory and pulsating dynamics of the drop at Re = 37.8 and Ca = 1.67 by illustrating its (a) trajectory, (b) deformation, and at Re = 10 and Ca = 10 by visualizing the corresponding (c) trajectory, and (d) deformation.

FIG. 8.

Evolution of the oscillatory and pulsating dynamics of the drop at Re = 37.8 and Ca = 1.67 by illustrating its (a) trajectory, (b) deformation, and at Re = 10 and Ca = 10 by visualizing the corresponding (c) trajectory, and (d) deformation.

Close modal

In addition, the directionality of the focal points and averages of deformation values obey the expected trend at each combination of Re and Ca in Fig. 8. For instance, the pulsating flow at Re = 37.8 and Ca = 1.67 has an equivalent frequency less than 0.1. Therefore, we can see a focal point closer to the channel center and a higher average deformation according to Figs. 8(a) and 8(b), respectively. Similarly, we can observe a focal point closer to the center at Re = 10 and Ca = 10 for the pulsating case with the larger steady portion. This is because this case has an equivalent frequency between 0.1 and 0.5 and less than the one with the lower steady portion. Figure 8(c) reflects this pattern. The deformation behavior in Fig. 8(d) is also qualitatively similar to that shown in Fig. 8(b). It is also momentous to pay attention to the 0.012W difference between the focal points of the pulsating cases at Re = 10 and Ca = 10 although the difference between the weights of their steady portions is only 0.01.

Tables I and II quantify the average migration velocities for the flow regimes discussed above. The observed trend in these tables is consistent with the information provided in the discussion of Fig. 6. The pulsating case at Re = 37.8 and Ca = 1.67 has the lowest equivalent frequency. Hence, it has the highest average migration velocity among others. Similarly, the v* at Re = 10 and Ca = 10 decreases as the equivalent frequency increases.

TABLE I.

Average migration velocities for different flow regimes at Re = 37.8 and Ca = 1.67.

Pressure gradient formAverage migration velocity
P0(0.93+0.07cos(0.1t)) 0.00 191 
P0cos(0.1t) 0.00 141 
P0cos(0.5t) 0.00 012 
Pressure gradient formAverage migration velocity
P0(0.93+0.07cos(0.1t)) 0.00 191 
P0cos(0.1t) 0.00 141 
P0cos(0.5t) 0.00 012 
TABLE II.

Average migration velocities for different flow regimes at Re = 10 and Ca = 10.

Pressure gradient formAverage migration velocity
P0cos(0.1t) 0.000 389 
P0(0.17+0.83cos(0.5t)) 0.000 297 
P0(0.07+0.93cos(0.5t)) 0.000 213 
P0cos(0.5t) 0.000 205 
Pressure gradient formAverage migration velocity
P0cos(0.1t) 0.000 389 
P0(0.17+0.83cos(0.5t)) 0.000 297 
P0(0.07+0.93cos(0.5t)) 0.000 213 
P0cos(0.5t) 0.000 205 

The effect of droplet size on its focal point in the oscillatory and pulsating flows at Re = 10 and Ca = 10 is summarized in Table III. We observe that reducing the drop size pushes its equilibrium location toward the wall.30,47,48,56–58 Furthermore, this size reduction appears to enhance the change in the focal point at different values of equivalent frequency.

TABLE III.

The effect of drop size on its focal point at Re = 10 and Ca = 10.

Pressure gradient formEquilibrium distance from center
aW=0.3aW=0.2
P0cos(0.1t) 0.138 0.193 
P0(0.17+0.83cos(0.5t)) 0.166 0.241 
P0cos(0.5t) 0.183 0.266 
Pressure gradient formEquilibrium distance from center
aW=0.3aW=0.2
P0cos(0.1t) 0.138 0.193 
P0(0.17+0.83cos(0.5t)) 0.166 0.241 
P0cos(0.5t) 0.183 0.266 

In this section, we evaluate the MFGP performance by having a dataset consisting of 29 low-fidelity and 22 high-fidelity observations. Figure 9 shows all these observations together. It can be seen that the required nested structure, as mentioned earlier, is satisfied. In other words, for any data point in the high-fidelity level, there is a corresponding point in the low-fidelity level. Figure 10 illustrates the predictions on the distance of the drop equilibrium position from the channel center over the range of 0–1 for the frequency and 0.09–1.67 for the Capillary number. The algorithm is trained over the entire available data. A few sanity checks have been done to ensure that these results make sense. First, Fig. 10(a) denotes that at any fixed value of the Capillary number, there is a global optimum frequency for which the equilibrium distance from the center (the z axis) is an extremum. This is also previously observed in the simulations for the input values in the dataset (Fig. 2). Second, the predictions made by the high-fidelity response are slightly higher than those of the low-fidelity response, which is also compatible with the simulations outcomes. Finally, we know from the underlying physics that in the steady flow (at a frequency of 0), the distance decreases by increasing the Capillary number.9,30 This can be seen in both figures of the mean low and high responses [Figs. 10(b) and 10(c)] as well as being confirmed quantitatively in the code. Figure 10(d) illustrates a relatively low variance for the high-fidelity response, being our main goal, in the entire domain. The red-colored region in this plot corresponds to places where the density of data points is less, and hence, the predicted outputs have more uncertainty.

FIG. 9.

Observations of the high and low fidelity levels.

FIG. 9.

Observations of the high and low fidelity levels.

Close modal
FIG. 10.

Prediction of the MFGP algorithm by illustrating the (a) surface plot of the high-fidelity response colored by red and low-fidelity response colored by green, (b) mean response of the low-fidelity level, (c) mean response of the high-fidelity level, and (d) variance of the high-fidelity response.

FIG. 10.

Prediction of the MFGP algorithm by illustrating the (a) surface plot of the high-fidelity response colored by red and low-fidelity response colored by green, (b) mean response of the low-fidelity level, (c) mean response of the high-fidelity level, and (d) variance of the high-fidelity response.

Close modal

Since there is no analytical solution to compare the predictions with, we decided to split the high-fidelity data into training and test sets in this section. This helps us evaluate the performance of the implemented MFGP. We assign 15 training and 7 test points at the high-fidelity level randomly. Then, we train the MFGP on the whole low-fidelity data (29 points) and only the training high-fidelity data, and evaluate the predictive responses on the high-fidelity test points, since those are our main targets. We do this entire procedure 500 times to eliminate the dependence of the results on the test points choice. This especially helps us examine the algorithm performance at the regions with less amount of data. After this procedure, the average of the mean squared error (MSE) was 0.00015 and the average of the R2 score was 0.9858. The successful reproduction of these results is also checked. Figure 11 shows this evaluation at the last (500th) test set. Figure 11(a) expresses the predictions of d* along with their 95% credible intervals and denotes that the observed data lies within the corresponding uncertainty bars. In this figure, the x axis label (i) denotes the index of each test point. Figure 11(b) visualizes the true or known value of d* vs its predicted value at each test point. This figure illustrates the strength of the model as the plotted points are very close to the line of y = x, which is an indication of the agreement between predictions and observations. The exact and predicted correlations between the high and low responses of d* are also very close to each other according to Fig. 11(c). All of these are evidence for the strong and successful performance of the implemented MFGP algorithm.

FIG. 11.

Evaluation of the MFGP performance on a random test set by illustrating the (a) observations and predictions with their uncertainties, (b) predictions vs observations, and (c) exact vs predicted correlations.

FIG. 11.

Evaluation of the MFGP performance on a random test set by illustrating the (a) observations and predictions with their uncertainties, (b) predictions vs observations, and (c) exact vs predicted correlations.

Close modal

Determination and control of the particle's equilibrium position in the microchannels are extremely crucial as it can help in a variety of microfluidics applications. This importance, as well as the need to overcome the issue of designing impractically long channels to work with sub-micron particles, led us to do some simulations to capture the dynamics of a single droplet suspended in an oscillatory flow within the channel. The drawback of the zero net flow rate and smaller throughputs of the oscillatory flow compared to steady one has been addressed by combining these two to create a pulsating flow regime. Both oscillatory and pulsating flows bring new drop equilibrium locations to the system. Pulsating flows also enable the presence of droplets at high Ca or Re that could breakup in the steady or a very low-frequency oscillatory regime. Moreover, fluctuations in the trajectory of the drop have been observed. It has been shown that the amplitude of these oscillations, the average of the oscillatory deformation, and the average migration velocity all decrease by increasing the frequency. The dependence of the drop focal point on the shape of the velocity profile has been investigated as well. It has been explored that this equilibrium position moves toward the wall in a plug-like profile, which is the case in very high Wo numbers. Due to the significant cost of these simulations, a recursive version of the Multi Fidelity Gaussian processes has been used to replace the numerous high-fidelity simulations that cannot be afforded numerically. The MFGP algorithm is used to predict the equilibrium distance of the drop from the channel center for a given range of the interplaying input parameters, namely the Capillary number and frequency, assuming a constant Reynolds number. In addition, its performance was evaluated by randomly shuffling the high-fidelity data 500 times and assigning 31.8% of it as the test set for an accurate quantitative comparison each time. The algorithm outputs high statistical scores, which is an indication of its reasonably accurate performance.

This work is supported by the National Science Foundation (Grant No. CBET-1705371) and by the USDA National Institute of Food and Agriculture (Hatch Project No. 1017342).

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

1.
D. R.
Gossett
,
W. M.
Weaver
,
A. J.
Mach
,
S. C.
Hur
,
H. T. K.
Tse
,
W.
Lee
,
H.
Amini
, and
D. D.
Carlo
, “
Label-free cell separation and sorting in microfluidic systems
,”
Anal. Bioanal. Chem.
397
,
3249
(
2010
).
2.
M.
Toner
and
D.
Irimia
, “
Blood-on-a-chip
,”
Annu. Rev. Biomed. Eng.
7
,
77
(
2005
).
3.
D. R.
Gossett
,
H. T. K.
Tse
,
J. S.
Dudani
,
K.
Goda
,
T. A.
Woods
,
S. W.
Graves
, and
D. D.
Carlo
, “
Inertial manipulation and transfer of microparticles across laminar fluid streams
,”
Small
8
,
2757
(
2012
).
4.
A.
Karimi
,
S.
Yazdi
, and
A.
Ardekani
, “
Hydrodynamic mechanisms of cell and particle trapping in microfluidics
,”
Biomicrofluidics
7
,
021501
(
2013
).
5.
D. D.
Carlo
,
D.
Irimia
,
R. G.
Tompkins
, and
M.
Toner
, “
Continuous inertial focusing, ordering, and separation of particles in microchannels
,”
Proc. Natl. Acad. Sci. U. S. A.
104
,
18892
(
2007
).
6.
H.
Zhao
,
E. S.
Shaqfeh
, and
V.
Narsimhan
, “
Shear-induced particle migration and margination in a cellular suspension
,”
Phys. Fluids
24
,
011902
(
2012
).
7.
R. L.
Marson
,
Y.
Huang
,
M.
Huang
,
T.
Fu
, and
R. G.
Larson
, “
Inertio-capillary cross-streamline drift of droplets in Poiseuille flow using dissipative particle dynamics simulations
,”
Soft Matter
14
,
2267
(
2018
).
8.
B. R.
Mutlu
,
J. F.
Edd
, and
M.
Toner
, “
Oscillatory inertial focusing in infinite microchannels
,”
Proc. Natl. Acad. Sci. U. S. A.
115
,
7682
(
2018
).
9.
H.
Lan
and
D. B.
Khismatullin
, “
A numerical study of the lateral migration and deformation of drops and leukocytes in a rectangular microchannel
,”
Int. J. Multiphase Flow
47
,
73
(
2012
).
10.
A.
Lafzi
,
A. H.
Raffiee
, and
S.
Dabiri
, “
Inertial migration of a deformable capsule in an oscillatory flow in a microchannel
,”
Phys. Rev. E
102
,
063110
(
2020
).
11.
P. G.
Erlandsson
and
N. D.
Robinson
, “
Electrolysis-reducing electrodes for electrokinetic devices
,”
Electrophoresis
32
,
784
(
2011
).
12.
K.
Chaudhury
,
S.
Mandal
, and
S.
Chakraborty
, “
Droplet migration characteristics in confined oscillatory microflows
,”
Phys. Rev. E
93
,
023106
(
2016
).
13.
S.
Pawłowska
,
P.
Nakielski
,
F.
Pierini
,
I. K.
Piechocka
,
K.
Zembrzycki
, and
T. A.
Kowalewski
, “
Lateral migration of electrospun hydrogel nanofilaments in an oscillatory flow
,”
PLoS One
12
,
e0187815
(
2017
).
14.
S.
Pawłowska
,
Institute of Fundamental Technological Research Polish of Academy Sciences
(
IPPT-PAN
,
2018
), Vol.
1
, see https://rcin.org.pl/Content/142950/WA727_178307_Pawlowska-Diffusion.pdf.
15.
L.
Zhu
,
J.
Rabault
, and
L.
Brandt
, “
The dynamics of a capsule in a wall-bounded oscillating shear flow
,”
Phys. Fluids
27
,
071902
(
2015
).
16.
K.
Sarkar
and
W. R.
Schowalter
, “
Deformation of a two-dimensional viscoelastic drop at non-zero Reynolds number in time-periodic extensional flows
,”
J. Non-Newtonian fluid Mech.
95
,
315
(
2000
).
17.
Y.
Zhang
,
Y.
Han
,
L.
Zhang
,
Q.
Chen
,
M.
Ding
, and
T.
Shi
, “
Dynamic mode of viscoelastic capsules in steady and oscillating shear flow
,”
Phys. Fluids
32
,
103310
(
2020
).
18.
Z.
Lu
,
E. D.
Dupuis
,
V. K.
Patel
,
A. M.
Momen
, and
S.
Shahab
, “
Ultrasonic oscillatory two-phase flow in microchannels
,”
Phys. Fluids
33
,
032003
(
2021
).
19.
P.
Perdikaris
,
M.
Raissi
,
A.
Damianou
,
N.
Lawrence
, and
G. E.
Karniadakis
, “
Nonlinear information fusion algorithms for data-efficient multi-fidelity modelling
,”
Proc. R. Soc. A
473
,
20160751
(
2017
).
20.
H.
Babaee
,
P.
Perdikaris
,
C.
Chryssostomidis
, and
G.
Karniadakis
, “
Multi-fidelity modelling of mixed convection based on experimental correlations and numerical simulations
,”
J. Fluid Mech.
809
,
895
(
2016
).
21.
S. O.
Unverdi
and
G.
Tryggvason
, “
A front-tracking method for viscous, incompressible, multi-fluid flows
,”
J. Comput. Phys.
100
,
25
(
1992
).
22.
S.
Dabiri
and
P.
Bhuvankar
, “
Scaling law for bubbles rising near vertical walls
,”
Phys. Fluids
28
,
062101
(
2016
).
23.
M.
Bayareh
,
A.
Doostmohammadi
,
S.
Dabiri
, and
A.
Ardekani
, “
On the rising motion of a drop in stratified fluids
,”
Phys. Fluids
25
,
103302
(
2013
).
24.
X.
Chen
,
J.
Lu
, and
G.
Tryggvason
, “
Numerical simulation of self-propelled non-equal sized droplets
,”
Phys. Fluids
31
,
052107
(
2019
).
25.
W.
Li
,
J.
Lu
,
G.
Tryggvason
, and
Y.
Zhang
, “
Numerical study of droplet motion on discontinuous wetting gradient surface with rough strip
,”
Phys. Fluids
33
,
012111
(
2021
).
26.
H.
Xu
and
B. F.
Bai
, “
Unsynchronized motion of inner and outer membranes of compound capsules in shear flow
,”
Phys. Fluids
32
,
127115
(
2020
).
27.
Z. Y.
Luo
,
X. L.
Shang
, and
B. F.
Bai
, “
Development of an impulsive model of dissociation in direct simulation Monte Carlo
,”
Phys. Fluids
31
,
087105
(
2019
).
28.
E. J.
Campbell
and
P.
Bagchi
, “
A computational model of amoeboid cell swimming
,”
Phys. Fluids
29
,
101902
(
2017
).
29.
B.
Dincau
,
E.
Dressaire
, and
A.
Sauret
, “
Pulsatile flow in microfluidic systems
,”
Small
16
,
1904032
(
2020
).
30.
D.-Y.
Pan
,
Y.-Q.
Lin
,
L.-X.
Zhang
, and
X.-M.
Shao
, “
Motion and deformation of immiscible droplet in plane Poiseuille flow at low Reynolds number
,”
J. Hydrodyn.
28
,
702
(
2016
).
31.
M.
Razi
and
M.
Pourghasemi
, “
Direct numerical simulation of deformable droplets motion with uncertain physical properties in macro and micro channels
,”
Comput. Fluids
154
,
200
(
2017
).
32.
M. C.
Kennedy
and
A.
O'Hagan
, “
Predicting the output from a complex computer code when fast approximations are available
,”
Biometrika
87
(
1
),
1
(
2000
).
33.
L.
Le Gratiet
and
J.
Garnier
, “
Recursive co-kriging model for design of computer experiments with multiple levels of fidelity
,”
Int. J. Uncertainty Quantif.
4
,
365
(
2014
).
34.
G.
Vishwanathan
and
G.
Juarez
, “
Generation and application of sub-kilohertz oscillatory flows in microchannels
,”
Microfluid. Nanofluid.
24
,
1
(
2020
).
35.
M. E.
Fay
,
D. R.
Myers
,
A.
Kumar
,
C. T.
Turbyfield
,
R.
Byler
,
K.
Crawford
,
R. G.
Mannino
,
A.
Laohapant
,
E. A.
Tyburski
,
Y.
Sakurai
 et al, “
Cellular softening mediates leukocyte demargination and trafficking, thereby increasing clinical blood counts
,”
Proc. Natl. Acad. Sci. U. S. A.
113
,
1987
(
2016
).
36.
J.
Zhang
,
S.
Yan
,
D.
Yuan
,
G.
Alici
,
N.-T.
Nguyen
,
M. E.
Warkiani
, and
W.
Li
, “
Fundamentals and applications of inertial microfluidics: A review
,”
Lab Chip
16
,
10
(
2016
).
37.
J.
Rivero-Rodriguez
and
B.
Scheid
, “
Bubble dynamics in microchannels: Inertial and capillary migration forces—Corrigendum
,”
J. Fluid Mech.
855
,
1242
(
2018
).
38.
K. W.
Seo
,
Y. J.
Kang
, and
S. J.
Lee
, “
Lateral migration and focusing of microspheres in a microchannel flow of viscoelastic fluids
,”
Phys. Fluids
26
,
063301
(
2014
).
39.
J. M.
Martel
and
M.
Toner
, “
Inertial focusing in microfluidics
,”
Annu. Rev. Biomed. Eng.
16
,
371
(
2014
).
40.
D.
Stoecklein
and
D. D.
Carlo
, “
Nonlinear microfluidics
,”
Anal. Chem.
91
,
296
(
2019
).
41.
B.
Ho
and
L.
Leal
, “
Inertial migration of rigid spheres in two-dimensional unidirectional flows
,”
J. Fluid Mech.
65
,
365
(
1974
).
42.
P.
Bagchi
and
S.
Balachandar
, “
Shear versus vortex-induced lift force on a rigid sphere at moderate Re
,”
J. Fluid Mech.
473
,
379
(
2002
).
43.
P.-H.
Chan
and
L.
Leal
, “
The motion of a deformable drop in a second-order fluid
,”
J. Fluid Mech.
92
,
131
(
1979
).
44.
C. A.
Stan
,
A. K.
Ellerbee
,
L.
Guglielmini
,
H. A.
Stone
, and
G. M.
Whitesides
, “
The magnitude of lift forces acting on drops and bubbles in liquids flowing inside microchannels
,”
Lab Chip
13
,
365
(
2013
).
45.
N.
Aggarwal
and
K.
Sarkar
, “
Deformation and breakup of a viscoelastic drop in a Newtonian matrix under steady shear
,”
J. Fluid Mech.
584
,
1
(
2007
).
46.
P.
Hadikhani
,
S. M. H.
Hashemi
,
G.
Balestra
,
L.
Zhu
,
M. A.
Modestino
,
F.
Gallaire
, and
D.
Psaltis
, “
Inertial manipulation of bubbles in rectangular microfluidic channels
,”
Lab Chip
18
,
1035
(
2018
).
47.
S.
Mortazavi
and
G.
Tryggvason
, “
A numerical study of the motion of drops in Poiseuille flow. Part 1. Lateral migration of one drop
,”
J. Fluid Mech.
411
,
325
(
2000
).
48.
D. D.
Carlo
,
J. F.
Edd
,
K. J.
Humphry
,
H. A.
Stone
, and
M.
Toner
, “
Particle segregation and dynamics in confined flows
,”
Phys. Rev. Lett.
102
,
094503
(
2009
).
49.
V.
O'Brien
, “
Pulsatile fully developed flow in rectangular channels
,”
J. Franklin Inst.
300
,
225
(
1975
).
50.
M.
Karbaschi
,
A.
Javadi
,
D.
Bastani
, and
R.
Miller
, “
High frequency oscillatory flow in micro channels
,”
Colloids Surf., A
460
,
355
(
2014
).
51.
H.
Noguchi
, “
Dynamic modes of red blood cells in oscillatory shear flow
,”
Phys. Rev. E
81
,
061920
(
2010
).
52.
U.
Torres-Herrera
, “
Dynamic permeability of fluids in rectangular and square microchannels: Shift and coupling of viscoelastic bidimensional resonances
,”
Phys. Fluids
33
,
012016
(
2021
).
53.
M.
Zhao
and
P.
Bagchi
, “
Dynamics of microcapsules in oscillating shear flow
,”
Phys. Fluids
23
,
111901
(
2011
).
54.
T.
Inamuro
,
R.
Tomita
, and
F.
Ogino
, “
Lattice Boltzmann simulations of drop deformation and breakup in shear flows
,”
Int. J. Mod. Phys. B
17
,
21
(
2003
).
55.
D.
Alghalibi
,
M. E.
Rosti
, and
L.
Brandt
, “
Inertial migration of a deformable particle in pipe flow
,”
Phys. Rev. Fluids
4
,
104201
(
2019
).
56.
Q.
Wang
,
D.
Yuan
, and
W.
Li
, “
Analysis of hydrodynamic mechanism on particles focusing in micro-channel flows
,”
Micromachines
8
,
197
(
2017
).
57.
S. R.
Bazaz
,
A.
Mashhadian
,
A.
Ehsani
,
S. C.
Saha
,
T.
Krüger
, and
M. E.
Warkiani
, “
Computational inertial microfluidics: A review
,”
Lab Chip
20
,
1023
(
2020
).
58.
X.
Chen
,
C.
Xue
,
L.
Zhang
,
G.
Hu
,
X.
Jiang
, and
J.
Sun
, “
Inertial migration of deformable droplets in a microchannel
,”
Phys. Fluids
26
,
112003
(
2014
).