Digital Sub-millimeter Bubble-Jets

Sub-millimeter bubble-jet is an uncommon phenomenon that shows great potential in fields such as drug delivery Rodríguez-Rodríguez et al. (2015), biochemical reactions Patel et al. (2021), and cell puncture Nakagawa et al. (2017). Many studies trigger cavitation through external excitations like ultrasound Lauterborn and Mettin (2023); Yusof et al. (2016), shock waves Deng et al. (2024); Huang et al. (2021), or lasers Chen et al. (2024); Agrež et al. (2024), thereby inducing sub-millimeter bubble-jet. However, due to the complex applied conditions, the stability and controllability of the jet are relatively too low to show clear practical prospects.Researches on bubble-jet above millimeters are common, typically focus on bubble emergence Brasz et al. (2018); Cheng et al. (2020); Berny et al. (2020), bubble coalescence Andredaki et al. (2021), and bubble stretching Séon and Antkowiak (2012); Li et al. (2019), or gravity-driven Cheng et al. (2020); Berny et al. (2020); Séon and Antkowiak (2012) and surface-wave focusing Li et al. (2019) bubble-jets, etc. Unlikely, these studies are macroscopically related to natural fields like oceanic atomization, aerosol generation, atmospheric bubble dynamics, cloud - formation processes, and air-pollution control Blanco-Rodríguez and Gordillo (2020); Lorenceau et al. (2004); Liger-Belair and Jeandet (2003); Cochran et al. (2017); Chingin et al. (2018), gravity often cannot be ignored or such phenomena occur in open areas with free boundaries. The sub-millimeter bubble-jets limit to the capillary scale ($l_{c}=\sqrt{\sigma/(\rho g)}\approx 2.7mm$), usually $Bo<1$ ($Bo=\rho gL^{2}/\sigma$). The gravitation is relatively fail to surface effects, and capillary force arouses to restrain gas-liquid deformation. Therefore, the sub-millimeter bubble-jet usually requires externalities stimulating as mentioned above. Ultrasound Nikolov and Wasan (2019); Tan et al. (2009) triggers bubble oscillation and cavitation through periodic pressure waves, leading to bubble expansion and rupture; Shock waves Luo and Niu (2019) drive the rapid compression and expansion through drastic pressure changes; Lasers Brown et al. (2011); Patrascioiu et al. (2014) rapidly heat the bubbles through the photothermal effect, generating rapid internal pressure differences. While in field of digital microfluidics, since usually $Re<1$, the bubble-jet phenomenon is rarely involved. Few research may be related, such as inward-rolling bubble-jets at relatively large Reynolds number ($Re\approx 400$) Pico et al. (2024), however, research on generating digital sub-millimeter bubble-jets has not reported yet.

This letter reports our inventions of generating digital sub-millimeter bubble-jet, as well as relevant findings and laws. We designed polymethylmethacrylate (PMMA) micro-chips with tapered channels Wang (2015, 2022), as shown in Fig.1(a). The unstable gas-liquid co-axial flow forms various flow patterns such as slug, jetting, dripping, and annular. The existence of the tapered zone breaks down the Galilean invariance to give two geometry variables: the taper angle $\alpha$ and the relative inner needle displacement $x$. In this work, $\alpha=5^{\circ},11^{\circ},17^{\circ},29^{\circ}$, and $x=500\mu\text{m},1000\mu\text{m}$, and $1500\mu\text{m}$, respectively. Their combinations yield 12 experimental data sets.The gas (argon) is injected as the dispersed phase from the middle needle, while the continuous phase (deionized water) is injected from the outer channel. The flow rates of the continuous phase (deionized water) and the dispersed phase (argon) provided by a peristaltic pump and a mass flow controller ranges of 5 $\sim$ 100 mL/min and 5 $\sim$ 50 mL/min. The interfacial tension is 72.8 mN/m, and the density and viscosity of water are 986.2 kg/m³ and 1.23 mPa·s (measured at room temperature).The dynamics of the bubbles and jets are illuminated by a high-intensity LED light, and snapshots are captured at 20,000 Hz (Phantom V611-16G-M). The dynamic process of jet excitation within monodisperse bubbles can be clearly captured, as shown in Fig.1(c). We also used the Volume of Fluid (VoF) to compare and verify the corresponding experimental results (Fig.1(b)).

Wang demonstrated that the co-axial liquid-liquid flow in the tapered region shows self-similarity and self-scaling. This was accomplished by introducing a similarity variable like the local virtual width, Wang (2022).

