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.

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 lysis11 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–Packard10 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/mm2 which causes rapid bubble expansion.15 Bubble collapse occurs approximately 10 μs later (in a 22 × 17 μm2 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/routers18 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 space20 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.

Here we present the 1D pump model19 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:

  1. Incompressible fluid.

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

  3. Surface tension forces negligible compared to viscous and pressure forces.

  4. Bubble nucleation physics reduced to an instantaneous pressure impulse, pv.

  5. 3D bubble and flow physics reduced to tracking 1D liquid/vapor interfaces.

  6. 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 × is placed at xo. The entire channel length is L. A voltage pulse lasting 1–10 μs creates a heat flux of around 500 W/mm2 causing rapid vaporization.15 Local boiling at the resistor's surface generates a high-pressure vapor bubble of pressure pv, Fig. 1(b), that pushes fluid out of the channel. pv can be modeled as

pv(t)=qoAδ(t)+pvr,
(1)

where qo is the mechanical momentum imparted to both fluid columns denoted as the bubble strength (on the order of 1010 kg m/s), A is the cross section channel area, and pvr 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 (po) and reach the saturation vapor pressure at ambient temperature (pvr0.3po).19 Additionally, we note that the bubble strength qo is a fitting parameter to reflect the initial momentum kick by the high-pressure vapor bubble. For a channel cross section of 22 × 17 μm2, 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 

FIG. 1.

Schematic of the inertial pump 1D model. (a) Depicts a microchannel with two fluid reservoirs of pressure p1b,2b and channel inlet/outlet pressure p1,2. At its simplest case, p1b,2b=po and the fluid starts at rest. (b)–(e) describes the pumping process. (b) Local boiling generates a high-pressure vapor bubble, pv, that pushes fluid out of the channel. (c) pv decreases due to bubble volumetric expansion and heat transfer and rapidly drops below atmospheric pressure to pvr. Bubble expansion is driven by inertia until the bubble reaches its maximum expansion volume. (d) The short arm bubble (to the left of the resistor) reverses direction while the long arm bubble (to the right of the resistor) continues expanding to the right. (e) The short arm bubble accelerates faster than the long arm bubble resulting in a momentum imbalance that drives fluid flow upon collapse. The point of collapse is offset from the location of bubble expansion giving rise to the primary pumping effect. The net momentum imparted to the fluid generates the secondary pumping effect, which is eventually brought to rest by viscous dissipation.

FIG. 1.

Schematic of the inertial pump 1D model. (a) Depicts a microchannel with two fluid reservoirs of pressure p1b,2b and channel inlet/outlet pressure p1,2. At its simplest case, p1b,2b=po and the fluid starts at rest. (b)–(e) describes the pumping process. (b) Local boiling generates a high-pressure vapor bubble, pv, that pushes fluid out of the channel. (c) pv decreases due to bubble volumetric expansion and heat transfer and rapidly drops below atmospheric pressure to pvr. Bubble expansion is driven by inertia until the bubble reaches its maximum expansion volume. (d) The short arm bubble (to the left of the resistor) reverses direction while the long arm bubble (to the right of the resistor) continues expanding to the right. (e) The short arm bubble accelerates faster than the long arm bubble resulting in a momentum imbalance that drives fluid flow upon collapse. The point of collapse is offset from the location of bubble expansion giving rise to the primary pumping effect. The net momentum imparted to the fluid generates the secondary pumping effect, which is eventually brought to rest by viscous dissipation.

Close modal

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 pulses19 

Aρx1x¨1+κx1ẋ1=(p1pv)A,
(2)
Aρ(Lx2)x¨2+κ(Lx2)ẋ2=(pvp2)A,
(3)

where A is the cross-sectional channel area, ρ is the fluid density, κ is the characteristic strength of the viscous force19 derived in the  Appendix (approximately 0.0184 Pa s for a 17 × 22 μm2 rectangular channel cross section), L is the channel length, p1 and p2 are the pressures at the channel inlet/outlet, and pv 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 function13,19

p1=p1bm2ρx1̇2H(x1̇),
(4)
p2=p2bm2ρx2̇2H(x2̇).
(5)

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 x1 and x2 define the liquid/vapor interface locations as depicted in Fig. 1 

ShortArmLongArmx1(0)=xox2(0)=xoẋ1(0)=qoρAxoẋ2(0)=qoρA(Lxo).
(6)

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̇ and x2̇ are discontinuous. That is

limt0ẋ1=v(1,0),
(7)
limt0+ẋ1=v(1,0+),
(8)
limt0ẋ2=v(2,0),
(9)
limt0+ẋ2=v(2,0+),
(10)

where t=0 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)=0 and v(2,0)=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 pv, x1̈, and x2̈,

