A new criterion of coalescence-induced microbubble detachment in three-dimensional microfluidic channel

Microbubble detachment is a common phenomenon in many industrial and medical processes, such as electrochemical reaction in microchannel,^1,2 thermal management of semiconductor,^3,4 self-cleaning for membrane biofouling,^5,6 and microbial colonization of dental surface,^7,8 to name a few. In general, external energy input and internal motion are two major mechanisms to drive the bubble detachment from a surface. External energy through heat transfer,⁹ acoustic effects,¹⁰ and light irradiation¹¹ causes the growing of microbubbles and further detaching from the surface. This type of detachment research has been on the phenomena of bubble detachment in specific industrial systems¹² and the optimization to balance between detachment and unclogging.^3,13 Whereas bubble detachment due to internal motion, e.g., microbubble coalescence,¹⁴ is more fundamental for exploring the underlying physics of bubble dynamics to support the design and fabrication of microfluidics with various applications. There have been efforts to reveal the effects of surface conditions,^15–18 initial conditions, and radius inequality^19,20 on the coalescence-induced droplet detachment. However, the mechanism of coalescence-induced microbubble detachment (CIMD) has been rarely addressed due to the fact that CIMD occurs extremely fast. The process of microbubble coalescence, from two (or more) touched parent microbubbles with their diameters from $10−3 to 10−6 m$ to eventually one single child bubble with a minimal surface area, takes a few hundreds of microseconds.²¹ A CIMD occurring at the last stage of the coalescence lasts only dozens of microseconds. Experimental investigation for such a fast activity requires sophisticated equipment, like an advanced fast camera with one million frames per second and microsecond exposure time. Therefore, experimentation of microbubble coalescence is usually more focused on the early stage of microbubble coalescence, i.e., neck growth.²² Whereas numerical simulation can be amenable to reveal the underlying physics related to the entire coalescence process, including CIMD through parameterization and classification.

The kinetic-based lattice Boltzmann method (LBM) has emerged for simulating a broad class of complex flows including pore-scale porous media flow,²³ multiphase/multicomponent flows,²⁴ combustion,²⁵ turbulence,²⁶ and others.²⁷ Comparing with the Navier–Stokes solvers, the main advantages of the LBM for this work are its amenability for modeling the intermolecular interactions at the two-phase interface to recover the appropriate multiphase dynamics without demanding computation cost and its suitability for scalable Graphics Processing Unit (GPU) parallelization²⁸ to achieve fast computation. In the past three decades, several multiphase models using the LBM have been developed, including the color fluid model,²⁹ the pseudo-potential model,³⁰ the interface tracking model,³¹ the free energy model,³² and the entropic LBM.³³ When a multiphase flow involves a large density ratio of liquid vs gas as well as a solid–liquid–gas interface, such as CIMD in the current study, the LBM models may become problematic. The major challenge is to overcome the numerical instability due to the occurrence of a parasitic current, a small-amplitude artificial/nonphysical velocity field arising from an imbalance of discretized forces across the interface. In general, in order to pursue a low parasitic current, either the accuracy of the contact angle in a certain region or the density ratio between two fluids has to be compromised. For example, insufficient accuracy at small, medium, and large contact angles has been reported when the pseudo-potential model,^34,35 the color fluid model,³⁶ and the free-energy based model^37,38 were employed, respectively. The interface tracking model³¹ can cover the full range of the contact angle but is limited to a low-density ratio of two fluids.^39,40 Detail discussions about these models are referred to in Huang's book.⁴¹ To overcome the drawbacks, the Lee group integrated the interface tracking and free energy models with the cubic polynomial boundary conditions.⁴² This approach has been continuously developed and refined in the last 10 years.^43,44 It has been demonstrated that this approach has eliminated the parasitic current without a compromise of inaccurate contact angle,^44,45 which well fits the need of the current study.