His data sets with different tapered configurations can be dimension-reduced to obtain a unified phase diagram and temporal and spatial scaling laws beyond physical mechanisms. We find that these self - similarity characteristics also hold in our gas-liquid system here.According to the corresponding conservation laws, the physical quantities are specially defined as follows: The superfacial velocity of the dispersed phase can be expressed as $u_{d}=Q_{d}/(\pi d^{2}/4)\cdot(W_{\text{local}}/W_{\text{out}})^{-1}$, while for the continuous phase is $u_{c}=Q_{c}/(W_{\text{local}}\cdot h-Q_{d}/u_{d})$. The overall velocity is defined as $u_{\text{TP}}=(Q_{d}+Q_{c})/(W_{\text{local}}\cdot h))$. The corresponding capillary number of the continuous phase, Weber number of the dispersed phase, and overall capillary number of the two phases are defined as

Under such definitions, since flow rate $Q_{c}$ and $Q_{d}$ ranges between 5 $\sim$ 100 mL/min and 5 $\sim$ 50 mL/min respectively, $u_{d}$ ranges 0.7 $\sim$ 4 m/s, and $u_{c}$ ranges 0.7 $\sim$ 7 m/s. Moreover, the range of $We_{d}$ is 0.01 $\sim$ 2, and that of $Ca_{c}$ is 0.01 $\sim$ 0.1.

After translating the data into self-scaling descriptions, we obtain a universal phase diagram in the $Ca_{c}\sim We_{d}$ space for all 12 different tapered configurations, as shown in Fig.2.Two types of phase charts are superimposed in Fig.2. One refers to flow patterns, and the other is the bubble-jet stimulating. Interestingly, these charts show obvious orthogonality and do not interfere, which indicating of independent mechanism.Since the large flow-rate pumps are used and dimensioned by $ml/min$ rather than $ml/h$, the Reynolds number moderately ranges $[20\sim 120]$ (Reynolds number is usually less than 1 in microfluidic research). Therefore, only two flow patterns, slug and annular, appear within this parameter range. The slug, known as Taylor flow, is regular periodic bubble flows, which relatively easy to identify (at the bottom of the phase diagram). While the annular is actually two types. One is the case where the gas-liquid flow is always stratified without disperse phase collapse (upper-left), where the inertia of the dispersed phase is relatively large while the viscosity of the continuous phase is relatively small. The other type of annular flow is formed during bubbles chase and coalescence, constantly merging and splitting into dynamic long bubbles (upper-right), where the inertia of the dispersed phase is very large, and the viscosity of the continuous phase is also relatively large. The dispersed phase will collapse.

And bubble-jets only occur when the dispersed-phase gas collapse.In fact, we found that there are three situations for bubble-jets to occur: 1) Dispersed-phase collapse near the nozzle in the tapered region; 2) Generated by the chasing and merging of two bubbles (coalescence in the annular snapshot in Fig.3(a)); 3) Inward-rolling of a long bubble during its movement Pico et al. (2024). This topic relates the first case.According to Fig.2, two types of threshold values exits for bubble-jets: 1) The bubble-jets excitation threshold; 2) The jet droplet emitting threshold. To generate an internal jet, a relatively large momentum transfer is required, and the shear force of the continuous phase needs to be large. However, the role of surface tension is rather contradictory. On the one hand, it suppresses surface deformation, on the other hand, it provides the original driving force for bubble-jet emitting.Actually, The threshold for the bubble-jets emitting has been discussed extensively in problems like free-surface emergence of bubble-jet, which depend on the critical Ohnesorge number ($Oh=\mu/\sqrt{\rho R\sigma}$, where $R$ is the initial bubble radius) and is related to the Bond number Gañán-Calvo (2017). As can be seen from our phase diagram Fig.2, this critical value also exists, manifested as the oblique dotted line $We_{d}\sim Ca_{c}^{-5.7}$ (with a lower intercept of $Ca_{c}\approx 0.03$). Similarly, the threshold value for the bubble-jet to generate droplets also exists, manifested as $We_{d}\sim Ca_{c}^{-4.2}$ (with an upper intercept of $Ca_{c}\approx 0.03$). We can define dimensionless numbers again to give the corresponding critical constants, but this will lead to problems such as inconsistent definitions and unclear physical concepts.Obviously, these existing thresholds are caused by the competitive relationships among various physical characteristics. Here, the dominant factor is the comparison among shear force, inertial force, and surface tension. However, this relationship can dominate the characteristics of the flow-pattern phase diagram. How does it transfer to the generation mechanism of bubble-jets? Obviously, the process of the bubble-jets generating droplets follows the well-studied Rayleigh-Plateau jet stabilities Yu et al. (2023); Anna (2016); Castillo-Orozco et al. (2015); Lhuissier et al. (2013). Therefore, the remaining most crucial issue is to clarify the initial bubble-jet excitation mechanisms and the laws of the intensity of bubble-jets.

