Investigations on the near-wall bubble dynamic behaviors in a diverging channel

Bubble- and microbubble-based processes have wide applications across various fundamental and industrial fields, including propulsion systems, petroleum exploitation, chemical engineering, nuclear reactors, wastewater treatment, and mineral flotation.^1–6 The diverging channel is a common geometric design in gas–liquid two-phase flow systems, found in components such as flowmeters, throttle valves, pumps, propulsion devices, bubble or cavitation generators, and pipe connections.^1,2,7–10 These applications frequently involve bubble dynamics in diverging channels where phenomena such as interphase energy or mass transfer and phase distribution are crucial. Studying bubbly flow through diverging channels is essential for developing more efficient industrial equipment and enhancing our understanding of complex physical phenomena.² 

In recent years, the study of bubble generators has increasingly emphasized engineering applications and design optimization. This research primarily focuses on enhancing bubble generation structures, evaluating device performance, and optimizing mass transfer applications. Innovative designs, including swirl-type,¹¹ self-priming,¹² and dual Venturi structures,¹³ enhance turbulence and significantly boost the production efficiency of microbubbles. These advanced structures facilitate the precise control of bubble size and distribution by intensifying fluid shear and turbulence. Performance evaluation focuses on key metrics such as bubble size distribution, generation frequency, oxygen mass transfer rate, and energy efficiency. Assessing mass transfer efficiency involves analyzing oxygen mass transfer coefficients, gas solubility, and device energy efficiency.^14,15 Future research on bubble generators is expected to prioritize reducing energy consumption, enhancing bubble generation precision, and achieving real-time monitoring and adjustment of bubble characteristics through intelligent and automated design.¹⁶ Accurately predicting bubble dynamics, including transport, breakup, and coalescence, will provide a solid theoretical foundation for optimizing device design and engineering applications. This research direction aims to advance bubble generation technology toward greater efficiency and intelligence through theoretical and experimental approaches.

In systems featuring Venturi or nozzle geometries with diverging sections, variables such as pressure variation, bubble acceleration and deceleration, coalescence and fragmentation, and bubble size distribution are highly coupled.^17,18 Paladino and Maliska¹⁹ focused on accelerated bubbly flow through Venturis, particularly in diverging sections, and found that added-mass forces were critical for predicting differential pressure, while lift forces significantly influenced void fraction distribution. Ahmadpour et al.² investigated bubbly flow through expansions and contractions using an Eulerian two-phase flow model and identified significant increases in drag and added-mass forces in expansions. They applied the standard k–ε model for simulations. However, for more accurate modeling of flow separation, many researchers have used more feasible turbulence models, such as the Shear Stress Transport (SST) k–ω, v²-f (four-equation model), or RNG k–ε.^20–22 These studies underscore the importance of accurately modeling the forces acting on bubbles within diverging channels, particularly in highly dynamic acceleration and deceleration regions.

To analyze the bubble dynamics in accelerated flow fields, many studies have focused on the motion of individual bubbles in nozzles or Venturis. Previous research has demonstrated that bubble-induced turbulence has a relatively minimal impact on the mean motion and pressure of bubbles and the surrounding liquid.^23–25 Consequently, many subsequent studies have overlooked turbulence effects, particularly in scenarios involving clean water or small bubbles. However, turbulence remains significant in situations involving contaminated water or larger bubbles. Kuo and Wallis²⁴ investigated the movement of individual bubbles through a converging–diverging nozzle, using a one-dimensional model to predict the bubble velocity. Their model accounted for steady and unsteady viscous drag, added mass, and pressure gradient forces. While their experimental data generally supported their predictions, they were unable to draw definitive conclusions regarding the added-mass coefficient. Kowe et al.²⁵ introduced a method of replacing the mean liquid velocity with a space-averaged velocity field, considering the superficial liquid velocity and void fraction to account for the added-mass effect in bubbly mixtures. Soubiran and Sherwood⁷ expanded on this by developing a nondimensionalized bubble motion balance, focusing on the steady drag and added-mass coefficients in an unbounded fluid. Although theoretical models incorporating these forces provide valuable approximations, they often fail to capture the full spectrum of bubble behaviors, especially under transient or turbulent conditions. This indicates a need for further refinement of models that consider dynamic interactions between bubbles and the surrounding fluid, potentially through the use of advanced experimental and computational techniques.