The current work is part of the continuous efforts^21,46,47 to reveal the underlying physics of microbubble coalescence to support the design and optimization of microfluidics where microbubbles are involved. Motivated by a recent experimental observation²¹ in the microbubble coalescence inside a polymer microfluidic gas generation device that CIMD occurs in some cases but not in others, we intend to explore the key factors that lead to CIMD and the corresponding mechanisms behind it. As the six experimental cases evidenced that CIMD occurs when the parent bubbles are equal-sized with a relatively large radius, we divided the study into three steps. First, we extend the in-house 2D GPU-parallelized LBM code to 3D flow involving the three-phase interaction among gas, liquid, and solid. Validations are carefully conducted via comparisons with analytical solutions and experimental results. Then, we perform a parametric study by simulating 12 cases. Together with the six experimental cases, we derived a new criterion for the occurrence of CIMD. At last, we explore the driving mechanism behind the CIMD and understand why the two conditions are needed.

The remainder of the paper is organized as follows. Section II formulates the lattice Boltzmann modeling for two-fluid flows and solid–fluid–gas interface. The experimental data and computational setup, validation, and numerical results are presented in Sec. III. Finally, Sec. IV provides a summary of the paper.

The formulation of the LBM for two-fluid flows on a solid surface consists of (1) governing equations, (2) lattice Boltzmann equations (LBEs), and (3) boundary conditions on a solid surface. The first parts have been presented in Ref. 24. Thus, we concisely introduce the main ideas and major equations here for the completion of this paper. The focus is on the introduction of the treatment of the solid surface to the two-fluid flow as a boundary condition. The two-fluid flow can be liquid–gas and liquid–liquid, involving bubbles and/or droplets. In our study, we refer to liquid–gas flows.

Two-fluid flows are governed by continuity, pressure evolution, and momentum equations. A diffuse interface is applied to separate gas and liquid as the modeling of the two-fluid flows; thus, the continuity equation can be written as the Cahn–Hilliard equation,⁴⁸ resulting in the following governing equations:

with $M(>0)$ the mobility, C the composition, μ the chemical potential, and p₁ the dynamic pressure.⁴⁹ The chemical potential is formulated as $μ=μ0−κ∇2C$ by minimizing the free surface energy on the two-fluid interface,

through the calculus of variations,⁵⁰ in which μ₀, E₀, and κ are the classical part of the chemical potential, bulk energy, and gradient parameter, respectively. In the vicinity of the critical point, van der Waals equation of state is simplified⁵¹ for the interface thickness and surface tension at equilibrium and the energy E₀ is formulated as $E0=βC2(C−1)2$⁠,⁵² with β a constant. Thus, $μ0=∂E0/∂C=2βC(C−1)(2C−1)$⁠. In an interface at equilibrium, the interface profile is $C(z)=0.5+0.5tanh(2z/D)$⁠, where z is the distance normal to the interface and D is the numerical interface thickness, which is chosen based on accuracy and stability. Given D and β, one can compute the gradient parameter $κ=βD2/8$ and the surface tension $σ=2κβ/6$⁠. The total pressure P is the sum of dynamic pressure p₁, the thermodynamic pressure p defined by $p=C∂E0/∂C−E0=βC2(C−1)(3C−1)$⁠, and the pressure caused by the curvature $−κC∇2C+κ|∇C|2/2$⁠, i.e., $P=p+p1−κC∇2C+κ|∇C|2/2$⁠.

The LBE (before the time discretization) for the two-fluid flow is given as follows:⁴⁴ 

where $F(=∇p−∇p1+C∇μ)$ is the intermolecular force between liquid and gas on the interface, $fα$ is the particle distribution function with discrete molecular velocity $eα$ along the α-th direction, and λ is the relaxation time related to the kinematic viscosity $ν=c2λ/3$⁠. The equilibrium distribution function is defined as $fαeq=ρωα[1+3(eα·u)/c2+9(eα·u)2/(2c4)−3u2/(2c2)]$⁠, where $ωα$ is the weight associated with a particular discretized velocity $eα$⁠, ρ and u are macroscopic density and velocity, respectively, and $c=δx/δt=1$ in lattice units (i.e., $δt=δx=1$⁠).