From an observational perspective, the choices of what we can measure are of great significance. Regarding the excitation intensity of bubble-jet, following common practice, the tip velocity $u_{\text{jet}}$ can be selected to represent Séon and Antkowiak (2012); Ghabache et al. (2014); Gañán-Calvo (2017, 2018); Deike et al. (2018); Li et al. (2020). Since the period for bubble-jet emitting is extremely short, less than 0.2 milliseconds, which makes the tip velocity change significantly and easy to capture, and suitable for observing morphological evolution processes. We measured $u_{\text{jet}}$ ranges $1\sim 13m/s$. Noted that, different from the jet triggered at the free surface, our bubble is moving, $u_{\text{jet}}$ we refer to is relative to the overall movement of the bubble. Thus, the actual jet velocity ranges approximately $3\sim 25m/s$, which the relative velocity $u_{\text{jet}}$ can better reflect the characteristics of the jet’s excitation intensity. Therefore, macroscopically, we can write the expression,

The dependent variables of $u_{\text{jet}}$ fall into three categories: physical property parameters, geometric parameters, and flow control parameters.This expression makes no difference from cavernous descriptions of digital bubble generating Wang (2022). Though representing the characteristics of the tapered region with the similarity variable $W_{\text{local}}$ has achieved significant dimensionality reduction, there are still quite a number of variables. As a result, the experimental cost is too high to withstand. Of course, many parameters can be regarded as constants, but doing so will lose the generalization of the research.

Let’s first analyze the bubble-jet emitting process, as shown in Fig.3(a). Through observation, we can find the following: 1) The break-up of the dispersed-phase occurs in the tapered region, where presents a “front-end stretching and rear-end squeezing” accelerated stretching flow Wang (2015). 2) The bubble-jet emitting of slug and annular are similar or identical. Thus, macroscopic bubble characteristics do not influence the mechanism and intensity of bubble-jet emitting. This explains why there is orthogonality between the flow-pattern map and the jet-emitting chart in Fig.2. 3) During the processes of collapse, capillary-wave propagation, focusing, and bubble-jet formation (the red line in Fig.3(a)), the interface vibration only takes place in the tail region of the bubble. These observations suggest that the action, formation, and triggering of bubble-jet should be a local process.If we develop a local model, it should be reduced to Fig.3(b). We only need to consider the region from the jet break-up point (Neck point) to the point of the bubble’s maximum diameter (Crest point), while the influence of the bubble’s head (the parts of slug and annular drawn with dotted lines) can be ignored. Here, we’ve identified another measurable quantity, $L_{\text{cone}}$, as seen in Fig.3(b). Through this intermediate variable $L_{\text{cone}}$, we may bridge the overall macroscopic bubble behavior with the local microscope of jet emitting.

$L_{\text{cone}}$ represents the extreme stretching state of the bubble tail and serves as the initial configuration for onset of ejections, thus it holds great significance. To find mechanisms in bubble tail stretching, we will look into the relevant energy processes. Taking into account the previous analysis, the substantial stretching of the bubble is associated with the structure of the tapered region for the “front-end stretching and rear-end squeezing” effects Wang (2015). Unlike the situation of a freely breaking jet, the stretching in this non-uniform shear flow field is much more pronounced. The work $W$ done by the fluid on the gas-liquid interface can be expressed as,