van der Geld et al.²⁶ explored single tap-water bubbles in a rectangular cross section Venturi, employing a modified drag model where the drag coefficient depended on both the bubble Reynolds number (Re_B) and the acceleration number (Ac). Through Laser Doppler Anemometry (LDA), they found that acceleration significantly impacted drag on moving bubbles and that wall effects were also considered in their model. Their predictions matched well with experimental data, particularly for small bubbles (less than 1 mm), where the added-mass coefficient was consistently 0.5. This stability in smaller bubbles highlights the dominance of acceleration effects. Simcik et al.²⁷ examined the Volume of Fluid (VOF) method for evaluating the added-mass coefficient in different bubble scenarios, such as ellipsoidal, spherical cap bubbles, and those near walls or in clusters. Their analysis, however, focused on scenarios where drag could be neglected, leaving some bubble dynamics underexplored. Song et al.²⁸ studied bubble motion and breakup in the diverging section of a Venturi, emphasizing that local flow structures significantly influence the bubble behavior. Using Particle Image Velocimetry (PIV),²⁹ they provided insights into the flow environments impacting bubble dynamics. Zhao et al.^8,9 investigated bubble deceleration near the Venturi wall, observing significant deceleration that prolonged bubble residence in turbulent areas, thus enhancing breakup efficiency. Despite these observations, the mechanisms behind this deceleration remained unclear.

Building on our previous research advancements and the analysis of bubble forces by predecessors, this study further elucidates the interaction mechanisms between bubbles and liquids through a combination of visualization experiments and numerical simulations. Thus, the main objectives of the present work are to analyze the movement characteristics and causes of the dramatic deceleration of individual bubbles as they enter a diverging channel near the wall. A visualization experiment was conducted to observe bubble motion in tap water, concentrating on deceleration behavior near the wall in the diverging section. To better understand the velocity and acceleration relationship between the liquid and the bubble, it is crucial to accurately predict the bubble’s velocity. This involves incorporating a source term in the momentum equation based on the theory by van der Geld et al.²⁶ By integrating the visualization experimental data with numerical simulations, this work aims to clarify the roles of forces acting on individual bubbles during deceleration.

The experimental setup, depicted in Fig. 1, includes water supply, air injection, and data acquisition. A horizontally mounted PMMA Venturi channel serves as the primary test section. Air is injected into the water stream through a 1 mm diameter orifice at the throat section. The high-velocity liquid flow generates bubbles that propagate along the channel wall. The air flow rate is varied from 0.008 to 0.08 m³/h, while the water flow rate ranges from 2.4 to 6.0 m³/h. Mass flowmeters, with ±0.5% and ±0.1% accuracy, measure the air and water flow rates, respectively. Figure 2 details the Venturi setup and dimensions, designed to minimize refraction errors during observations for accurate flow visualization. The diverging section has a 20.0° expansion angle, with inlet and outlet widths of 25 and 50 mm, respectively, maintaining a height of 10 mm throughout. Static pressure is measured at the Venturi’s inlet, throat, and outlet using three pressure sensors (±0.04% accuracy), enabling detailed analysis of pressure distribution. High-speed imaging is performed with a FASTCAM Mini UX100 camera, capturing 4000–10 000 fps, facilitating detailed observation of rapid bubble formation, coalescence, and breakup processes, and offering insights into the micro-scale physics of gas–liquid interactions in accelerating and decelerating flows.

Considering the interference from neighboring bubble wakes, where the wake length ranges from ∼0.9 to 3.3 times of bubble diameter (D_B) within the bubble Reynolds number range of 3.6–44 (obtained from simulation results), approximately calculating by O(0.5D_BRe_B^0.5),³⁰ we selected bubbles that are at least 3D_B away from their nearest neighbors for the collection of bubble dynamics parameters and for the subsequent studies on bubble transport characteristics in this research.