Defining a new particle distribution function $gα=fαc2/3+(p1−ρc2/3)Γα(0)$ in which $Γα(u)=fαeq/ρ$ and taking the total derivative $Dt=∂t+eα·∇$ of $gα$ and along with characteristics over the time step δt result in

where $∇MD$ and $∇CD$ are referred to as mixed difference approximation and central difference approximation, respectively,⁴⁴ and $τ(=λ/δt)$ is the non-dimensional relaxation time.

The hydrodynamic pressure and momentum are the zeroth and first-order moment of $g¯α$⁠, computed as $p1=∑g¯α+δt6u·∇CDρc2$ and $ρu=3c2∑eαg¯α−δt2C∇CDμ$⁠.

For the transformation of composition C, a second distribution function is introduced in a simple format of $hα=(C/ρ)fα$ and $hαeq=(C/ρ)fαeq$⁠. Similarly, taking the total derivative D_t of $hα$ yields

The composition C is the zeroth-order moment of $h¯α$ computed as $C=∑αh¯α+0.5δtM∇2μ$⁠. The density ρ and the dimensionless relaxation frequency (⁠$1/τ$⁠) are taken as linear functions of the composition by $ρ(C)=Cρ1+(1−C)ρ2$ and $1/τ(C)=C/τ1+(1−C)/τ2$⁠.

Take the solid boundary effects into consideration for lattice Boltzmann modeling of the gas–liquid gas flow. There are four different types of boundary conditions as follows:

1. Unknown particle distribution functions at boundary nodes are obtained from the bounce-back scheme, i.e., $gα(xs,t)=gα(xf,t)$ and $hα(xs,t)=hα(xf,t), xs$ is the node on the boundary, $xf=xs−eαδt$ is the node in the fluid.

2. The boundary for $∇ϕ$ is applied to prevent unphysical mass and momentum transfer through the boundary nodes. $ϕ$ is a macroscopic variable such as p₁, C, and μ. The no flux condition is presented as $eα·∇ϕ|s=0$⁠. After discretizing the no flux condition, it can be simplified as $ϕ(xs+eαδt)=ϕ(xs−eαδt), ϕ(xs+ 2eαδt)=ϕ(xs−2eαδt)$⁠, where the points $(xs+eαδ), (xs+ 2eαδ)$ are in the fluid domain, the points $(xs−eαδ),(xs− 2eαδ)$ are out of the fluid domain.

3. The boundary for $∇2μ$ ensures no mass flux in the normal direction of the solid wall, i.e.,

   $n·∇μ|s=0.$

   (8)
4. The boundary for $∇2C$ can be derived from minimizing the free surface energy

   $Ψs=∫S(ϕ0−ϕ1Cs+ϕ2Cs2−ϕ3Cs3+⋯)dS$

   (9)
   caused by the interactions between the liquid–gas interface and solid surface. Minimizing the total energy by calculus of variations⁵⁰ leads to

   $n·∇C|s=dΨs/dCs,$

   (10)
   where n is the unit vector in the normal direction of the solid wall. Retain higher-order terms in $Ψs$ can avoid numerical instability. Hence, the cubic boundary condition is selected⁴² in this work. To construct a cubic boundary condition, the parameters in Eq. (9) can be selected as $ϕ0=ϕ1=0, ϕ2=1/2ϕc$⁠, and $ϕ3=1/3ϕc$⁠, where $ϕc$ is a constant to be determined.⁴⁴ Thus, the condition for the solid wall, i.e., Eq. (10) can be simplified as

   $n·∇C|s=−ϕcκ(Cs−Cs2),$

   (11)

   where $ϕc=Ωc2βκ$⁠. Here Ω_c is the wetting potential which related to the contact angle α, i.e., $Ωc=cos α$⁠.