Among them, the work $W_{sq}$ done on the gas-liquid interface by “rear-end squeezing” should come from the work done in the direction of interface deformation by the normal pressure difference $\Delta p\sim\sigma\kappa+\rho_{c}aW_{\text{local}}$ across the interface and the shear stress $\tau\sim\mu_{c}u/W_{\text{local}}$ along the interface tangent. They can be expressed as $W_{sq1}=\int_{S}\Delta p(\mathbf{u}\cdot\mathbf{n})\mathrm{d}A$ and $W_{sq2}=\int_{S}\tau(\mathbf{u}\cdot\mathbf{n})\mathrm{d}A$ respectively. Here, $\kappa\sim W_{\text{local}}$ is the interface curvature, $u\sim(Q_{d}+Q_{c})/W_{\text{local}}^{2}$ is the flow field velocity, $a=du/dt\sim u^{2}/W_{\text{local}}$ is the fluid acceleration caused by tapering, and $S\sim W_{\text{local}}^{2}$ is the area of the bubble-tail interface on which work is done.The work $W_{st}$ on the interface of the bubble head described by the “front-end stretching” effect is essentially the same as $W_{sq}$, so the work expressions are similar, with the difference lying in the specific flow-field characteristics. However, since the bubble-jet we observed exhibits the characteristics of local flow behavior, and the occurrence time of the entire bubble-jet process is much shorter than the flow time, the work of the overall “front-end stretching” of the bubble should be negligible in the triggering process of the bubble-jet, unless the bubble-jet occurs in the case of relatively small droplets (not belonging to the slug and annular flow patterns). Then we have $W_{st}\approx 0$, and thus $W\approx W_{sq}=f(Q_{d}+Q_{c},W_{\text{local}};\rho_{c},\mu_{c},\sigma)$.As an observable quantity representing the result of work-done, the expression of the intermediate variable $L_{\text{cone}}$ should also contain the same parameter variables. Therefore, we write $L_{\text{cone}}$ as the following formula,

Taking $W_{\text{local}}$, $\sigma$, and $Q_{d}+Q_{c}$ as independent characteristic variables, we obtain the dimensionless expression of Eq.5, that is $L_{\text{cone}}/W_{\text{local}}\sim\Phi(Oh_{\text{local}},Ca_{\text{local}})$. Among them, $Oh_{\text{local}}=\mu_{c}/(\rho_{c}\sigma W_{\text{local}})^{0.5}$ and $Ca_{\text{local}}=\mu_{c}(Q_{d}+Q_{c})/(\sigma\cdot W_{\text{local}}^{2})$. These two dimensionless numbers jointly characterize the local behavior of bubble stretching. $Oh_{\text{local}}$ emphasizes the competitive relationship between viscous force, inertial force, and surface tension, while $Ca_{\text{local}}$ describes the competitive relationship between fluid inertia and surface tension. They play a crucial role in the morphological deformation of the bubble, the jet velocity, and the final jet intensity. Through experimental data analysis, we obtain the relationship: $L_{\text{cone}}/W_{\text{local}}\sim Oh_{\text{local}}^{-2}\cdot Ca_{\text{local}}^{1.6}$. Restoring it to the form of original variables, we get

Formula 6 is in excellent agreement with our experimental data. For details, you can refer to the supplementary materials. Now that we’ve clearly analyzed the initial morphological characteristics in the onset of ejections process, let’s move on to discuss the intensity and mechanism of bubble-jet. Since the bubble-jet emitting process is a transformation from potential energy to kinetic energy, this transformation mainly involves two parts. One is the restoring potential energy of the stretched bubble, and the other is the energy obtained by the jet from the change in inertia of the two-phase fluid. Therefore, we have

Among them, $K_{\text{in}}=\int_{V}1/2\rho_{c}u_{\text{jet}}^{2}dV$ represents the kinetic energy of the jet, where $u_{\text{jet}}$ is the field velocity which is impractical for experiments. In fact, we usually substitute the field velocity in the integral with measured tip velocity of the jet, thus $K_{in}\sim\rho_{c}u_{\text{jet}}^{2}$. Since we are discussing the initial energy balance of the bubble-jet here, the intermediate variable $L_{\text{cone}}$ is an excellent local spatial characteristic quantity to reflect the energy situation during the bubble-jet triggering stage. Therefore, the capillary restoration energy is $\delta E_{ca}\sim\sigma L_{\text{cone}}^{2}$, and the inertial energy change is $\delta E_{in}\sim\rho_{c}u_{\text{TP}}^{2}$. Because the time of bubble rebound is shorter compared to the acceleration effect in the tapered region of the flow field, the viscous dissipation $E_{dis}\sim\mu_{c}(u_{\text{TP}}/L_{\text{cone}})^{2}$ should be taken into account during the rebound process and the bubble-jet emitting process. In this way, from Eq.7, the expression of $u_{\text{jet}}$ can be written as follows,