A Digital Image Analysis (DIA) method was developed using MATLAB to measure the bubble motion parameters from high-speed camera images. The DIA process, illustrated in Fig. 3, encompasses three main steps: image segmentation and enhancement, bubble edge detection and identification, and bubble parameter calculation. In the initial step, the original image is divided into several blocks to enhance measurement accuracy. Each segmented image is filtered using a median filter³¹ to remove noise, followed by sharpening with top-hat and bottom-hat filters³¹ to improve image quality. In the second step, the target bubble is isolated from the background using gray level histograms according to Otsu’s method.³² The Canny algorithm³³ is then employed to detect bubble edges, generating a binary image. The identified bubbles are subsequently filled with pixels to complete their shapes. Finally, bubble parameters are calculated from the binary image. This image analysis technique is detailed in the work by Zhao et al.^9,31 Notably, when bubbles merge, the watershed approach⁹ is applied to separate them. Only individual bubbles, maintaining a distance greater than the size of the bubbles (typically 1.7–3.8 mm), are considered in the analysis.

Figure 5 presents the experimental errors in bubble size and position measurements. The image resolution was set to 50 pixels/mm to capture the bubble deceleration process. To enhance bubble edge detection accuracy, local thresholding and the Canny algorithm were employed in the DIA method. The uncertainty in bubble edge detection is typically less than 1.0 pixel. Experimental results indicate uncertainties of ∼5.0% and 3.0% in bubble volume and surface area measurements, respectively. This leads to a 3.1% uncertainty in the SMD of the bubble, as shown in Fig. 5(a). The uncertainty in determining the bubble’s centroid is less than 1.4 pixels using the DIA method. The high-speed camera introduces a 10⁻⁶ s system error, and with a minimum frame rate of 4000 fps over 1000 frames, the time uncertainty is 0.45%. Consequently, the uncertainties in bubble position, velocity, and acceleration between adjacent frames are 2.8% [Fig. 5(b)], 4.0%, and 5.7%, respectively.

To better understand the flow field in the divergent channel, CFD were conducted using the CFD code ANSYS Fluent. The rationality of 2D modeling for studying bubble rise in large-scale stilling pools has been well established. For example, Gaston et al.³⁶ investigated cavitation bubbles in a Venturi flow at the centerline with 1D bubble dynamic equations, and Li et al.³⁷ simulated ventilated bubbles in a micro-scale Venturi using the Volume of Fluid (VOF) model. In the current study, the bubbles deviate from the centerline of the narrow channel, necessitating a 3D modeling approach, as shown in Fig. 6(a). A structured meshing scheme was employed. To accurately capture the transitional flow, the SST k–ω turbulent model was used. The near-wall grids were further refined with an initial grid point spacing (y⁺ < 1) for more precise near-wall calculations. It is crucial to note that fine mesh schemes are required for y⁺ calculations near the wall and for the VOF model. However, the cells near the wall are often too narrow to meet the interface capturing requirements of the VOF model. To address this, an adaptive-boundary method was developed, as illustrated in Fig. 6(a). This approach enables precise calculations for both wall conditions and the VOF model.

Adaptive mesh technology was utilized to capture clear interfaces, which is critical for calculating slip velocity and interactions between the liquid and the bubble. At the velocity inlet, fully developed flow velocity profiles were defined using a User Defined Function (UDF) to achieve high computational efficiency, as depicted in Fig. 6(b). The outlet boundary was set as a pressure outlet, based on experimental measurements.

In defining the drag force model with UDS, the instantaneous bubble interface position was captured, and the VOF model’s prediction deviations were corrected by introducing source terms at the interface. This approach yielded more accurate information on bubble transport and the flow field around it without disturbing the liquid phase flow field.³⁹ 

Second-order upwind discretization is employed for momentum equations, PRESTO! for pressure, and first-order upwind for turbulent quantities. Table I summarizes the UDF and UDS implemented in this work. Figure 8 shows the simulation verifications for both cell number (N_cell) and pressure drop (ΔP), defined as the differential pressure between the inlet and outlet sensors. A mesh scheme comprising ∼6.3 × 10⁶ cells was adopted for this study, as further increases in cell number yielded negligible improvements in simulation accuracy. The numerical predictions demonstrate good agreement with experimental measurements.

Figures 9–11 present the simulation results for bubble deformation, trajectory, and velocity, respectively, for Q_L = 4.54 m³/h. Discrepancies between the conventional VOF model and experimental results are noted. The conventional VOF model, utilizing the k–ε turbulence model without momentum source terms, underestimates lateral forces on the bubble, resulting in a trajectory further from the wall and a higher velocity. This model fails to accurately capture the flow structure, leading to significant differences in bubble shape compared to experimental data. The modified VOF simulation considers the source term S_Φ in the momentum equation [Eq. (8)], incorporates our experimental findings, and employs mesh adaptation and the SST k–ω model. The momentum source term is applied only during periods of dramatic deceleration and is deactivated when this deceleration ceases, as per experimental observations.