The experiment to study microbubble transport was carried out inside a polymer microfluidic gas generation device, fabricated by aligning and sequentially stacking/thermophoresis bonding multiple layers of patterned polystyrene (PS) film,⁵³ as shown in Fig. 1(a). Figure 1(b) schematizes the cross-sectional view of the device. Hydrogen peroxide (⁠$H2O2$⁠) enters the inlet and performs a catalytic decomposition producing water (⁠$H2O$⁠) and oxygen (O₂) on the layer of Pt catalyst. The oxygen bubbles travel from left to right, driven by imbalanced capillary force due to unequal channel heights, and then trapped at the gas vent through a membrane. Right before the gas vent, the inner walls of the microchannel change from hydrophilic to hydrophobic surfaces, whereas water, driven by the moving bubbles, leaves at the outlet. The coalescence event was captured by high-speed, high-resolution synchrotron x-ray microscopy at beamline 32 ID-B of the Advanced Photon Source at the Argonne National Laboratory from the cross-sectional view of the reaction microchannel of the microfluidic gas generator device. In this experiment, 30% hydrogen peroxide (⁠$H2O2$⁠) was used as a reactant solution. The catalytic reaction is triggered by contacting the Pt catalyst, and $H2O2$ produces water (⁠$H2O$⁠) and oxygen (O₂). Microbubbles are adhesion on the Pt catalyst with 36°. As shown in Fig. 1(c), the phase contrast image captured by the high-speed camera at a frame rate of 100 kHz (10 μs/frame) clearly shows the profiles of microbubbles and other substances with 2 μm pixel size; the height of the channel is $144 μm$⁠; the detailed description of the design, fabrication, and imaging acquisition of the experiment is referred to as published papers.^21,53

Eight archived videos of bubble transport recorded from the experiment are available for this study. We identified six cases clearly showing the unequal-sized bubble coalescence. We use ImageJ to measure the radii and locations of the initial bubbles and coalescing bubbles. The density and viscosity of oxygen and hydrogen peroxide are $ρO2=1.3 kg/m3, ρH2O2=1060 kg/m3, ηO2=1.92×10−5 kg/(m s),ηH2O2=1.06×10−3 kg/(m s)$⁠, respectively.

Concurrent numerical simulation is conducted in a 3D domain schematized in Fig. 2. Two-parent oxygen bubbles (shaded) are aligned along the Z-direction and placed on the X–Z plane at Y = 0, where the solid interface is located in the hydrogen peroxide solution. The radii and center locations of the father (larger) and mother (smaller) bubbles are r_F, r_M, O_F, and O_M, respectively. The coalesced child bubble (red circle) with minimal surface area sits at O_e with a radius of R_e. We define the size inequality γ as the radius ratio of father vs mother bubbles, i.e., $γ=rF/rM$⁠. The length of the channel is L, which is determined by the size of microbubbles and varies from $100 μm$ to $264 μm$⁠. And the width of the channel W is the same as the length fixed at $144 μm$⁠. The six cases identified from the experiment are seen in Table I. The equilibrium contact angle is set as 36° for all the cases.

We impose periodic boundary in x and z directions. In y-direction, the no-slip boundary condition is used on the top wall; thus, the unknown particle distribution functions at the wall nodes are obtained from the bounce-back scheme. On the bottom wall where the solid–gas–liquid interface locates, the macroscopic variables in the solid are assumed to take the mirror image of in the fluid reflecting no flux at the interface, i.e., $ϕ(xs+eαδt)=ϕ(xs−eαδt),ϕ(xs+2eαδt)=ϕ(xs−2eαδt)$⁠, which $ϕ$ can be any of the physical variables, including composition C, pressure p₁, and chemical potential μ.

To capture the 3D coalescence process with CIMD, which is time demanding, we extended the 2D LBM code (with no solid phase involved) developed in our previous study⁴⁶ to 3D with a solid interface and implemented Compute Unified Device Architecture (CUDA) parallelization on Intel Broadwell E5–2683 v3@ 2.30 GHz and NVIDIA Tesla P100 with optimization. Table II exhibits the GPU acceleration with five different spatial resolutions. Here we use Million Lattice Updates Per Second (MLUPS) to indicate the computation performance.