Select $\rho_{c}$, $\sigma$, and $u_{\text{TP}}$ as independent scaling variables, we derive $u_{\text{jet}}/u_{\text{TP}}\sim f(We_{\text{TP}},Ca_{\text{TP}})$. Among them, $We_{\text{TP}}=\rho_{c}L_{\text{cone}}u_{\text{TP}}^{2}/\sigma$ represents the relative importance between the inertial force and the surface tension determined by $u_{\text{jet}}$ and $L_{\text{cone}}$ during the local similarity. Through experimental data analysis, we obtain the relationship $u_{\text{jet}}/u_{\text{TP}}\sim We_{\text{TP}}\cdot Ca_{\text{TP}}^{-2.6}$. Finally, when expanded into the primitive variables, we have

It is evident that the intermediate variable $L_{\text{cone}}$ exhibits a remarkable connective effect. Upon substituting Eq. 6 into Eq. 9, the final relationship of the bubble-jet velocity $u_{\text{jet}}$ can be expressed as,

Thus, the global expression of $u_{\text{jet}}$ in Eq.3 is solved, and answers the jet intensity problem. Plot Eq.10 in Fig.4, it can be seen that the 12 datasets are basically on the formula (black straight line),which is quite accurate in characterizing the bubble-jet intensity features.

Eq.10 is fascinating for: 1) Incorporating inertial, surface tension, viscous, and the tapered geometry effects, which is consistent with viewpoints that the magnitudes of various mechanical quantities at the bubble-jet triggering stage are similar Gañán-Calvo (2017); Ismail et al. (2018). 2) Comparing $W_{\text{local}}^{1.6}$ and $H^{0.4}$, the effect of the tapered scale is to the fourth power of the pipe width, indicating that the fluid acceleration has significant influences on jet emitting, agreeing with literature Wang (2015, 2022). 3) If transformed, we can have $u_{\text{jet}}\propto Re^{3}/La$, where $La=\frac{\rho U^{2}L}{\mu\sigma}$, which shows fluid inertia is crucial in bubble-jet occurrence; also for bubble-jet digital emitting not being discovered in microfluidic researches due to small $Re$. 4) Surface tension is still positive, which is basically consistent with works on bubble emergence, bubble chasing, etc. Brasz et al. (2018); Berny et al. (2020); Andredaki et al. (2021). 5) Note that the formula holds for all the tapered configurations, and its generality can be ensured, that is, the tapered gas-liquid flow maintains self-scaling flow structures. 6) Such regular exponents from experimental data under compound physical mechanisms is rare, which indicates that the competing among various mechanisms is evenly matched, acting clearly, so the data can be well-posed for physical laws. 7) Due to the limits of instruments in ranges ($We_{d}\sim[0.01,2]$, $Ca_{c}\sim[0.004,0.1]$), outer border of the bubble-jetting is open, and critical boundary is unclear, when bubble-jet velocity of magnitude of $O(10^{2})$ appears in our simulations, which leads to another interesting topic.

Another thing that needs to be clarified is the locality issue. For the local model in Fig.3(b), it draws the forward capillary wave propagating. So how is the local nature of the model ensured? Three clues: 1) According to free space bubble-jetting, where capillary waves move towards the bubble bottom and focusing, called inertial-capillary focusing Ghabache et al. (2014); Andredaki et al. (2021); Berny et al. (2020), our waves should meet at the bubble head and excite jets. But here the bubble-jetting occurs at the bubble tail, so the jet must be triggered by reflected waves. When the capillary wave reaches crest points, it encounters the boundary resistance and quickly reflects, then focuses at the channel center and triggers the jet. Therefore, the reflected wave is dominant. 2) Capillary wave speed $c=\sqrt{(\sigma/\rho)k}$($k$ wave wavelength $\sim 1mm$), $\sim 0.3m/s$. The time for wave to reflect and focus after reaching the crest point is $W_{local}/(2c)$ or $h/(2c)$, $\sim 1\mu s$. For the forward wave, this time is insufficient to travel out of tail region. 3) In terms of intensity, the capillary wave may be relatively weak to give limited energy. But, for exciting a fast jet, the role of the inertial wave is large, which is shown in Eq.7 and Eq.10. Usually, the inertial wave speed is much greater than that of the capillary wave, which shorts the triggering time, also the chance of forward wave leaving.

This work was supported by the National Natural Science Foundation of China - Yunnan Joint Fund Key Project, U2002214. The authors would like to express their sincere gratitude to Professors Jin-Song Zhang, Guo-Hui Hu, Zhe-Wei Zhou, and Qian Xu for valuable assistance.