The results indicate that the modified simulation predicts more closely results with experimental data than the conventional VOF model. This approach better reflects bubble dynamics in complex flow fields, accurately capturing the interactions between the bubble, the surrounding fluid, and the channel geometry. The enhanced predictive capability of this method provides crucial insights into the factors contributing to sudden bubble deceleration in divergent channels.

Accurate prediction of the continuous phase flow field distribution is crucial for analyzing bubble forces. To verify the continuous phase flow field calculation, this study simulated a 2D turbulent flow in a diffuser based on the experimental data of Obi et al. (1993)⁴⁶ using the same numerical simulation scheme. The geometry and mesh of the diffuser are shown in Fig. 12(a), with detailed parameters referenced from Ref. 46. The inlet bulk Reynolds number Re_bulk,in is 9000, similar in order of magnitude to the continuous phase Reynolds number in this study. Figures 12(b) and 12(c) demonstrate that both the v²-f and SST k–ω models can effectively predict the mainstream and separation flow structures in the diverging channel, with the latter showing better agreement with experimental data. This validates that the adopted turbulence models and computational schemes can accurately predict the velocity distribution of the continuous phase flow field.

The detailed observation of bubble behavior, including shape deformation and breakup, is critical for understanding and optimizing multiphase flow systems. As air phase enters the high-speed water flow in the throat region near the wall, it forms discrete bubbles, as shown in Fig. 13. These bubbles maintain a relatively steady shape and size in the throat until they reach the diverging section, while the higher Q_L and lower Q_G lead to smaller bubbles forming in the throat. When the flow rates exceed a certain threshold, the bubbly flow transitions into a plug flow. Upon entering the diverging section, the bubbles start deforming into ellipsoidal shapes and eventually break up at specific locations.³¹ 

Figure 14 shows the overall bubbly flow through the diverging section. Prior to entering the diverging section, bubbles keep approximate spheroidal or ellipsoidal shapes. Upon entering the diverging section, significant deformation occurs, typically resulting in elongated shapes inclined at specific angles to the mainstream flow (x-axis). Higher Q_L results in more pronounced deformation. Furthermore, the distance between adjacent bubbles decreases gradually due to a reduction in velocity. The decelerated bubbles often stack together rather than spreading to the channel center, attributed to the high-speed mainstream flow at the channel’s center. As a result, the bubbles break up into two or more smaller bubbles in the strong turbulent flow.

Figure 15 presents the simulation results of single-phase liquid flow structures in the diverging section, including velocity distribution, streamlines, static pressure, and turbulent kinetic energy (TKE). A local reduced pressure region is generated at the junction of the throat and diverging section walls (x = 0 mm, y = 0 mm), as indicated by the circled area in Fig. 15(c). This pressure drop induces the streamlines to suddenly swerve along the divergent wall, known as the Coanda effect.⁴⁷ Analyzing the bubble motion, indicated by blue arrows in Fig. 13, reveals that a large-scale vortex structure is evident in the diverging section [Fig. 15(b)]. On comparing Fig. 14 with the distribution of TKE in Fig. 15(d), it can be found that the vortex traps bubbles in the strong turbulent flow for extended periods, leading to complete fragmentation.⁴⁸ This suggests that the force constraining bubbles near the wall, namely lift, plays a significant role in the efficient preparation of bubbles by the bubble generator. According to bubble breakup mechanisms,⁴⁹ viscous shear force and turbulent fluctuations are the main mechanisms driving bubble breakup.⁵⁰ The breakup process is not a standalone phenomenon but is associated with the complete history of interactions between the liquid and the bubble prior to breakup. Studying the movement characteristics of individual bubbles provides valuable insights into bubbles behaviors before breakup.