It is noticed that the acceleration is much less significant than those we have achieved in our other projects^23,54,55 due to the much more complicated LBEs in the current two-phase flow. Nevertheless, the GPU parallelism has significantly reduced the computation time. It has already accelerated more than 1400 times compared with serial code. In the majority resolution of this parametric study $220×120×220$⁠, the wall-clock time of each job is reduced from 43 days (without GPU parallelism) to 0.5 h (with GPU parallelism). It has already significantly enhanced our research capability for this study.

Convergence check to find an appropriate spatial resolution and verification to demonstrate the accuracy of the computation are done concurrently in what follows. As shown in Fig. 3, a half-sphere bubble (in gray) with radius $R=60 μm$ and contact angle $α=90°$ is initially sitting on the solid wetting surface. Due to the unbalanced capillary force, the bubble will evolve toward an equilibrium state with an expected contact angle $αeq=60°$ determined by the wetting potential Ω_c for $ϕc$ in Eq. (11). The simulated contact angle, corresponding to the tangent direction of a contour line (black line) with contour level C = 0.5 to the solid surface at the three-phase contact point, can be calculated by Eq. (12). The height of the highest liquid–gas interface a and the length of the contact line of the bubble on the solid surface b are measured when the bubble reaches the final equilibrium state,

We now vary the spatial resolution from 44 × 24 × 44 to 264 × 144 × 264 with six different levels in between and calculated the simulated α^eq at each resolution. In Table III, the accuracy (the fourth column) is expressed by the relative error of simulated to expected α^eq, whereas the convergence (the fifth column) represents the relative difference between the two successive resolutions. To balance the computational cost and accuracy need, we select the spatial resolution of 220 × 120 × 220, as both the accuracy and convergence are smaller than 1%, to conduct the remaining study unless otherwise indicated.

Next, we vary the expected α^eq by changing the wetting potential Ω_c for $ϕc$ in Eq. (11) and further demonstrate the accuracy of LBM simulation. Maintaining the same initial geometrical conditions and the physical parameters as those in Fig. 3 except the wetting potential, we vary the expected α^eq from 40° and 150° and simulate the α^eq correspondingly. In Fig. 4(a), the line shows the analytical relationship between the expected contact angle and the wetting potential Ω_c. The corresponding simulated equilibrium contact angles are marked as dots. The nearly identical contact angles demonstrate the reliability of the numerical simulation. Meanwhile, we plot the contour lines of $C=0.2,0.5,0.8$ vertical to the solid surface from the simulation in Figs. 4(b) and 4(c) for $αeq=150°$ and $αeq=45°$⁠, respectively. The analytical equilibrium bubble shape is well aligned with the simulated bubble shape in each case, further confirming the reliability of the free-energy based LBM with cubic boundary conditions.

From the six experimental cases shown in Table I, there are two distinct phenomena in the late stage of the bubble coalescence—with and without CIMD. The CIMD only occurs in case 5 and only be seen in the movie taken from the experiment. In the snapshots of the plane pictures, the CIMD is blocked by the small bubbles adjacent to the solid interface. To demonstrate the with/without CIMD behavior, we numerically simulated the complete processes of bubble coalescence for cases 5 and 6 with the initial parent bubble settings as close as possible to the experimental environments, respectively. Figure 5 shows the time evolution of bubble coalescence from two attached bubbles initially to a single child bubble with and without CIMD in (a) case 5 and (b) case 6, respectively. It is noted that the small bubbles seen in the images (bottom row in each case) have no contacts to the coalescing bubbles, and thus no influences on the coalescence process. The overlapping was due to the imaging view in the experiment. In each case, six snapshots at representative time instants for both experiment and computation show the corresponding coalescing bubbles in the simulation (top row) and experiment (bottom row). The simulation results with (a) and without (b) CIMD are in good agreement with the experimental measurement, which demonstrates the reliability of our computation model. Additionally, the uniqueness of case 5 includes nearly unit size inequality and a relatively large radius of the father bubble. Hence, hypothesize CIMD occurs when the two-parent bubbles are equal and with relatively large size. Due to the fast coalescence and uncontrollable bubble size in the experiment, we use the validated LBM simulation to test the hypothesis and derive a criterion, if it exists, for CIMD.

