In this study, the nonlinear effect of contactless bubble–bubble interactions in inertial micropumps is characterized via reduced parameter one-dimensional and three-dimensional computational fluid dynamics (3D CFD) modeling. A one-dimensional pump model is developed to account for contactless bubble-bubble interactions, and the accuracy of the developed one-dimensional model is assessed via the commercial volume of fluid CFD software, FLOW-3D. The FLOW-3D CFD model is validated against experimental bubble dynamics images as well as experimental pump data. Precollapse and postcollapse bubble and flow dynamics for two resistors in a channel have been successfully explained by the modified one-dimensional model. The net pumping effect design space is characterized as a function of resistor placement and firing time delay. The one-dimensional model accurately predicts cumulative flow for simultaneous resistor firing with inner-channel resistor placements (0.2L < x < 0.8L where L is the channel length) as well as delayed resistor firing with inner-channel resistor placements when the time delay is greater than the time required for the vapor bubble to fill the channel cross section. In general, one-dimensional model accuracy suffers at near-reservoir resistor placements and short time delays which we propose is a result of 3D bubble-reservoir interactions and transverse bubble growth interactions, respectively, that are not captured by the one-dimensional model. We find that the one-dimensional model accuracy improves for smaller channel heights. We envision the developed one-dimensional model as a first-order rapid design tool for inertial pump-based microfluidic systems operating in the contactless bubble–bubble interaction nonlinear regime.

## I. INTRODUCTION

In the last few decades, lab-on-a-chip devices have utilized a variety of passive and active flow control systems to move fluid through microchannels such as capillary action,^{1} magneto-hydrodynamics,^{2} electrophoresis,^{3} electroosmosis,^{4} mechanical peristaltic motion,^{5} centrifugal forces,^{6} hydrostatic forces,^{7} and most commonly external syringe pumps.^{8} While the microfluidics devices themselves are small, the size of the flow control system is generally significantly larger due to external fluidic pumps or power supplies giving rise to the current adage of “lab-on-a-chip” or “chip-in-a-lab.” Successful commercialization of microfluidic devices requires systems integration, scalability, and standardization of all microfluidic components including the flow control system.^{9} Thus, there is a need for an integrated, scalable, and standardized microfluidic pump source.

Inertial micropumps are an emerging micropump technology able to be integrated directly into microfluidic channels, have no moving parts, are scalable, and leverage standard semiconductor mass fabrication techniques.^{10} Unlike other pump sources, inertial micropumps can be thought of as a more general electromechanical actuator. Inertial micropumps can be used as a microfluidic pump source but also can be used for other microfluidic applications such as cell lysis^{11} and fluid mixing.^{12} This application versatility and the ability to directly integrate inertial micropumps into microchannels using existing semiconductor mass fabrication techniques hold great promise for microfluidic commercialization. Inertial micropumps were first theorized and demonstrated by Prosperetti *et al.*^{13,14} in 2000 and commercialized by Hewlett–Packard^{10} in 2010. These pumps consist of a thermal inkjet (TIJ) microresistor that locally boils liquid in a microchannel generating a vapor bubble that performs mechanical work. A voltage pulse of around 1–10 *μ*s is applied generating a heat flux in excess of 500 W/mm^{2} which causes rapid bubble expansion.^{15} Bubble collapse occurs approximately 10 *μ*s later (in a 22 × 17 *μ*m^{2} channel) as the vapor bubble loses energy due to volume expansion and thermal losses.^{16} When placed asymmetrically in a microchannel with reservoirs at either end, momentum imbalances upon bubble collapse result in a net fluid pumping effect.^{10,13} In addition to simply pumping fluid, inertial pumps were successfully used as micromixers,^{12} fluid jets,^{17} and fluid sorters/routers^{18} demonstrating the wide application of this technology.