Bubble trajectories are determined by capturing the centroid and displacement of bubbles using the DIA method. In experiments, more than 15 bubbles within the SMD range of 1.6–2.8 mm were measured for each selected case of Q_L. Figure 16(a) illustrates the trajectories of individual bubbles from the throat to the diverging section near the wall. The results indicate that bubbles follow a roughly zigzag path rather than maintaining a fixed distance from the walls. Bubble rotation and deformation, liquid turbulence, bubble-induced flow (such as bubbles wakes), bubble–bubble interactions, wall effects and boundary layer dynamics (such as the slight variations in a bubble’s initial position relative to the boundary layers can lead to divergent paths over time), etc., all lead to the unique trajectories of each bubble. However, some common characteristics can be observed: as the liquid flow rate increases, bubbles tend to move closer to the walls; in the diverging section, bubbles travel at a closer distance h_B,W compared to their motion in the throat, as shown in Fig. 16(b). The minimum distances generally occur within the first 0–8 mm of the diverging section, suggesting that bubbles suddenly swerve upon entering the diverging section, later shifting toward the opposite side to form zigzag trajectories.

Bubble velocities are shown in Fig. 16, with velocities along the x-axis represented by U_B,x and those along the y-axis represented by U_B,y. U_L,ave denotes the cross-sectional mean velocity of the liquid. Figure 17(a) shows that the bubbles undergo two significant processes in the diverging section. In process 1, illustrated by the gray region in Fig. 17(a), the bubbles experience rapid deceleration over short distances. Subsequently, in process 2, they undergo normal deceleration akin to the liquid’s deceleration, due to the gradually increasing cross-sectional area. The rapid deceleration of process 1 typically occurs at x ≈ 0.1, 0.0, and −2.0 mm and concludes at x ≈ 7.0, 4.2, and 2.2 mm for Q_L = 2.44, 4.54, and 6.04 m³/h, respectively. Previous studies by Kuo and Wallis²⁴ and van der Geld et al.²⁶ indicated bubbles moving near channel centerlines, thus missing this pronounced deceleration process. By comparing the sideward velocity U_B,y in Fig. 17(b) with bubble trajectories in Fig. 16, it becomes clear that sudden trajectory changes occur concurrently with dramatic deceleration, suggesting a strong association with the Coanda effect. A peak in bubble velocity appears earlier as Q_L increases, and if the liquid flow rate is large enough (Q_L > 4.54 m³/h), dramatic deceleration may occur even before the bubbles enter the diverging section.

In our previous studies,^8,9 sideward motion was often overlooked due to simplified analysis requirements. To investigate whether this motion contributes to or influences dramatic deceleration, we plotted the velocity magnitudes of the bubbles U_B = (U_B,x2 + U_B,y2)^0.5 as well as the ratios of U_B,y to U_B,x in Figs. 17(c) and 17(d), respectively. The velocity magnitudes U_B mirror U_B,x, indicating the reality of the two distinct deceleration processes. From a kinetic energy perspective, as shown in Fig. 17(d), sideward velocities U_B,y contribute about 15%–30% to bubble deceleration under the current flow conditions.

Figure 18 illustrates the dramatic deceleration process along the bubble trajectory, where the local accelerations U_B,x∂U_B,x/∂x are obtained from the averaged values in Fig. 17(a). As mentioned earlier, the larger the Q_L, the more significant the deceleration. Peaks in the acceleration magnitude occur at around x = 4.0, 2.0, and −1.0 mm for Q_L = 2.44, 4.54, and 6.04 m³/h, respectively. With an increase in Q_L, the peak appears earlier. Similarly, when the liquid flow rate exceeds ∼6.04 m³/h, the maximum deceleration occurs even before the bubble enters the diverging section. The maximum magnitudes of the accelerations reach about 400–4000 m/s², which is 40–400 times the acceleration due to gravity, for Q_L ranging from 2.44 to 6.04 m³/h. This underscores the need for a deeper understanding of the intricate interactions between the bubble and the liquid.

The intense deceleration and acceleration observed in bubble motion before breakup significantly impact bubble dynamics and breakup mechanisms. This pre-breakup phase marks a critical transition from mean flow dominance to influence by local turbulent structures. Zhao et al.⁵² demonstrated how rapid accelerations cause bubble deformation and rotation, altering surface energy and internal pressure, and increasing breakup susceptibility. Rapid deceleration can enhance local turbulence by creating small-scale eddies,²⁹ which may reduce the energy barrier for turbulent-induced breakup. This study focuses on the zigzag motion, which keeps bubbles in areas of high TKE, affecting the breakup frequency and size during downstream transport. In addition, sideward motion contributing 15%–30% of the bubble’s kinetic energy during deceleration indicates a non-uniform stress distribution. This section quantitatively analyzes these forces.