ρxov1,0+=qoA,
(11)
ρ(Lxo)v2,0+=qoA.
(12)

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

At collapse t = tc, the bubble is at position x = xc 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

AρLx¨+κLẋ=(p1p2)A,
(13)

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:

PostCollapsex(tc)=xcẋ(tc)=ρAxcv1c+ρA(Lxc)v2cρAL.
(14)

Again, we present a generalized approach to finding the postcollapse initial conditions. Similar to the initial momentum impulse qo during bubble expansion, the collapse phase can be represented as a momentum impulse qc provided by momentum imbalance at collapse. As such, postcollapse is driven by a pressure spike pc at t = tc, which is the secondary pumping effect

pc(t)=qcAδ(ttc).
(15)

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

AρLx¨+κLẋ=(pc+p1p2)A.
(16)

Integrating infinitesimally about t = tc and applying conservation of momentum gives the following initial conditions:

ρLvtc+ρLtcvtc=qcA,
(17)
ρAxcv1c+ρA(Lxc)v2c=ρALvtc+,
(18)

where ρLtcvtc=0 since the combined fluid column did not exist before the collapse and can be thought of as having fluid length Ltc = 0. Quantities qc 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 qc reflects the momentum imparted to the fluid column upon collapse.

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

·V=0,
(19)
ρDVDt=P+μ2V+ρF,
(20)
ρcpDTDt=k2T+qgen,
(21)

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, cp is the heat capacity at constant pressure, k is the thermal conductivity, T is the temperature, and qgen 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 Pvap < Po causing bubble collapse.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, tr. We approximate the time-dependent resistor temperature as the following rectangle function:

TR(t)={0t<tonTR,maxtontτ0t>τ.
(22)

Heating occurs during a ton = 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.

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.

FIG. 2.

FLOW-3D CFD validation. FLOW-3D bubble dynamics compared to experimental data at discrete time steps. Experimental images adapted with permission from Govyadinov et al., Microfluid. Nanofluid. 20, 73 (2016). Copyright 2016 Author(s), licensed under a Creative Commons Attribution (CC BY) license.16 The U-shaped channel geometry was matched with that used for experimental bubble dynamics data16 where L = 403 μm, W = 22 μm, and H = 17 μm. No-slip boundary conditions were applied to all boundaries except the bottom of the reservoir in which a static pressure boundary condition was applied.

FIG. 2.

FLOW-3D CFD validation. FLOW-3D bubble dynamics compared to experimental data at discrete time steps. Experimental images adapted with permission from Govyadinov et al., Microfluid. Nanofluid. 20, 73 (2016). Copyright 2016 Author(s), licensed under a Creative Commons Attribution (CC BY) license.16 The U-shaped channel geometry was matched with that used for experimental bubble dynamics data16 where L = 403 μm, W = 22 μm, and H = 17 μm. No-slip boundary conditions were applied to all boundaries except the bottom of the reservoir in which a static pressure boundary condition was applied.

Close modal
FIG. 3.

Cumulative flow model validation illustrates the cumulative volume displaced per pulse of FLOW-3D model data in comparison to experimental data16 (republished with permission from Govyadinov et al., Microfluid. Nanofluid. 20, 73 (2016). Copyright 2016 Author(s), licensed under a Creative Commons Attribution (CC BY) license) and 1D model fit. A τ = 1.5 μs firing pulse was used throughout this study. L = 403 μm, w = 22 μm, h = 17 μm, μ = 8.9 × 104 Pa s, and ρ = 1000 kg/m3. 1D model fit parameters were qo = 3.52 × 1010 kg m/s and xo = 73 μm.

FIG. 3.

Cumulative flow model validation illustrates the cumulative volume displaced per pulse of FLOW-3D model data in comparison to experimental data16 (republished with permission from Govyadinov et al., Microfluid. Nanofluid. 20, 73 (2016). Copyright 2016 Author(s), licensed under a Creative Commons Attribution (CC BY) license) and 1D model fit. A τ = 1.5 μs firing pulse was used throughout this study. L = 403 μm, w = 22 μm, h = 17 μm, μ = 8.9 × 104 Pa s, and ρ = 1000 kg/m3. 1D model fit parameters were qo = 3.52 × 1010 kg m/s and xo = 73 μm.

Close modal
FIG. 4.

Mesh analysis. Mesh analysis study using 0.50, 0.75, 1, and 1.5 μm grid cell resolutions showing convergence to experimental data16 (republished with permission from Govyadinov et al., Microfluid. Nanofluid. 20, 73 (2016). Copyright 2016 Author(s), licensed under a Creative Commons Attribution (CC BY) license).

FIG. 4.

Mesh analysis. Mesh analysis study using 0.50, 0.75, 1, and 1.5 μm grid cell resolutions showing convergence to experimental data16 (republished with permission from Govyadinov et al., Microfluid. Nanofluid. 20, 73 (2016). Copyright 2016 Author(s), licensed under a Creative Commons Attribution (CC BY) license).

Close modal

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

FIG. 5.

Contactless bubble–bubble interaction model development describes the model layout consisting of (a)–(c) stage 1, precollapse and (d)–(g) stage 2, postcollapse. The precollapse stage is subdivided to account for simultaneous resistor firing (τ = 0) and delayed resistor firing (τ > 0). The postcollapse stage is subdivided to account for the most general case where the model calculates which bubble collapses first and applies the correct dynamic equations. t = tc,1 and t = tc,2 refer to the time of collapse for bubbles 1 and 2, respectively.

FIG. 5.

Contactless bubble–bubble interaction model development describes the model layout consisting of (a)–(c) stage 1, precollapse and (d)–(g) stage 2, postcollapse. The precollapse stage is subdivided to account for simultaneous resistor firing (τ = 0) and delayed resistor firing (τ > 0). The postcollapse stage is subdivided to account for the most general case where the model calculates which bubble collapses first and applies the correct dynamic equations. t = tc,1 and t = tc,2 refer to the time of collapse for bubbles 1 and 2, respectively.

Close modal

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,

p1v=qoAδ(t)+pvr,
(23)
p2v=qoAδ(tτ)+pvr.
(24)

We first analyze the case when τ = 0 and the resistors fire simultaneously. The second vapor bubble introduces two additional interfaces to the dynamic equations, x3 and x4. Applying a momentum balance gives the following dynamic equations for precollapse:

(p1p1v)Aκx1x1̇=ρAx1x1̈,
(25)
(p1vp2v)Aκ(x3x2)x2̇=ρA(x3x2)x2̈,
(26)
(p1vp2v)Aκ(x3x2)x3̇=ρA(x3x2)x3̈,
(27)
(p2vp2)Aκ(Lx4)x4̇=ρA(Lx4)x4̈.
(28)

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

qo=ρAxo,1v1,0+,
(29)
v2,0+=0,
(30)
v3,0+=0,
(31)
qo=ρA(Lxo,2)v4,0+.
(32)

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

v1,τ+=v1,τ,
(33)
qo=ρA(xo,2x2,τ)(v2,τ+v2,τ),
(34)
qo=ρA(xo,2x2,τ)(v3,τ+v3,τ),
(35)
qo=ρA(Lxo,2)(v4,τ+v4,τ).
(36)

We denote the velocities at the end of stage 1 as v(14),τ where v(34),τ=v2,τ. Initial positions are x(1,2),τ+ = x(1,2),τ 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

pc,1=qcAδ(ttc,1).
(37)

The dynamic equations become

(pc,1+p1p2v)Aκx3x3̇=ρAx3x3̈,
(38)
(p2vp2)Aκ(Lx4)x4̇=ρA(Lx4)x4̈,
(39)

with initial velocities described by the following system of equations where v(1,4),tc,1 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

qc=ρAx3(tc,1)(v3,tc,1+v3,tc,1),
(40)
v4,tc,1+=v4,tc,1,
(41)
ρAx3(tc,1)v3,tc,1++ρA(Lx4(tc,1))v4,tc,1+=ρAxc,1v1,tc,1+ρA(x3(tc,1)xc,1)v2,tc,1+ρA(Lx4(tc,1))v4,tc,1.
(42)

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,

pc,2=qcAδ(ttc,2).
(43)

The dynamic equation becomes

(pc,2+p1p2)AκLẋ=ρALx¨,
(44)

with the following initial conditions:

qc=ρALvtc,2+,
(45)
ρAxc,2v3,tc,2+ρA(Lxc,2)v4,tc,2=ρALvtc,2+,
(46)
x(tc,2)=xc,2.
(47)

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

(p1p1v)Aκx1x1̇=ρAx1x1¨,
(48)
(pc,2+p1vp2)Aκ(Lx2)x2̇=ρA(Lx2)x2¨,
(49)

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

v1,tc,2+=v1,tc,2,
(50)
qc=ρA(Lx2(tc,2))(v2,tc,2+v2,tc,2),
(51)
ρAx1(tc,2)v1,tc,2++ρA(Lx2(tc,2))v2,tc,2+=ρAx1(tc,2)v1,tc,2+ρA(xc,2x2(tc,2))v3,tc,2+ρA(Lx4(tc,2))v4,tc,2.
(52)

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

(pc,1+p1p2)AκLẋ=ρALx¨,
(53)

with the following initial conditions:

qc=ρALvtc,1+,
(54)
ρAxc,1v1,tc,1+ρA(Lxc,1)v2,tc,1=ρALvtc,1+,
(55)
x(tc,1)=xc,1.
(56)

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.

FIG. 6.

Contactless bubble–bubble interaction τ = 0, ξo,1 = 0.2, ξo,2 = 0.7, L = 500 μm. (a) and (b) describe one-dimensional model data overlayed on FLOW-3D CFD data. (a) depicts the bubble edges over time and (b) shows the cumulative flow through the microchannel over time. Resistor placement and firing delay were ξo,1 = 0.2, ξo,2 = 0.7, and τ = 0 with a one-dimensional model bubble strength of qo = 4.26 × 1010 kg m/s. (c)–(f) illustrates 3D CFD flow structures and bubble dynamics. (c) highlights both axial and transverse bubble growth. We note flow along the axial direction until the vapor bubble plugs the channel >4 μs after initial expansion (d). Postcollapse vortex structures are shown in (e) and (f). Bubble edges are labeled x14 and net flow was measured at the flux plane.

FIG. 6.

Contactless bubble–bubble interaction τ = 0, ξo,1 = 0.2, ξo,2 = 0.7, L = 500 μm. (a) and (b) describe one-dimensional model data overlayed on FLOW-3D CFD data. (a) depicts the bubble edges over time and (b) shows the cumulative flow through the microchannel over time. Resistor placement and firing delay were ξo,1 = 0.2, ξo,2 = 0.7, and τ = 0 with a one-dimensional model bubble strength of qo = 4.26 × 1010 kg m/s. (c)–(f) illustrates 3D CFD flow structures and bubble dynamics. (c) highlights both axial and transverse bubble growth. We note flow along the axial direction until the vapor bubble plugs the channel >4 μs after initial expansion (d). Postcollapse vortex structures are shown in (e) and (f). Bubble edges are labeled x14 and net flow was measured at the flux plane.

Close modal
FIG. 7.

Contactless bubble–bubble interaction τ = 0, ξo,1 = 0.2, ξo,2 = 0.8, L = 500 μm. (a) and (b) describe one-dimensional model data overlayed on FLOW-3D CFD data. (a) depicts the bubble edges over time and (b) depicts the cumulative flow through the microchannel over time. Resistor placement and firing delay were ξo,1 = 0.2, ξo,2 = 0.8, and τ = 0 with a one-dimensional model bubble strength of qo = 4.26 × 1010 kg m/s. (c) and (d) illustrates 3D CFD flow structures and bubble dynamics. (c) shows symmetric bubble expansion with little internal flow while (d) shows postcollapse vortices occurring approximately 10 μs after bubble collapse. Bubble edges are labeled x14 and net flow was measured at the flux plane.

FIG. 7.

Contactless bubble–bubble interaction τ = 0, ξo,1 = 0.2, ξo,2 = 0.8, L = 500 μm. (a) and (b) describe one-dimensional model data overlayed on FLOW-3D CFD data. (a) depicts the bubble edges over time and (b) depicts the cumulative flow through the microchannel over time. Resistor placement and firing delay were ξo,1 = 0.2, ξo,2 = 0.8, and τ = 0 with a one-dimensional model bubble strength of qo = 4.26 × 1010 kg m/s. (c) and (d) illustrates 3D CFD flow structures and bubble dynamics. (c) shows symmetric bubble expansion with little internal flow while (d) shows postcollapse vortices occurring approximately 10 μs after bubble collapse. Bubble edges are labeled x14 and net flow was measured at the flux plane.

Close modal
FIG. 8.

Contactless bubble–bubble interaction τ = 10 μs, ξo,1 = 0.2, ξo,2 = 0.8, L = 500 μm. (a) and (b) describe one-dimensional model data overlayed on FLOW-3D CFD data. (a) depicts the bubble edges over time and (b) depicts the cumulative flow through the microchannel over time. Resistor placement and firing delay were ξo,1 = 0.2, ξo,2 = 0.8, and τ = 10 μs with a one-dimensional model bubble strength of qo = 4.26 × 1010 kg m/s. (c)–(f) illustrates 3D CFD flow structures and bubble dynamics. (c) shows initial bubble expansion, (d) shows the first postcollapse vortex formation, (e) showcases the momentum imbalance between fluid legs giving rise to the secondary pumping effect, and (f) shows the second postcollapse vortex formation. Bubble edges are labeled x14 and net flow was measured at the flux plane.

FIG. 8.

Contactless bubble–bubble interaction τ = 10 μs, ξo,1 = 0.2, ξo,2 = 0.8, L = 500 μm. (a) and (b) describe one-dimensional model data overlayed on FLOW-3D CFD data. (a) depicts the bubble edges over time and (b) depicts the cumulative flow through the microchannel over time. Resistor placement and firing delay were ξo,1 = 0.2, ξo,2 = 0.8, and τ = 10 μs with a one-dimensional model bubble strength of qo = 4.26 × 1010 kg m/s. (c)–(f) illustrates 3D CFD flow structures and bubble dynamics. (c) shows initial bubble expansion, (d) shows the first postcollapse vortex formation, (e) showcases the momentum imbalance between fluid legs giving rise to the secondary pumping effect, and (f) shows the second postcollapse vortex formation. Bubble edges are labeled x14 and net flow was measured at the flux plane.

Close modal

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 ξo = xo/L. In these case studies, the channel cross-sectional area matched our validation case (22 × 17 μm2) and channel lengths are all L = 500 μm. In order to accurately predict net cumulative flow, the bubble strength qo 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 μm channel following the curve fitting approach of Kornilovitch et al.16 The bubble strength was found to be qo = 4.26 × 10−10 kg m/s and used throughout this study.

In Fig. 6, two resistors are placed asymmetrically (ξo,1 = 0.20 and ξ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 x1 and x2 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 μm2 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 (ξo,1 = 0.20 and ξ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 x1 and x2 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 qo 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 (ξo,1 = 0.20 and ξ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.

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, ξo,1 and ξ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 ξo,2 = 1 − ξ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 ξo,1 = 0.20 and the placement of resistor 2 (ξ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.

FIG. 9.

Resistor placement design space analysis using the developed one-dimensional model with H = 17 μm. (a) illustrates the cumulative flow resulting from varying xo,1 and xo,2 when τ = 0 and L = 300 μm. (b) illustrates the cumulative flow resulting from varying xo,1 and xo,2 when τ = 5 μs and L = 300 μm. (c) illustrates the cumulative flow resulting from varying xo,1 and xo,2 when τ = 0 and L = 500 μm. (d) illustrates the cumulative flow resulting from varying xo,1 and xo,2 when τ = 5 μs and L = 500 μm. Dashed black line is the condition for no net flow and separates forward from reverse flow. Black diamonds are points used in the accuracy analysis comparing predicted cumulative flow from the one-dimensional model to the FLOW-3D CFD model shown in (e) and (f).

FIG. 9.

Resistor placement design space analysis using the developed one-dimensional model with H = 17 μm. (a) illustrates the cumulative flow resulting from varying xo,1 and xo,2 when τ = 0 and L = 300 μm. (b) illustrates the cumulative flow resulting from varying xo,1 and xo,2 when τ = 5 μs and L = 300 μm. (c) illustrates the cumulative flow resulting from varying xo,1 and xo,2 when τ = 0 and L = 500 μm. (d) illustrates the cumulative flow resulting from varying xo,1 and xo,2 when τ = 5 μs and L = 500 μm. Dashed black line is the condition for no net flow and separates forward from reverse flow. Black diamonds are points used in the accuracy analysis comparing predicted cumulative flow from the one-dimensional model to the FLOW-3D CFD model shown in (e) and (f).

Close modal
FIG. 10.

Firing delay design space using the developed one-dimensional model with H = 17 μm. (a) Illustrates the cumulative flow resulting from varying τ and ξo,2 when ξo,1=0.2, L = 500 μm, and H = 17 μm. Dashed black line is the condition for no net flow and separates forward from reverse flow. Black diamonds and triangles are sample points used in the accuracy analysis comparing predicted cumulative flow from the one-dimensional model to the FLOW-3D CFD model shown in (b) and (c), respectively.

FIG. 10.

Firing delay design space using the developed one-dimensional model with H = 17 μm. (a) Illustrates the cumulative flow resulting from varying τ and ξo,2 when ξo,1=0.2, L = 500 μm, and H = 17 μm. Dashed black line is the condition for no net flow and separates forward from reverse flow. Black diamonds and triangles are sample points used in the accuracy analysis comparing predicted cumulative flow from the one-dimensional model to the FLOW-3D CFD model shown in (b) and (c), respectively.

Close modal

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 μm2 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.

FIG. 11.

Firing delay design space using the developed one-dimensional model with H = 8.5 μm. (a) Illustrates the cumulative flow resulting from varying τ and ξo,2 when ξo,1=0.2, L = 500 μm, and H = 8.5 μm. Dashed black line is the condition for no net flow and separates forward from reverse flow. Black diamonds and triangles are sample points used in the accuracy analysis comparing predicted cumulative flow from the one-dimensional model to the FLOW-3D CFD model shown in (b) and (c), respectively.

FIG. 11.

Firing delay design space using the developed one-dimensional model with H = 8.5 μm. (a) Illustrates the cumulative flow resulting from varying τ and ξo,2 when ξo,1=0.2, L = 500 μm, and H = 8.5 μm. Dashed black line is the condition for no net flow and separates forward from reverse flow. Black diamonds and triangles are sample points used in the accuracy analysis comparing predicted cumulative flow from the one-dimensional model to the FLOW-3D CFD model shown in (b) and (c), respectively.

Close modal

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.

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.

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.

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

A [m2]

channel cross-sectional area

L [m]

channel length

m [{0, 1}]

discrete model pressure model select parameter

po [Pa]

atmospheric pressure

pc [Pa]

postcollapse pressure spike

pv [Pa]

bubble pressure

pvr [Pa]

vapor saturation pressure at room temperature

p1b,2b [Pa]

bulk reservoir pressure

p1,2 [Pa]

pressure at channel entrance region

qc [kg m/s]

postcollapse momentum impulse

qo [kg m/s]

bubble strength

tc [s]

time at bubble collapse

xo [m]

resistor placement

xc [m]

location of bubble collapse

x1,2,3,4 [m]

location of liquid/vapor interface

κ [Pa s]

viscous dissipation coefficient

vc [m/s]

velocity at bubble collapse

ρ [kg/m3]

fluid density

τ [s]

resistor firing offset

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 

Vx(y,z)=ΔPμL16abn,m=0sin(p2n+1y)sin(q2m+1z)p2n+1q2m+1(p2n+12+q2m+12),
(A1)
Q=ΔPμL64abS1(a,b),
(A2)

where

p2n+1=π(2n+1)a,
(A3)
q2m+1=π(2m+1)b,
(A4)
S1(a,b)=n,m=01p2n+1q2m+1(p2n+12+q2m+12).
(A5)

The velocity profile can be re-written in terms of the average velocity, v = Q/A

Vx(y,z)=vab4S1n,m=0sin(p2n+1y)sin(q2m+1z)p2n+1q2m+1(p2n+12+q2m+12).
(A6)

Calculation of the viscous stress tensor and integrating over the wall perimeter gives the total viscous force which can be rewritten as F = κvL to agree with the one-dimensional model. In the case of a rectangular channel cross section26 

κ=abμS2(a,b)S1(a,b),
(A7)

where

S2(a,b)=n,m=01(p2n+12q2m+12).
(A8)
1.
S.
Hassan
and
X.
Zhang
, “
Design and fabrication of capillary-driven flow device for point-of-care diagnostics
,”
Biosensors
10
,
39
(
2020
).
2.
Q.
Shizhi
and
H.
Bau
, “
Magneto-hydrodynamics based microfluidics
,”
Mech. Res. Commun.
36
,
10
(
2009
).
3.
N.
Mishchuk
,
T.
Heldal
,
T.
Volden
,
J.
Auerswald
, and
H.
Knapp
, “
Micropump based on electroosmosis of the second kind
,”
Electrophoresis
30
,
3499
(
2009
).
4.
J.
Snyder
,
J.
Getpreecharsawas
,
D.
Fang
,
T.
Gaborski
,
C.
Striemer
,
P.
Fauchet
,
D.
Borkholder
, and
J.
McGrath
, “
High-performance, low-voltage electroosmotic pumps with molecularly thin silicon nanomembranes
,”
Proc. Nat. Acad. Sci. U. S. A.
110
,
18425
18430
(
2013
).
5.
K.
Vinayakumar
,
G.
Nadiger
,
V.
Shetty
,
S.
Dinesh
,
M.
Nayak
, and
K.
Rajanna
, “
Packaged peristaltic micropump for controlled drug delivery application
,”
Rev. Sci. Instrum.
88
,
015102
(
2017
).
6.
D.
Duffy
,
H.
Gillis
,
J.
Lin
,
N.
Sheppard
, and
G.
Kellogg
, “
Microfabricated centrifugal microfluidic systems: Characterization and multiple enzymatic assays
,”
Anal. Chem.
71
,
4669
(
1999
).
7.
V.
Gnyawali
,
M.
Saremi
,
M.
Kolios
, and
S.
Tsai
, “
Stable microfluidic flow focusing using hydrostatics
,”
Biomicrofluidics
11
,
034104
(
2017
).
8.
J.
Lake
,
K.
Heyde
, and
W.
Ruder
, “
Low-cost feedback-controlled syringe pressure pumps for microfluidics applications
,”
PLoS One
12
,
e0175089
(
2017
).
9.
M. I.
Mohammed
,
S.
Haswell
, and
I.
Gibson
, “
Lab-on-a-chip or chip-in-a-lab: Challenges of commercialization lost in translation
,”
Procedia Technology
20
,
54
–59 (
2015
), proceedings of The 1st International Design Technology Conference, DESTECH2015, Geelong.
10.
E.
Torniainen
,
A.
Govyadinov
,
D.
Markel
, and
P.
Kornilovitch
, “
Bubble-driven inertial micropump
,”
Phys. Fluids
24
,
122003
(
2012
).
11.
H.
Hoefemann
,
S.
Wadle
,
N.
Bakhtina
,
V.
Kondrashov
,
N.
Wangler
, and
R.
Zengerle
, “
Sorting and lysis of single cells by bubblejet technology
,”
Sens. Actuators, B
168
,
442
445
(
2012
).
12.
B.
Hayes
,
A.
Hayes
,
M.
Rolleston
,
A.
Ferreira
, and
J.
Kirsher
, “
Pulsatory mixing of laminar flow using bubble-driven micro-pumps
,” in
Proceedings of the ASME 2018 International Mechanical Engineering Congress and Exposition
(
2018
), Vol.
7
.
13.
E.
Ory
,
H.
Yuan
,
A.
Prosperetti
,
S.
Popinet
, and
S.
Zaleski
, “
Growth and collapse of a vapor bubble in a narrow tube
,”
Phys. Fluids
12
,
1268
(
2000
).
14.
Z.
Yin
and
A.
Prosperetti
, “‘
Blinking bubble’ micropump with microfabricated heaters
,”
J. Micromech. Microeng.
15
,
1683
(
2005
).
15.
M.
Einat
and
M.
Grajower
, “
Microboiling measurements of thermal-inkjet heaters
,”
J. Microelectromech. Syst.
19
,
391
(
2010
).
16.
A.
Govyadinov
,
P.
Kornilovitch
,
D.
Markel
, and
E.
Torniainen
, “
Single-pulse dynamics and flow rates of inertial micropumps
,”
Microfluid. Nanofluid.
20
,
73
(
2016
).
17.
E.
Sourtiji
and
Y.
Peles
, “
A micro-synthetic jet in a microchannel using bubble growth and collapse
,”
Appl. Therm. Eng.
160
,
114084
(
2019
).
18.
B.
Hayes
,
A.
Govyadinov
, and
P.
Kornilovitch
, “
Microfluidic switchboards with integrated inertial pumps
,”
Microfluid. Nanofluid.
22
,
15
(
2018
).
19.
P.
Kornilovitch
,
A.
Govyadinov
,
D.
Markel
, and
E.
Torniainen
, “
One-dimensional model of inertial pumping
,”
Phys. Rev. E
87
,
023012
(
2013
).
20.
H.
Yuan
and
A.
Prosperetti
, “
The pumping effect of growing and collapsing bubbles in a tube
,”
J. Micromech. Microeng.
9
,
402
413
(
1999
).
21.
J.
Zou
,
B.
Li
, and
C.
Ji
, “
Interactions between two oscillating bubbles in a rigid tube
,”
Exp. Therm. Fluid Sci.
61
,
105
(
2015
).
22.
C.
Hirt
and
B.
Nichols
, “
Volume of fluid (vof) method for the dynamics of free boundaries
,”
J. Comput. Phys.
39
,
201
225
(
1981
).
23.
C.
Borgnakke
and
R. E.
Sonntag
,
Fundamentals of Thermodynamics
, 8th ed. (
Wiley
,
1999
).
24.
O. E.
Ruiz
, “
CFD model of the thermal inkjet droplet ejection process
,”
in
Proceeding of Heat Transfer Summer Conference
(
2007
), Vol.
3
.
25.
T.
Theofanous
,
L.
Biasi
,
H.
Isbin
, and
H.
Fauske
, “
A theoretical study on bubble growth in constant and time-dependent pressure fields
,”
Chem. Eng. Sci.
24
,
885
897
(
1969
).
26.
S.
Timoshenko
and
J.
Goodier
,
Theory of Elasticity
, 3rd ed. (
McGaw-Hill, Inc
.,
1970
).