Prosperetti *et al.* first developed a reduced parameter one-dimensional model of the pumping effect of a single resistor in a channel by considering a momentum balance at the liquid/vapor interface accounting for surface tension, pressure, and viscous forces.^{13} Later, Kornilovitch *et al.* described the pumping dynamics of a single resistor in a channel in greater detail using a similar reduced parameter one-dimensional model with an emphasis on precollapse and postcollapse pumping dynamics.^{19} Current theoretical and experimental works in inertial pump technology primarily focus on linear operating regimes, that is, when both bubble dynamics and fluid flow are fully extinguished before subsequent resistor firings. Little work to date has studied the nonlinear regime when multiple bubble events overlap in time but do not physically contact (henceforth called contactless bubble–bubble interactions). It is to be noted that this study deals exclusively with bubbles that interact *at a distance* in which one bubble modifies the liquid environment affecting another bubble but the bubbles do not occupy the same physical space. In a realistic lab-on-a-chip device, a microfluidic circuit may consist of thousands of inertial pumps in which contactless bubble–bubble interactions are unavoidable or even desirable. Thus, it is important to understand the nonlinear contactless bubble–bubble interaction regime of inertial pumps. Yuan and Prosperetti first proposed a theoretical one-dimensional model accounting for contactless bubble–bubble interactions in a three resistor system for three example resistor placement locations but did not validate the model's accuracy nor map the potential system design space^{20} and Zou *et al.* experimentally characterized laser-induced contactless bubble–bubble interactions and net flow in a two bubble system but in much larger circular channels on the order of d = 5 mm.^{21} In this study, we extensively model bubble and flow physics during the contactless bubble–bubble interactions in realistic two inertial pump systems via both one-dimensional and three-dimensional computational fluid dynamics (3D CFD) modeling. We validate our one-dimensional model using commercial volume of fluid (VOF) CFD software, FLOW-3D, and showcase full 3D bubble and flow dynamics not captured by previous and current one-dimensional models. Finally, we discuss the developed one-dimensional model's accuracy and limitations with emphasis on its predictive capabilities for future application as a rapid design tool for inertial pump-based microfluidic systems operating in the contactless bubble–bubble interaction nonlinear regime. To our knowledge, this is the first work to date that characterizes the inertial pump design space for contactless bubble–bubble interactions as well as validates the accuracy and predictive capability of a one-dimensional model for contactless bubble–bubble interactions via 3D CFD.

Inertial micropumps hold the potential to be to microfluidics what the transistor is to modern electronics. Cascades of thousands of inertial micropumps in a microfluidic circuit may one day be commonplace for commercial microfluidic devices. Interaction between bubbles will likely be unavoidable or even desirable in such a system. As such, this work provides the framework to understand contactless bubble–bubble interactions as well as formulate a one-dimensional model to quickly and accurately describe system performance without computationally expensive full 3D CFD modeling. We envision the developed one-dimensional model in this study as a future tool to aid microfluidic designers using inertial pumps in the contactless bubble–bubble interaction nonlinear regime.

## II. INERTIAL PUMP PHYSICS AND THE 1D PUMP MODEL

Here we present the 1D pump model^{19} described by Kornilovitch *et al.* which we build upon in this paper to account for contactless bubble–bubble interactions. Primary model assumptions are as follows:

Incompressible fluid.

Infinite speed of sound giving rise to instantaneous pressure wave propagation.

Surface tension forces negligible compared to viscous and pressure forces.

Bubble nucleation physics reduced to an instantaneous pressure impulse,

*p*._{v}3D bubble and flow physics reduced to tracking 1D liquid/vapor interfaces.

Bubble instantaneously occupies the channel cross section when formed.

Consider the geometry of Fig. 1(a). A microchannel of cross-sectional area *A* is connected to two fluid reservoirs with initial pressure $p1b,2b$. A thin film resistor with dimensions *w* × $\u2113$ is placed at *x _{o}*. The entire channel length is

*L*. A voltage pulse lasting 1–10

*μ*s creates a heat flux of around 500 W/mm

^{2}causing rapid vaporization.

^{15}Local boiling at the resistor's surface generates a high-pressure vapor bubble of pressure

*p*, Fig. 1(b), that pushes fluid out of the channel.

_{v}*p*can be modeled as

_{v}where *q _{o}* is the mechanical momentum imparted to both fluid columns denoted as the bubble strength (on the order of $10\u221210$ kg m/s), A is the cross section channel area, and

*p*is the vapor saturation pressure at room temperature. Here, the delta function accounts for the initial momentum kick by the high-pressure vapor bubble at t = 0 [Fig. 1(b)] after which volumetric expansion and heat transfer cause the bubble pressure to quickly drop below atmospheric pressure (

_{vr}*p*) and reach the saturation vapor pressure at ambient temperature ($pvr\u22480.3po$).

_{o}^{19}Additionally, we note that the bubble strength

*q*is a fitting parameter to reflect the initial momentum kick by the high-pressure vapor bubble. For a channel cross section of 22 × 17

_{o}*μ*m

^{2}, the bubble expansion and collapse process take around 10

*μ*s after which flow is brought to rest by viscous forces in roughly 50

*μ*s.