Figure 17 illustrates the relative velocity distribution between liquid and bubble in the x-direction (U_L,x − U_B,x), showing the dramatic deceleration processes experienced by the bubbles. The bubble surface is delineated by black outlines (void fraction α = 0.5), with blue regions in the contour indicating liquid velocities lower than the bubble’s and red regions denoting the opposite. The simulation results reveal a characteristic relative velocity distribution around the bubble, characterized by the upstream main flow exceeding bubble velocity and downstream recirculation flow lagging behind. This result provides clearer flow field information regarding the forces acting on the bubbles.

The zigzag motion reflects a complex interplay of hydrodynamic forces and fluid–structure interactions, governed by the Coanda effect, as illustrated in Fig. 20. In the diverging channel, velocity profiles skew due to viscous effects and adverse pressure gradients. Near the “convex” wall, boundary layer deceleration can lead to flow separation, forming a low-pressure region downstream. Main flow acceleration around the curve further reduces pressure according to Bernoulli’s principle. As bubbles enter this gradient, they experience a lateral force toward the wall, resulting from pressure gradients, lift forces due to near-wall velocity and rotation, and wall lubrication or repulsion forces. These forces determine the bubble’s position and trajectory. As bubbles approach the wall, acceleration between the bubble and wall reduces local pressure, reinforcing wall-ward motion, until opposing forces such as wall lubrication push the bubble away.

Pressure recovery downstream creates a more uniform pressure field, reducing the lateral force on the bubble. In the diverging section, the added mass force increases due to pressure recovery, resisting the bubble’s motion, while drag shifts from opposing to propelling the bubble as its speed drops below the surrounding fluid velocity.

In summary, the interaction of forces such as drag, lift, pressure gradients, added mass, and wall lubrication results in the zigzag trajectory. Understanding these forces allows for better prediction and control of bubble behavior in complex flows.

To elucidate the drag force evolution, Fig. 22 depicts the history of U_L,x − U_B,x distribution and bubble movement process from t = 0 ms to t = 1.2 ms, using Q_L = 4.54 m³/h as a representative case. In the throat section, the bubble typically encounters resistance due to negative U_L,x − U_B,x. Upon entering the divergent section, the bubble is simultaneously propelled by the mainstream and restrained by the downstream liquid. As the bubble undergoes dramatic deceleration and the relative velocity increases, the drag force transitions from resistance to propulsion, contributing to significant variable acceleration processes.

An intriguing phenomenon, highlighted by elliptical dashed lines in Figs. 22 and 23, is the presence of a liquid block at the front of the bubble moving faster than the bubble itself. This suggests that during the bubble’s deceleration, liquid with smaller acceleration is displaced, creating a higher pressure region at the bubble’s front. This pressure difference causes the surrounding liquid to push forward, creating a “block” that moves faster than the bubble and momentarily acts as a cushion against further deceleration. The faster-moving liquid implies that the bubble’s deceleration is influenced not only by pressure gradient and viscous drag but also by the added mass effect. The added mass force, representing the inertia of fluid moving with the bubble, is crucial here. The relative acceleration between the bubble and the liquid block is affected by the dynamic changes in the added mass coefficient, especially near walls or in complex flows.

The results of F_L–W are shown in Fig. 25 illustrating the interaction between the lift and wall lubrication forces. The total magnitudes of the lateral force approximately range from ±0.6 to ±3.1 × 10² N/m³ under various flow rates. The lift force significantly increases in regions where the Coanda effect occurs and decreases as the bubble moves away from the throat. This increase causes the bubble to approach the wall, inducing an opposing wall lubrication force. Consequently, the bubble is pushed away from the wall until equilibrium is reached between these forces. This lateral motion accounts for ∼15% to 30% of the bubble’s forward displacement [Fig. 17(d)] and directs the bubble into a low-velocity region, where the lift force substantially contributes to the dramatic deceleration of the bubble. To summary the Coanda effect, which amplifies the lift force, plays a significant role in the lateral movement, pushing the bubble toward the wall. However, the increasing wall lubrication force provides resistance, preventing the bubble from adhering to the wall and maintaining a balance that ultimately influences the bubble’s trajectory and deceleration.

