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.

## I. INTRODUCTION

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 Schowalter^{16} 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.

## II. METHODOLOGY

### A. The fluid dynamics simulations

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 method^{21–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:

where $\rho $ is the density of the fluid, *p* represents the pressure, **u** is the velocity vector, *t* is the time, $\tau =\mu (\u2207u+\u2207uT)$ 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

where Δ is the constant Eulerian grid size.

To generate the oscillatory flow, a cosine wave of pressure gradient with a constant amplitude [in the form of $P0\u2009cos(\omega 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(\omega 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 *U*_{0} (channel centerline velocity corresponding to the steady case) are used as the characteristic length and velocity, respectively. In other words, $x*=xW,\u2009u*=uU0,\u2009t*=tWU0,\u2009P*=P\mu U0W,\u2009T*=TWU0$ (where *T* is the period), and $\omega *=2\pi T*$. Three dimensionless parameters describe and affect the motion of the drop: (i) Reynolds number $Re=\rho U02W\mu $, expressing the ratio between inertial forces to viscous ones; (ii) Capillary number $Ca=\mu U0\gamma $, 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 ($\omega *$). The effect of the frequency can also be embedded in the Womersley number $Wo=(\omega W2\nu )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\xd7114\xd7114$ (high fidelity level) and a coarse grid of $128\xd776\xd776$ (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\xd7152\xd7152$ 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 ($\omega *$).

### B. Multi fidelity GP

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:

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:

In which, *v* is known as the signal strength, and $\u2113i$ 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,\u2026,xn)$ and $y1:n=(y1,\u2026,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,

We then model our likelihood as $p(y1:n|f1:n)\u223cN(f1:n,\sigma 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,

where the posterior mean function is

and the posterior covariance function is

With $k(x,x1:n)=(k(x,x1),\u2026,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}

where $f*t\u22121(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 *g _{t}* follows the description of GP:

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., $xt\u2286xt\u22121$) (ii) The Markov property

which translates into assuming that given the nearest point $ft\u22121(x)$, we can learn nothing more about $ft(x)$ from any other model output $ft\u22121(x\u2032)$, for any $x\u2260x\u2032$.

Since $f*t\u22121$ 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}

where $kt\rho ,\u2009ktf$, and $kt\delta $ are valid covariance functions, and $\theta t\rho ,\u2009\theta tf$, and $\theta t\delta $ are their hyperparameters, respectively. For our application, we have chosen the SE kernel function

where $\sigma 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 $t\u22652$, we have to perform predictions given uncertain inputs, where the uncertainty is propagated along each recursive step. Thus, the posterior distribution is given by^{19}

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

## III. RESULTS AND DISCUSSION

### A. Oscillatory flow

Dynamics of the droplet in an oscillatory flow in the microchannel has been studied, and the effects of *Re*, *Ca*, and $\omega *$ have been investigated. *Re* ranges between 10 and 100, *Ca* ranges between 0.09 and 10, and $\omega *$ values are chosen such that for a channel with a square cross section of $100\u2009\mu 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 by^{44}

where $Cap=\mu U0\gamma aW$ is the drop Capillary number, *V _{avg}* 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 $\omega *=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.0003\u2009W$, and the one between those of $L=4$ and $L=8$ is $0.0005\u2009W$. The results have also been shown to be grid independent, by comparing the equilibrium positions for the case of *Re* = 10, *Ca* = 1, and $\omega *=0.1$ with two different grids of $196\xd7114\xd7114$ and $256\xd7152\xd7152$. The difference between their focal distances from the center is $0.0009\u2009W$.

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 force^{9,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, $\omega *$ 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 *V _{avg}* 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

*V*determine the magnitude of the shear gradient force. For non-steady flows, including oscillatory and pulsating ones,

_{avg}*V*is time-dependent. The average of this

_{avg}*V*in each half of a corresponding periodic cycle decreases as the $\omega *$ (or

_{avg}*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 $\omega *$,

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

Stokes number in the oscillatory and pulsating flows is defined as $St=W\omega \nu $. O'Brien^{49} 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.^{49} Figure 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 ($\omega *=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 ($\omega *=8$). The average of *V _{avg}* 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

*V*is small. Consequently, the focal point at

_{avg}*Wo*= 6.3 does not follow the trend observed for the lower

*Wo*numbers and is pushed toward the wall.

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

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 $\omega *=8$ has a deformation of close to 0, which confirms that it has a very low *V _{avg}*. In fact, the droplet in this case travels about only $0.02\u2009W$ along the flow direction and remains almost spherical though the

*Ca*= 1 corresponds to a very deformable drop. Nevertheless, it migrates around $0.08\u2009W$ from the initial location to its focal point.

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, $\omega *=1$, and $\omega *=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 $\omega *$ 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.

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.015\u2009W$ 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}

Viscosity ratio in the range of $1\u2264\lambda <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.

### B. Pulsating flow

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 $\omega *$ 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.

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.012\u2009W$ 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.

Pressure gradient form . | Average migration velocity . |
---|---|

$P0(0.93+0.07\u2009\u2009cos(0.1t))$ | 0.00 191 |

$P0\u2009\u2009cos(0.1t)$ | 0.00 141 |

$P0\u2009\u2009cos(0.5t)$ | 0.00 012 |

Pressure gradient form . | Average migration velocity . |
---|---|

$P0(0.93+0.07\u2009\u2009cos(0.1t))$ | 0.00 191 |

$P0\u2009\u2009cos(0.1t)$ | 0.00 141 |

$P0\u2009\u2009cos(0.5t)$ | 0.00 012 |

Pressure gradient form . | Average migration velocity . |
---|---|

$P0\u2009\u2009cos(0.1t)$ | 0.000 389 |

$P0(0.17+0.83\u2009\u2009cos(0.5t))$ | 0.000 297 |

$P0(0.07+0.93\u2009\u2009cos(0.5t))$ | 0.000 213 |

$P0\u2009\u2009cos(0.5t)$ | 0.000 205 |

Pressure gradient form . | Average migration velocity . |
---|---|

$P0\u2009\u2009cos(0.1t)$ | 0.000 389 |

$P0(0.17+0.83\u2009\u2009cos(0.5t))$ | 0.000 297 |

$P0(0.07+0.93\u2009\u2009cos(0.5t))$ | 0.000 213 |

$P0\u2009\u2009cos(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.

Pressure gradient form . | Equilibrium distance from center . | . |
---|---|---|

$aW=0.3$ . | $aW=0.2$ . | |

$P0\u2009cos(0.1t)$ | 0.138 | 0.193 |

$P0(0.17+0.83\u2009\u2009cos(0.5t))$ | 0.166 | 0.241 |

$P0\u2009cos(0.5t)$ | 0.183 | 0.266 |

Pressure gradient form . | Equilibrium distance from center . | . |
---|---|---|

$aW=0.3$ . | $aW=0.2$ . | |

$P0\u2009cos(0.1t)$ | 0.138 | 0.193 |

$P0(0.17+0.83\u2009\u2009cos(0.5t))$ | 0.166 | 0.241 |

$P0\u2009cos(0.5t)$ | 0.183 | 0.266 |

### C. MFGP

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.

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 *R*^{2} 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.

## IV. CONCLUSIONS

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.

## ACKNOWLEDGMENTS

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

## DATA AVAILABILITY

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