^{10}Bubble expansion is driven by inertia until the bubble reaches its maximum expansion volume. The short arm bubble then reverses direction while the long arm bubble continues expanding to the right, Fig. 1(d). The short arm bubble accelerates faster than the long arm bubble creating a momentum imbalance upon collapse that drives fluid flow, Fig. 1(e). The point of bubble collapse is offset from its point of expansion resulting in the primary pumping effect. The net momentum imparted to the fluid upon collapse results in the secondary pumping effect which is eventually brought to rest by viscous dissipation.

^{10}

Kornilovitch *et al.* analyzed the momentum gained/lost during the bubble expansion and collapse cycle to derive the following 1D model which was validated against experimental data for *nonoverlapping* firing pulses^{19}

where A is the cross-sectional channel area, *ρ* is the fluid density, *κ* is the characteristic strength of the viscous force^{19} derived in the Appendix (approximately 0.0184 Pa s for a 17 × 22 *μ*m^{2} rectangular channel cross section), L is the channel length, *p*_{1} and *p*_{2} are the pressures at the channel inlet/outlet, and *p _{v}* is the bubble pressure. Pressures immediately outside the channel $p1,2$ are different from bulk reservoir pressures $p1b,2b$ due to source flow from the reservoir. Prosperetti

*et al.*and Kornilovitch

*et al.*accounted for the pressure drop due to source flow by a Bernoulli term added to the bulk reservoir pressure which we adopt in this paper where m = {0, 1} is a discrete parameter to select between pressure models and H(x) is the Heaviside function

^{13,19}

In a simple 1 bubble and 1 channel network, initial positions start at the resistor center and initial velocities are found through momentum balances where *x*_{1} and *x*_{2} define the liquid/vapor interface locations as depicted in Fig. 1

However, in more complex networks, finding initial conditions may be nontrivial. Here, we describe a general procedure to identity consistent initial conditions for such nontrivial systems. At t = 0, velocities $x1\u0307$ and $x2\u0307$ are discontinuous. That is

where $t=0\u2212$ refers to the time immediately before bubble actuation and $t=0+$ refers to the time immediately after bubble actuation. In this case, $v(1,0\u2212)=0$ and $v(2,0\u2212)=0$ since the fluid starts at rest. Integrating Eqs. (2) and (3) infinitesimally about t = 0 results in only the contribution from delta-like components. Namely, from integration of *p _{v}*, $x1\u0308$, and $x2\u0308$,

Thus, we arrive at a statement of conservation of momentum which yields the same initial conditions described by Kornilovitch *et al.*

At collapse t = *t _{c}*, the bubble is at position x =

*x*with short and long arm velocities $v(1,2)c$. Here, we use the notation $v(1,2)c$ to refer to the velocity of interfaces 1 and 2, respectively, at the time of the collapse. The end precollapse positions and velocities define the starting initial conditions of the postcollapse phase. After the bubble collapse, the dynamic equations become

_{c}where a single interface *x* now models the fluid motion. From conservation of momentum, Kornilovitch *et al.* showed the initial conditions for the postcollapse bubble regime as follows:

Again, we present a generalized approach to finding the postcollapse initial conditions. Similar to the initial momentum impulse *q _{o}* during bubble expansion, the collapse phase can be represented as a momentum impulse

*q*provided by momentum imbalance at collapse. As such, postcollapse is driven by a pressure spike

_{c}*p*at t =

_{c}*t*, which is the secondary pumping effect

_{c}Thus, the dynamic equation [Eq. (13)] is modified to account for the added collapse pressure impulse

Integrating infinitesimally about t = *t _{c}* and applying conservation of momentum gives the following initial conditions:

where $\rho Ltc\u2212vtc\u2212=0$ since the combined fluid column did not exist before the collapse and can be thought of as having fluid length $Ltc\u2212$ = 0. Quantities *q _{c}* and $vtc+$ are unknown variables that are found through solving the system of equations generated by infinitesimal integration of the dynamic equations [Eq. (17)] and applying conservation of mass to the entire system [Eq. (18)]. In this simple 1 bubble and 1 channel network, $vtc+$ matches with initial postcollapse velocity given by Kornilovitch

*et al.*in Eq. (14) and

*q*reflects the momentum imparted to the fluid column upon collapse.

_{c}## III. FLOW-3D COMPUTATIONAL MODEL

### A. Model description