Figures 22 and 23 highlight the non-uniform nature of the velocity and acceleration fields around the bubble. The differential acceleration observed indicates that the bubble is navigating through a stratified velocity field, where different layers of fluid exhibit varying levels of inertia. Such a stratified flow environment can lead to interesting dynamics, including the stretching and deformation of the bubble shape, as the forces acting on different parts of the bubble surface are not uniform.^38,39

Figure 26 illustrates the ratio of shear stress at the bubble–liquid interface to that at the bubble top-point, elucidating the complex stress distributions experienced by bubbles as they traverse the diverging section. Two distinct peaks in shear stress are observed: one acting on the bubble bottom (τ_bot,x) and the other acting on the bubble top (τ_top,x). For flow rates Q_L = 2.44 m³/h and Q_L = 4.54 m³/h, both τ_bot,x and τ_top,x are negative, but their effects are opposing. The negative τ_bot,x propels the bubble forward, while the negative τ_top,x restrains its motion. This opposing stress configuration results in the stretching of the bubble, leading to elongated deformation and ultimately binary breakup, as depicted in Fig. 14. Interestingly, at Q_L = 6.04 m³/h, a different stress regime emerges with τ_top,x > 0 and τ_bot,x < 0, indicating that the bubble experiences hindrance from both its top and bottom. This unique stress distribution may contribute to more complex deformation patterns and potentially different breakup mechanisms at higher flow rates. The observed shear stress distributions reveal a fundamental mechanism driving bubble deformation and shear-induced breakup in Venturi-type bubble generators. The non-uniform stress distribution arises from the bubble’s movement between a high-speed mainstream and low-speed recirculation flow near the wall. Future research could explore shear stress distributions and other forces such as pressure gradients and surface tension for a comprehensive understanding of bubble dynamics in complex flows.

This paper studied the dynamics of bubbles moving through a diverging channel near the wall using both visualization experiments and numerical simulations with a modified VOF model. The main conclusions are as follows:

1. Integrating the van der Geld et al. drag model through UDS methods marks a significant improvement in simulating bubble dynamics in accelerating flow fields. The modified approach has improved the accuracy of drag force predictions, yielding simulation results that match closely with experimental observations on bubble trajectory and velocity.

2. The deceleration of the bubble is closely tied to its zigzag trajectory, governed by the interplay of pressure gradient, drag (accounting for the acceleration flow effect), added mass (estimated from the bubble motion balance), and lateral forces (results of the lift against the wall lubrication). The largest magnitudes of these forces are −10⁶, 10⁶, −10⁶, and ±10² N/m³, respectively.

3. The zigzag trajectories and significant changes in rotational perturbation of single bubble wake are due to alterations in near-wall flow structures and wake dynamics, reflecting the instability in bubble transport influenced by non-uniform vortex structures and shear stresses.

4. Forces acting on the bubble have distinct roles in its rapid deceleration. The pressure gradient and added mass forces are the primary drivers of deceleration. Initially, drag force provides slight resistance but later acts as a propellant during deceleration, while the added mass force opposes this effect. Lift forces induce lateral motion, steering the bubble into slower-moving liquid and thus significantly contributing to its deceleration.

5. A novel simulation result reveals a faster-moving liquid region downstream of the bubble during rapid deceleration. This behavior underscores the critical role of added mass force in bubble dynamics, indicating that momentum from the surrounding fluid regions also influences deceleration.

The deceleration and deformation of the bubble contribute to the formation of fine bubbles in the diverging section, as prolonged residence time in the regions of strong shear and turbulence facilitates this process. Future studies should explore modified closure relationships for the forces in the bubble motion equation, particularly for systems containing multiple bubbles.

This work was supported by the National Natural Science Foundation of China (Grant No. U2333210) and the Scientific Research Foundation of Civil Aviation Flight University of China (Grant Nos. J2021-105 and J2022-002). The authors acknowledge the Department of Energy and Power Engineering, College of Water Resources and Hydropower, Sichuan University, for providing the experimental facilities and equipment. Special thanks are extended to the faculty and students for their support throughout the course of this research.

The authors declare that they have no competing interest.

Liang Zhao: Formal analysis (lead); Investigation (lead); Methodology (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead). Jiang Huang: Formal analysis (equal); Investigation (equal); Methodology (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal). Caozhi Chen: Formal analysis (equal); Writing – review & editing (equal). Jianan Gao: Formal analysis (equal); Writing – review & editing (equal).

No data were used for the research described in this article.