We extend Table I with six experiment cases to Table IV with a total of 18 cases by adding 12 computational cases. The 18 cases vary γ from 1 to 3, r_F from 19 μm to 62 μm, and $rF/H$ from 0.13 to 0.43 Among the 18 cases, cases 10, 14, and 17 ended with CIMD, of which the first two are from computation, and the third one was from the experiment. These three cases meet the conditions of the (nearly) unit γ and relatively large, $rF/H$ (0.27 and beyond). Those cases with the unit γ but relatively small, $rF/H$ (0.25 and below), i.e., cases 3 and 7, do not have CIMD, whereas the other cases with relatively large, $rF/H$ (0.27 and beyond), but non-unit γ, i.e., cases 11, 12, 13, and 15–18, do not develop CIMD either. Figure 6 shows the phase diagram of $rF/H$ vs γ in the ranges of $0.13−0.43$ vs $0−3$⁠, respectively. Cases with CIMD are marked as dots. It is confirmed that CIMD occurs when parent bubbles are (nearly) equal (⁠$γ≈1$⁠) with relatively large father radius (⁠$rF/H>0.27$⁠, i.e., $rF>40 μm$⁠). This criterion for CIMD is consistent with the observation of bubble detachment driven by heart transfer (external energy).⁵⁶ 

After the hypothesis is confirmed, a pertinent question about why CIMD occurs when parent bubbles are (nearly) equal with a relatively large father radius arises. In this section, we aim to address this question by exploring the underlying physics behind CIMD from the following two aspects.

In Fig. 7, we compare the time evolution of the instantaneous velocity vector fields on the vertical (X–Y) plane at the center of the bubble in the Z-direction at (1) early (top row), (2) intermediate (middle row), and (3) late (bottom row) stages of coalescence between unequal (left column, $γ=$ 1.94) and equal (right column, $γ=$ 1.0) initial parent bubbles. The color map represents the level of the velocity magnitude. The radius of the father bubble is the same, i.e., $rF=$ 48 μm or $rF/H=0.33$⁠, in the two cases. The distinct evolution of the velocity field illustrates why CIMD occurs in equal bubble coalescence but not in the unequal case. The coalescence process of the unequal-sized case (left column of Fig. 7) is described as follows. At the early stage (a-1), right after the neck bridge is formed, the father and mother bubbles have unbalanced capillary force due to their different curvatures, larger at the mother bubble side than the father bubble side. Two opposite but unequal flow streams squeeze at the neck area and induce a pair of unequal flow stream vertical to the neck interface, which tends to flatten the neck bridge. When the neck bridge becomes flattened at the intermediate stage (a-2), forming a single child bubble, the top neck bridge is disappeared, and the mother bubble is detached. Big kinetic energy is supplied to push the low part of the mother bubble toward the father bubble. At the late stage (a-3), when the child bubble oscillates till a minimal area reaches, the kinetic energy spreads in the entire field, and the bubble sticks on the solid interface. In the case of the equal-sized case (right column of Fig. 7), the evolution of the velocity vector field clearly illustrates why and how the CIMD is developed. At the early stage (b-1), due to balanced capillary forces from equal-sized bubbles, two balanced flow streams squeeze from the opposite direction and push the fluid leaving from the top and bottom neck bridge. At the intermediate stage (b-2), the top neck bridge disappears. Two strong horizontal flow streams continue to squeeze and push the fluid to leave vertically. Due to the blockage of the downward stream caused by a solid surface, the intensity of upward streams is promoted. A pair of vortices symmetric to the vertical centerline is formed. At the late stage (b-3), the pair of the symmetric vortices continues to drive the fluid moving upward, leading to the detachment of the bubble. Figure 8 exhibits the time evolution of the instantaneous vorticity fields on the same plane and at the same time stages. On the left, symmetric and paired vortices appear at the early stage (a-1). Still pairs of vortices are inclined along either northwest to southeast or southwest to northeast, parallel or vertical to the tangential direction of the neck bridge. At later stages (b) and (c), no mechanism drives the single bubble to move upward, whereas in the equal-sized case on the right, paired and symmetric vortices continue to drive the fluid moving vertically away from the gas–fluid interface all the time. Due to the blockage of the downward flow stream, the upward stream dominates and causes bubble detachment at the end of the coalescence process. Overall, the CIMD occurs due to the symmetric properties, including the capillary forces, kinetic energy, and vortices.