Thermal bubble nucleation and collapse can be represented mathematically as a complex system of coupled multiphysics equations which cannot be analytically solved. High surface heat flux, fluid-solid heat transfer coupling, rapid vapor interface expansion, surface tension effects, and small-time timescales (on the order of *μ*s's), add significant complexity to modeling the thermal bubble nucleation and expansion.

The governing equations for mass, momentum, and energy transport for an incompressible fluid with constant material properties are described by the Navier–Stokes equations

where **V** is the velocity field, *ρ* is the fluid density, P is the scalar pressure field, *μ* is the fluid viscosity, **F** is an external body force, *c _{p}* is the heat capacity at constant pressure, k is the thermal conductivity, T is the temperature, and $q\u2033\u2032gen$ is the volumetric heat generation. In this study, FLOW-3D has been used to discretize and numerically solve the Navier–Stokes equations via the VOF method. Liquid/vapor interfaces are determined following the approach developed by Hirt and Nichols.

^{22}Fractional Area/Volume Obstacle Representation (FAVOR) allows complex mesh generation taking into account both solids and liquids within a domain. Homogeneous bubble nucleation is modeled by approximating the vapor pressure inside the bubble as a function of bubble temperature through the Clausius–Clapeyron equation.

^{23}Fluid heats to the superheat temperature upon which explosive bubble nucleation occurs.

^{24}The mass flux due to evaporation/condensation is modeled using the kinetic theory described by Theofanous

*et al.*

^{25}Thus, the Clausius–Clapeyron equation maps from bubble temperature to pressure, and the kinetic theory mass flux maps mass transport from liquid to vapor phases. During bubble nucleation, energy is lost due to (a) overcoming surface tension forces during nucleation, (b) bubble volumetric expansion, and (c) heat transfer between the vapor bubble and surrounding fluid. This yields a rapid pressure drop bringing

*P*<

_{vap}*P*causing bubble collapse.

_{o}^{24}

In this model, we simplify full transient heat transfer between the resistor and surrounding fluid by assuming a prescribed time-dependent resistor temperature, $TR(t)$. This eliminates the need for a detailed material description of resistor film stacks which is proprietary in thermal microbubble systems and not available in the experimental literature data used for model validation. Using this approach, the film stack can be reduced to an effective single film of thickness, *t _{r}*. We approximate the time-dependent resistor temperature as the following rectangle function:

Heating occurs during a *t _{on}* = 1.5

*μ*s rectangular firing pulse with $TR,max$ = 873.15 K. $TR(t)$ is a reasonable approximation as, in reality, Joule heating yields a sharp exponential temperature rise which is reflected by the much faster temperature rise than other timescales of the $TR(t)$ rectangle function. Thus, full transient heat transfer can be greatly simplified. We used 300 nm grid cells at the resistor surface to properly resolve bubble nucleation dynamics. We demonstrate in Sec. III B that these model approximations yielded simulated net flows that are in excellent agreement with experimental data.

### B. Model validation

Figure 2 compares experimental bubble dynamics to simulation results in FLOW-3D at discrete time steps. FLOW-3D bubble dynamics were in good agreement with experimental images. To further validate simulation results, the net cumulative flow (meaning volume displaced) was extracted from FLOW-3D models and was in excellent agreement with experimental data in Fig. 3. To assess the mesh dependency of the FLOW-3D models, we performed a mesh analysis study for 0.50, 0.75, 1, and 1.5 *μ*m grid cell resolutions. Figure 4 shows mesh convergence at 0.50 and 0.75 *μ*m grid cell resolutions. 1 *μ*m grid cells were utilized in this study to best balance accuracy and computational time.

## IV. SINGLE CHANNEL CONTACTLESS BUBBLE–BUBBLE INTERACTION 1D MODEL DEVELOPMENT

Here, we further develop the one-dimensional pump model to account for contactless bubble–bubble interaction dynamics. We decompose the model into two main stages as shown in Fig. 5: (1) precollapse which accounts for bubble expansion and collapse and (2) postcollapse which accounts for fluid motion due to bubble collapse. During precollapse, the boundary conditions for simultaneous resistor firing (*τ* = 0) are slightly different than in the case of delayed resistor firing (*τ* > 0). During postcollapse, we develop the most general model by calculating which bubble collapses first and applying the required dynamic equations

Precollapse. *τ* = 0. Simultaneous firing.

Consider Figs. 5(a)–5(c). The two resistors are assumed to have the same dimensions and firing parameters such that the resulting bubble strengths are the same, $q1=q2=qo$. In the general case, resistor dimensions and firing parameters can vary resulting in distinct bubble strengths for each resistor. Resistor 1 fires at t = 0 and resistor 2 fires at t = *τ* resulting in bubble pressures $p1v$ and $p2v$,

We first analyze the case when *τ* = 0 and the resistors fire simultaneously. The second vapor bubble introduces two additional interfaces to the dynamic equations, *x*_{3} and *x*_{4}. Applying a momentum balance gives the following dynamic equations for precollapse:

Initial velocities are then found through infinitesimal integration about t = 0,

Initial positions are $x1(0)$ = $x2(0)$ = $xo,1$ and $x3(0)$ = $x4(0)$ = $xo,2$. Notice that firing at *τ* = 0 with equal strength resistors in a fluid at rest results in no internal flow during precollapse as $v2,0+$ = $v3,0+$ = 0 regardless of resistor placement in the channel. It should also be noted that the one-dimensional model assumption of the infinite speed of sound gives rise to instantaneous pressure wave propagation in the model. In a microfluidic channel of L = 500 *μ*m and speed of sound in the water of c = 1480 m/s, the pressure wave propagation delay would be approximately 300 ns. Thus, we define delayed resistor firing as a time delay greater than 300 ns.

Precollapse. *τ* ≥ 0. Delayed firing.

Now, consider the precollapse case where the firing of resistor 2 is delayed by some *τ* where *τ* > 0. Unlike previously, fluid is moving at all points in the channel when resistor 2 fires. Precollapse is then modeled by (a) bubble 1 expansion and (b) bubble 2 expansion. Dynamic equations and initial conditions for stage 1 are the same as the single bubble case in Eqs. (2), (3), (11), and (12). During bubble 2 expansion, the dynamic equations now contain all 4 liquid/vapor interfaces as denoted in Eqs. (25)–(28). Infinitesimal integration about t = *τ* gives the initial velocities for the system taking into account nonstationary flow due to bubble 1 expansion

We denote the velocities at the end of stage 1 as $v(1\u22124),\tau \u2212$ where $v(3\u22124),\tau \u2212=v2,\tau \u2212$. Initial positions are $x(1,2),\tau +$ = $x(1,2),\tau \u2212$ and $x3,4$ = $xo,2$ which is the resistor center point.

Postcollapse. Condition I. Bubble 1 Collapses First.

The first bubble to collapse depends on resistor placements and firing time delays. As such, there exist two conditions when modeling postcollapse: (I) bubble 1 collapses first or (II) bubble 2 collapses first. Consider condition I depicted in Figs. 5(d) and 5(e). Bubble 1 collapses at t = $tc,1$ and x = $xc,1$ with a corresponding pressure impulse $pc,1$ upon the collapse

The dynamic equations become

with initial velocities described by the following system of equations where $v(1,4),tc,1\u2212$ are precollapse velocities and $x(1,4),tc,1+$ and $v(1,4),tc,1+$ are postcollapse conditions. Initial velocities are found by infinitesimal integration about t = $tc,1$ and applying conservation of momentum

Initial positions are $x(3,4),tc,1+$ = $x(3,4),tc,1$. Consider Fig. 5(e). Bubble 2 collapses at t = $tc,2$ and x = $xc,2$ which results in a corresponding pressure impulse $pc,2$,

The dynamic equation becomes

with the following initial conditions:

Postcollapse. Condition II. Bubble 2 collapses first.

Now, we repeat the process for condition II. Figures 5(f) and 5(g) show bubble 2 collapsing first. Bubble 2 collapses at t = $tc,2$ and x = $xc,2$ with a corresponding pressure impulse $pc,2$ upon collapse. The dynamic equations become

with initial velocities described by the following system of equations where $v(1,4),tc,2\u2212$ are precollapse velocities and $x(1,4),tc,2+$ and $v(1,4),tc,2+$ are postcollapse conditions:

Initial positions are $x(1,2),tc,2+$ = $x3,4(tc,2)$. Consider Fig. 5(g). Bubble 1 collapses at t = $tc,1$ and x = $xc,1$ which results in a corresponding pressure impulse $pc,1$. The dynamic equation becomes

with the following initial conditions:

## V. RESULTS AND DISCUSSION

In this section, the developed one-dimensional model is applied to predict bubble and flow dynamics during contactless bubble–bubble interaction. Three case studies varying resistor placement and firing time delay are presented in Figs. 6–8 comparing both one-dimensional model predictions and 3D CFD results. 3D bubble and flow structures not captured by the one-dimensional model such as postcollapse vortices, bubble-reservoir interactions, and transverse bubble growth are shown and discussed as fundamental limitations to the one-dimensional model's accuracy. We conclude this section with a discussion on the system design space in terms of resistor placement and firing time delay as well as a discussion on proposed one-dimensional model accuracy regimes and constraints.

### A. Case studies

In this subsection, we apply the developed one-dimensional and FLOW-3D CFD model to describe three distinct contactless bubble–bubble interaction regimes: (1) simultaneous firing with asymmetric resistor placement, (2) simultaneous firing with symmetric resistor placement, and (3) delayed firing with symmetric resistor placement. The resistor placement is described by the normalized distance *ξ _{o}*, where

*ξ*= $xo/L$. In these case studies, the channel cross-sectional area matched our validation case (22 × 17

_{o}*μ*m

^{2}) and channel lengths are all L = 500

*μ*m. In order to accurately predict net cumulative flow, the bubble strength

*q*must first be characterized. The bubble strength is a function of fluid heat of vaporization, resistor dimensions, and resistor firing temperature/voltage, which are not changed by varying resistor placement and/or firing time delay. As such, the bubble strength is a constant so long as the fluid and resistor dimensions remain the same, as in this study. The bubble strength was extracted from FLOW-3D CFD data of a single resistor placed at

_{o}*ξ*= 0.20 in a L = 500

_{o}*μ*m channel following the curve fitting approach of Kornilovitch

*et al.*

^{16}The bubble strength was found to be

*q*= 4.26 × 10

_{o}^{−10 }kg m/s and used throughout this study.

In Fig. 6, two resistors are placed asymmetrically ($\xi o,1$ = 0.20 and $\xi o,2$ = 0.70) in a channel of length L = 500 *μ*m. Both resistors are fired simultaneously at t = 0 *μ*s. As shown in Fig. 6(a), simultaneous firing causes little movement of bubble interfaces *x*_{1} and *x*_{2} consistent with one-dimensional model physics. Figure 6(b) shows the predicted cumulative flow. The one-dimensional model is in excellent agreement with FLOW-3D CFD. The one-dimensional model predicts no flow during the precollapse phase as both resistors have the same bubble strength so internal flow integrated across the flux plane shown is 0. However, FLOW-3D CFD shows cumulative flow during precollapse. Note in Figs. 6(c) and 6(d) that the bubble requires greater than 4 *μ*s to fully fill the channel cross section. As such, during realistic bubble expansion, fluid can flow over the bubble unlike the assumption in the one-dimensional model where the bubbles occupy the entire channel cross section instantaneously. Nevertheless, the predicted cumulative flow is in excellent agreement with 3D CFD. Figures 6(e) and 6(f) show postcollapse fluid vortices that occur when the vapor bubble fully collapses. This is the first time such behavior has been modeled and we hypothesize that these postcollapse 3D fluid vortexes are the mechanism for bubble-based mixing observed in past experimental studies.^{12} Additionally, we observe that it takes approximately 6 *μ*s for the bubble to fill the 22 × 17 *μ*m^{2} channel cross section during its transverse growth phase. The bubble rebound effect upon collapse shown in Fig. 6(f) agrees with previous experimental and numerical studies.^{16}

In Fig. 7, two resistors are placed symmetrically ($\xi o,1$ = 0.20 and $\xi o,2$ = 0.80) in a channel of length L = 500 *μ*m. Both resistors are fired simultaneously at t = 0 *μ*s. Similar to before, simultaneous firing causes little movement of bubble interfaces *x*_{1} and *x*_{2} as shown in Fig. 7(a). Theoretically, the predicted cumulative flow is 0 due to symmetric resistor placement, as the one-dimensional model shows. However, FLOW-3D CFD results show negligible flow during precollapse (t < 18 *μ*s) and then a small spike at the point of bubble collapse (t = 18 *μ*s) giving rise to a nearly negligible cumulative flow of 0.21 pl. We conjecture that the cumulative flow spike arises from slight numerical errors which cause the ideal resistor symmetry condition to be broken. We note that although symmetric resistor placement does not generate zero net flow as theoretically predicted, the net flow is significantly smaller than typical displacements in the nonsymmetric conditions. In general, we suggest that 3D postcollapse vortices, transverse bubble growth, and bubble-reservoir interactions represent fundamental limitations to one-dimensional model accuracy. Such differences could be accounted for in the one-dimensional model by using experimental correlations of the bubble strength *q _{o}* with resistor placements, channel lengths, and channel cross sections (similar to experimental correlations commonly done in heat transfer studies) but such an analysis is beyond the scope of this paper and will be the focus of future work.

In Fig. 8, two resistors are placed symmetrically ($\xi o,1$ = 0.20 and $\xi o,2$ = 0.80) in a channel of length L = 500 *μ*m. Now, the resistors are fired asynchronously where resistor 1 fires at t = 0 *μ*s and resistor 2 fires at t = 10 *μ*s. The *τ* = 10 *μ*s delay results in net cumulative flow whereas in the simultaneous firing case no flow was expected. Later, we show that the time delay can be used to obtain greater flow control. Consider Figs. 8(a) and 8(b), bubble dynamics and cumulative flow were accurately predicted by the one-dimensional model until bubble 1 collapses. After bubble 1 collapses, the flow rebounds from reverse flow to forward flow but at a lesser extent than predicted by CFD results. As shown in Figs. 8(c) and 8(d), we conjecture that it is the 3D transverse bubble growth and postcollapse vortices that cause the system to deviate from theoretical predictions resulting in a faster bubble collapse time predicted by CFD than the one-dimensional model. Specifically, the CFD model accounts for 3D transverse bubble growth which dissipates the bubble's mechanical momentum resulting in a faster bubble collapse than the one-dimensional model. Since the one-dimensional model assumes all mechanical momentum imparted by the bubble nucleation goes into the fluid flow, the vapor bubble takes longer to collapse. Figure 8(c) shows flow over the second bubble during transverse growth. Figure 8(d) illustrates 3D fluid vortices which dissipate energy that the one-dimensional model assumes goes into the forward flow. Figure 8(e) highlights the secondary pumping effect from CFD simulations.

### B. System design space and 1D model accuracy

Here, we map the design space of contactless bubble–bubble interaction for two resistors in a channel using the developed one-dimensional model and discuss the one-dimensional model's accuracy and limitations. Figures 9(a)–9(d) show the net cumulative flow as a function of normalized resistors 1 and 2 placement, $\xi o,1$ and $\xi o,2$, respectively, for two-time delays and channel lengths with the condition of zero cumulative flow marked. For each time delay, decreasing the channel length caused the magnitude of cumulative flow to increase. Additionally, the band of max cumulative flow resistor placement positions became more linearly distributed in smaller channel lengths than in the longer channel. For each channel length, delayed resistor firing resulted in increased nonlinearity of the system as can be seen by the condition of zero cumulative flow. In the simultaneous firing cases, the condition for zero cumulative flow is symmetric resistor placement or $\xi o,2$ = 1 − $\xi o,1$. In addition to assessing the resistor placement design space, we assess the effect of firing time delay on cumulative flow in Fig. 10. Here, resistor 1 is placed at $\xi o,1$ = 0.20 and the placement of resistor 2 ($\xi o,2$) and the firing time delay (*τ*) is allowed to vary. Figure 10(c) shows that the firing time delay can, in specific cases, cause the net flow to shift from forward to reverse flow and thus can be used as a flow controller to achieve a desired cumulative flow rate.

Next, we discuss the accuracy and limitations of the one-dimensional model. An accuracy analysis was performed using FLOW-3D CFD results as the predictive standard upon which to compare one-dimensional model results. Figures 9(e), 9(f), 10(b), 10(c), 11(b), and 11(c) depict both FLOW-3D CFD and one-dimensional model predicted cumulative flows for sample points spanning different contactless bubble-bubble interaction regimes. Consider Fig. 9(e). The one-dimensional model was the most accurate in predicting cumulative flow for simultaneous resistor firing (*τ* = 0) with inner-channel resistor placements (0.2L < x < 0.8L where L is the channel length). As the resistor placement moved closer to the reservoirs, bubble-reservoir 3D flow interactions became significant causing the one-dimensional model's accuracy to diverge. These trends are also observed in Fig. 9(f) when *τ* = 5 *μ*s for delayed firing. However, with delayed resistor firing in a H = 17 *μ*m channel, the bubble takes approximately 6 *μ*s to fully fill the 22 × 17 *μ*m^{2} channel cross section. Thus, fluid moves around the second bubble during its expansion accounting for observed net flow. This transverse growth effect is not captured by the one-dimensional model causing accuracy to degrade in the delayed resistor firing regime. Consider Figs. 10(b) and 10(c). By taking sample points with increasing firing delays, the effect of transverse bubble growth is more pronounced. For small firing delays, the one-dimensional model fails to accurately predict the cumulative flow. As the firing delay increases, the model accuracy improves. The one-dimensional model reaches improved accuracy beyond *τ* = 6 *μ*s which is when the bubble can be safely assumed to fully take up the cross section of the channel and thus transverse bubble growth effects are minimized. Figure 11 further emphasizes the impact of transverse bubble growth effects on model accuracy. In this case, the channel height is halved to H = 8.5 *μ*m upon which the vapor bubble plugs the channel cross section in approximately 2 *μ*s. Here, one-dimensional model accuracy is improved due to minimizing the impact of transverse bubble growth effects.

We conclude this section by discussing the importance of developing reduced parameter one-dimensional models for inertial pumps. Design space analysis provides a way to predict resistor placements as well as firing time delays needed to achieve a required net cumulative flow for systems operating in the contactless bubble–bubble interaction regime. To use inertial pumps in any application, one must first understand and control fluid flow by characterizing the design space. In the design space analysis examples presented in this study, each contour plot consisted of 40 × 40 = 1600 data points. FLOW-3D CFD required 160 core hours per sample run to compute. As such, a single parameter space mapping of 1600 points would require 256 000 core hours or 667 days (with a 16 core computer) to compute in comparison to <1 min run time per contour plot for the one-dimensional model. By developing accurate reduced parameter one-dimensional models, first-order systems analysis can be rapidly performed to inform design decisions of inertial pump-based microfluidic systems operating in the contactless bubble–bubble interaction regime.

## VI. CONCLUSIONS

The present study deals with a fundamental understanding of contactless bubble-bubble interactions in inertial pump microfluidic systems through both one-dimensional modeling and 3D CFD. A reduced parameter one-dimensional model was developed and validated through 3D CFD. One-dimensional model accuracy was demonstrated for inner-channel resistor placements (0.2L < x < 0.8L) with simultaneous resistor firing as well as for delayed resistor firing with inner-channel resistor placements where the firing delay was greater than 6 *μ*s. Model accuracy was found to improve with a decrease in channel height. It was suggested that model limitations were due to 3D fluid structures such as bubble-reservoir interactions and postcollapse vortices as well as 3D bubble effects such as transverse bubble growth not captured by the one-dimensional model. Delayed resistor firing was found to enhance the nonlinearity of the system and, in specific cases, cause the net flow to shift from forward to reverse flow. The developed one-dimensional model provided significant time savings over CFD where a 1600 data point sweep took less than 1 min for the 1D model compared to 256 000 core hours for CFD. Furthermore, it was proposed that the observed 3D postcollapse vortices are the mechanism behind inertial pump-based micromixing which was experimentally demonstrated in previous work.

Cascades of thousands of inertial micro-pumps in a microfluidic circuit may one day be commonplace for commercial microfluidic devices in which interaction between bubbles will likely be unavoidable. As such, this work provides the framework to understand contactless bubble–bubble interactions as well as formulate a one-dimensional model to quickly and accurately describe system performance without computationally expensive full 3D CFD modeling. We envision that this one-dimensional model will be an essential tool for microfluidic designers using inertial pumps in the contactless bubble–bubble interaction nonlinear regime.

## ACKNOWLEDGMENTS

The authors would like to thank Pavel Kornilovitch for insights and discussions on mathematical modeling as well as Alex Govyadinov for fruitful discussions of this work.

This material is based upon work supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE 1650115. Any opinion, findings, and conclusions, or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.

## DATA AVAILABILITY

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

## NOMENCLATURE

*A*[m^{2}]channel cross-sectional area

*L*[m]channel length

*m*[{0, 1}]discrete model pressure model select parameter

*p*[Pa]_{o}atmospheric pressure

*p*[Pa]_{c}postcollapse pressure spike

*p*[Pa]_{v}bubble pressure

*p*[Pa]_{vr}vapor saturation pressure at room temperature

- $p1b,2b$ [Pa]
bulk reservoir pressure

- $p1,2$ [Pa]
pressure at channel entrance region

*q*[kg m/s]_{c}postcollapse momentum impulse

*q*[kg m/s]_{o}bubble strength

*t*[s]_{c}time at bubble collapse

*x*[m]_{o}resistor placement

*x*[m]_{c}location of bubble collapse

- $x1,2,3,4$ [m]
location of liquid/vapor interface

*κ*[Pa s]viscous dissipation coefficient

*v*[m/s]_{c}velocity at bubble collapse

*ρ*[kg/m^{3}]fluid density

*τ*[s]resistor firing offset

### APPENDIX: CALCULATION OF κ FOR A RECTANGULAR CHANNEL OF CROSS-SECTIONAL AREA a × b

Here, we follow the derivation for the viscous stress dissipation factor *κ* put forth by Kornilovitch *et al.*^{19} and apply it to a rectangular channel cross section. Let a and b be the rectangle dimensions along the y axis and z axis, respectively. The series solution for the velocity profile and flow rate is given by the following expressions:^{26}

where

The velocity profile can be re-written in terms of the average velocity, $\u27e8v\u27e9$ = Q/A

Calculation of the viscous stress tensor and integrating over the wall perimeter gives the total viscous force which can be rewritten as F = $\kappa \u27e8v\u27e9L$ to agree with the one-dimensional model. In the case of a rectangular channel cross section^{26}

where