To understand why the CIMD occurs when the radius of the equal-sized parent bubbles is relatively large, we look into the development of the total kinetic energy for the coalescence process. As shown in Fig. 9, four equal-sized (⁠$γ=$ 1.0) cases are included with different radii, r_F = 48 μm, 42 μm, 36 μm, and 30 μm corresponding to case # 14, 10, 6, and 3 in Table IV, respectively. The total kinetic energy is normalized by the initial free surface energy contributed by the liquid–gas interface, Eq. (4) for $Ψb$⁠, and solid–liquid–gas interface, Eq. (9) for $Ψs$ and the evolution time is normalized by the inertial timescale $ti(=ρhrF3/σ)$⁠. The CIMD appears in the cases with the relatively large size cases (solid lines) but not in relatively small size cases (dashed lines). The kinetic energy development follows the same tendency in either with or without CIMD bubble coalescence. Initially, the two touched parent bubbles are static with zero kinetic energy. The coalescence is initiated due to the surface tension, converting the free surface energy to kinetic energy. At this early stage of coalescence, inertia dominates, resulting in an increase in kinetic energy till a peak at about the inertial timescale (t = t_i). Then, kinetic energy decays due to the viscous dissipation till the end of the coalescence. The difference between with and without CIMD is the amount of kinetic energy transferred from the initial free surface energy. Figure 9 implies that the CIMD occurs when the kinetic energy is large enough, i.e., more than 1.5% of the initial free energy. The detachment occurs when $2.5<t/ti<3.0$ in the large size cases. The large size cases can reserve enough kinetic energy to drive the detachment in the post-coalescence stage after the energy decay. This analysis explains why the CIMD does not occur in the small size cases (dashed lines in Fig. 9).

We have systematically studied the bubble coalescence in a 3D microfluidic channel with and without CIMD combining the experimental observation and numerical parameterization. A new criterion for the occurring of CIMD has been derived from 18 experimental and computational cases, varying the size inequality (the radius ratio of father vs mother bubble) and the radius of the father bubble. The in-house 2D GPU-parallelized LBM code for two-fluid flows has been extended to 3D two-fluid flows on a solid surface. The interaction among liquid, gas, and solid is modeled as the boundary conditions on the surface. The GPU acceleration enables the massive parameterization, from which the new criterion of CIMD is derived. Validations via comparisons with analytical solutions and experimental results have demonstrated the reliability of the computation. The new criterion is that when the two-parent bubbles are nearly equal (⁠$γ≈1$⁠) with a relatively large radius, e.g., $rF>40 μm,rF/H>0.27$⁠, a CIMD is expected at the end of the coalescence. The mechanisms behind this CIMD criterion are explored. When the two-parent bubbles are (nearly) equal, the two-fluid flow has symmetric properties along the horizontal direction, including the capillary forces, kinetic energy, vortices, and the flow is blocked at the solid surface, forming a flow streaming away from the solid surface. When the kinetic energy of the flow stream is large enough, the CIMD occurs. We will continue to study the various effects on the CIMD criterion derived in this work. The near future work includes the effects of initial on the CIMD. We expect the new criterion will impact the design and optimization of microfluidics in various applications.

This research was supported by the National Science Foundation under Grant No. 1264739, the National Natural Science Foundation of China No.11872062, Fundamental Research Funds for the Provincial Universities of Zhejiang No. 2020YW05, and the Graduate Scholarship from the Department of Mechanical and Energy Engineering at IUPUI. And this work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation Grant No. ACI-1053575.⁵⁷ 

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