Issue 
Mechanics & Industry
Volume 25, 2024



Article Number  21  
Number of page(s)  16  
DOI  https://doi.org/10.1051/meca/2024016  
Published online  05 July 2024 
Original Article
Numerical simulation of single bubble motion fragmentation mechanism in Venturitype bubble generator
^{1}
North China Oil&Gas Company of SINOPEO, Zhengzhou 450007, China
^{2}
College of New Energy, China University of Petroleum (East China), Qingdao 266580, China
^{3}
Shandong Provincial Key Laboratory for Treatment of Oilfield Produced Water and Environmental Pollution, Dongying 257026, China
^{4}
Sinopec Petroleum Engineering Design Co., Ltd., Dongying 257026, China
^{*} email: liq@upc.edu.cn
Received:
8
December
2023
Accepted:
14
May
2024
Microbubbles have been widely used in power, chemical, mining and petroleum applications to improve energy efficiency. The venturetype bubble generator can improve reaction efficiency in chemical engineering as an efficient bubble generation method. Studying the flow field and bubble bursting mechanism can enhance the bubble generation performance. The volume of fluid (VOF) multiphase flow model was used in the Open Field Operation and Manipulation (OpenFOAM) framework to study the deformation and fragmentation behaviour of a single bubble in a Quasi threedimensional. The relationship between the mechanism of bubble breakup and the flow field in the diverging section of a Venturitype bubble generator was revealed. The reasons for the uneven distribution of bubble size were analyzed and discussed. The numerical simulation result shows that the fragmentation of a single bubble injected by the converging section is more substantial than that injected by the throat in the axial direction. Bubble fragmentation occurs mainly in the diverging section. The bubble's trajectory is highly similar to the vortex trajectory of the diverging section. The size of subbubbles generated by single bubble fragmentation decreases with the diverging angle and liquid flow rate increase. The differential distribution of turbulent kinetic energy in the radial position directly leads to the uneven distribution of the bubble size of the broken bubbles. As the liquid Reynolds number increases, static and dynamic erosion breakup are more prominent. The size of the subbubbles is much smaller.
Key words: Venturitype bubble generator / single bubble breakup / vortex movement / OpenFOAM
© J. Chen et al., Published by EDP Sciences 2024
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Abbreviations
OpenFOAM: Open Field Operation and Manipulation
CFD: Computational fluid dynamics
R^{2} : Curve fitting determination coefficient
c: Controllable compression factor
Fσ: Surface tension model momentum source term
σ: Surface tension coefficient
1 Introduction
Microbubbles have been widely used in power, chemical, mining and petroleum applications. In the petroleum industry, the high viscosity of heavy and residual oils results in incomplete combustion and a tremendous waste of energy. Therefore, researchers have come up with the idea of injecting microbubbles into heavy and residual oils so that a considerable number of bubbles can be formed, thus increasing the combustion area of the liquid fuel and improving combustion efficiency. In addition, microbubbles are also the basis for enhanced mass transfer by increasing the interface contact area and phase interface renewal frequency. Microbubbles have a large specific surface area, vital concomitant, enormous internal pressure, long residence time in water, etc., and are an excellent carrier to strengthen the mass transfer between gas and liquid in industrial production [1]. Efficient bubble generation technology can improve energy efficiency in crude oil refining, liquid fuel combustion, and other processes. It is also a hot topic in the research of petrochemical systems.
Up to now, the structure design and performance of the Venturitype bubble generator have been extensively studied. However, less research has been done on the bubble breakup process and mechanism of Venturitype bubble generators. The general view is that bubble transport in a Venturitype bubble generator is a complex unsteady flow phenomenon [2,3], including many complex flow characteristics, such as turbulence, twophase flow, and cavitation [4]. Early research on bubble breakup showed that bubble breakup in turbulent fields results from the collision of turbulent vortices of similar size to bubbles [5]. In experimental studies [6], bubble breakage in the diverging section caused by a pressure surge was experimentally captured. Zhao et al. [7,8] measured the bubble velocity and estimated the bubbling force. It was concluded that the rapid deceleration of the bubble in the diverging section had an essential effect on the process of bubble breaking. Uesawa et al. [9] experimentally analyzed the bubble fragmentation process and concluded that the pressure difference on the bubble surface causes the bubble to breakup. Fujiwara et al. [10] believed that the instability of the bubble interface caused bubble fragmentation due to the rapid increase of pressure in the asymptotic section. Song et al. [11] attribute bubble fragmentation to a combination of turbulent vortextobubble collisions, bubble surface instability, and viscous shear. Ding et al. [12] used a fluid volume model coupled with a level set method to splice virtual bubbles into a steady flow in a venturi and validated the numerical model through visualization experiments. It is found that the bubbles have binary crushing mode and multiple crushing mode. The axial pressure gradient force of the binary fragmentation mode mainly exists in the centre region. The radial shear force leads to multiple fragmentation modes in the sidewall region.
From the experimental point of view, highspeed photography technology has helped study the Venturitype bubble generator. However, there are still some problems with the validity of data processing and measurement of pulsating velocities and pulsating pressures in venturi turbulent fields. The current understanding of the shedding, breaking, motion, diverging, and regulation law of bubbles under complex working conditions is minimal, which further limits the scaledup design and industrial application of Venturitype bubble generators.
Numerical simulation has become an effective means of studying the complex multiphase flow characteristics in venturi channels with the development of computational fluid dynamics (CFD) technology. Ju et al. [13] concluded from numerical simulations that there is enormous shear stress at the entrance of the dilated end, and the bubbles show complete breakup. Ding et al. [14] calculated the bubblebreaking behaviour based on the VOF model in Fluent for different initial bubble positions and found that the turbulent dissipation rate is negatively correlated with the diameter of the bubble breakup. Zhao et al. [15] calculated the diverging section and the force state of the bubble. The results show that the rotational motion of the bubble is caused by shear stress, and pressure gradient force plays a crucial role in the deformation of the bubble. Li et al. [16] investigated the bubble fragmentation mechanism in the Venturitype bubble generator under a low flow rate based on numerical simulation and found that vortex plays a significant role in the transition of different fragmentation mechanisms. The proportion of the disturbed broken area expands with the movement of the vortex upstream in the diverging section, and the efficiency of bubble generation increases accordingly. Bie et al. [17] studied the dynamics and rupture of bubbles in a swirlventuri bubble generator (SVBG) using a highspeed camera system and a bubble image velocimeter. The research shows that the two trajectories of bubbles in the axial and radial directions are closely related to foam bursting. The vortex region expands for relatively high liquid Reynolds numbers, and the bubble size decreases.
In the current work, the breakage process of a single bubble in a Venturitype bubble generator was studied with the multiphase flow model VOF and the Large eddy simulation (LES) turbulence model in the OpenFOAM framework. The deformation and breakup behaviour of the initial bubble in a twodimensional flow field were studied by numerical simulation. The influence of initial bubble position, vortex motion, flow velocity, diverging angle and other factors on bubble fragmentation characteristics was investigated. The causes of uneven bubble size generated by the Venturitype bubble generator were discussed and analyzed. The above research provides a theoretical reference for improving the foaming performance of Venturitype bubble generators, producing microbubbles more efficiently and improving the energy utilization efficiency of the petrochemical industry.
2 Numerical methodology
2.1 Phase equation
The mainstream interface capturing methods mainly include the LS (Level Set) [18] method and the VOF [19] method. As an interface tracing method built under Eulerian mesh, the VOF method provides a simple and economical way to trace free boundaries in 2D or 3D meshes. It ensures fluid mass conservation during the computational process and handles the nonlinear phenomena of free surface reconstruction. It is more flexible and efficient than the finite difference method in dealing with complex free surface configuration problems.
The phase equation of multiphase flow is the continuity equation for each phase. In the VOF method, a variable α is introduced to represent the phase fraction of the fluid in the grid cell. The phase equation can be expressed as
$$\frac{\partial \alpha}{\partial t}+U\cdot \nabla \alpha =0.$$(1)
where U is the velocity of the mixed phase in the grid cell, the phase fraction α has distinct numerical characteristics in a gasliquid twophase flow system, where the liquid phase fraction. The α_{L} and the gas phase fraction α_{G} satisfy α_{L } + α_{G} =1. If the grid cell is filled with liquid phase α_{L} = 1, if the grid cell is filled with gas α_{G} = 1. If the value is between 0 and 1, the grid cells are mixed with gas and liquid.
The method used in OpenFOAM is to add an artificial compression term to the phase equation of interFoam to squeeze the phase fraction near the phase interface to counteract the blurring of the phase interface due to this numerical dissipation. The artificial compression term in the phase equation is solved based on the MULES (Multidimensional Universal Limiter with Explicit Solution) algorithm of the FCT (Flux Corrected Transport) method [20, 21]. Therefore, the phase equation (1) can be rewritten in the following form:
$$\frac{\partial \alpha}{\partial t}+\nabla \cdot (\alpha U)+\nabla \cdot \left(\alpha (1\alpha ){\text{U}}_{\text{c}}\right)=0.$$(2)
The U_{c} is the compression velocity to be modelled at the phase interface, and its direction is perpendicular to the interface. Specifying c as the controllable compression factor, U_{c} may be expressed as
$${\text{U}}_{c}=c\leftU\right\frac{\nabla \alpha}{\left\nabla \alpha \right}.$$(3)
From equation (2), it can be seen that as the value of c increases, the larger the value of U_{c}, the more pronounced the interface compression effect is. Numerical stability may decrease with increasing phase interface resolution as c increases. The physically unrealistic high speeds at the phase interface are parasitic currents [22]. Deshpande et al. analyzed and tested the problem and control of parasitic currents in interFoam. By testing, the generation of parasitic currents was found to be excessive or decaying with the coulomb number in the unsteady state, whereas almost no droplet displacement occurs in the steady state calculations [21]. This is also similar to Harvie's findings that the generation of parasitic currents is only secondary to the moving interface [22]. Considering the limitations of the size of the venturi model and the calculation conditions in this paper, the value of c is set to 1 [15]. The final form of the phase equation of the interFoam solver is as follows.
$$\frac{\partial \alpha}{\partial t}+\nabla \cdot (\alpha U)+\nabla \cdot \left(\alpha (1\alpha )c\leftU\right\frac{\nabla \alpha}{\left\nabla \alpha \right}\right)=0.$$(4)
2.2 Governing equations
For the gas–liquid twofluid model, the multinomial fluid continuity and momentum equations considering gravity and source terms can be expressed as
$$\frac{\partial \rho}{\partial t}+\nabla \cdot (\rho \text{U})=0,$$(5)
$$\frac{\partial \rho \text{U}}{\partial t}+\nabla \cdot (\rho \text{U}\text{U})\nabla \cdot \left(\mu (\nabla \text{U}+\nabla {\text{U}}^{T})\right)=\nabla p+\rho \text{g}+{\text{F}}_{\mathrm{\sigma}},$$(6)
where ρ is the mixed phase density, µ is the mixed phase viscosity, and F_{σ} is the surface tension model momentum source term. The continuous surface tension model CSF (Continuum Surface Force) method can convert the surface tension at the gas–liquid interface into an equivalent bulk force using an additional pressure gradient. F_{σ} is
$${F}_{\sigma}=\sigma \nabla \cdot \left(\frac{\nabla \alpha}{\left\nabla \alpha \right}\right)\nabla \alpha ,$$(7)
where σ is the surface tension coefficient, the final expression of the momentum equation in the interFoam solver takes the form.
$$\frac{\partial \rho U}{\partial t}+\nabla \cdot (\rho UU)\nabla \cdot \left(\mu (\nabla U+\nabla {U}^{\text{T}})\right)=\nabla {p}_{\text{rgh}}g\cdot \text{h}\nabla \rho +\sigma k\nabla \alpha .$$(8)
2.3 Turbulence model
The gas phase is transformed into a multiscale cluster of bubbles through a complex crushing process that occurs in only a few milliseconds in a venturetype bubble generator. It is crucial to select a turbulence model that can capture the largescale structural evolution and transient turbulence characteristics of the flow field. LES is applied in many related studies [23]. In this work, LES obtained the largescale turbulent structure by directly solving the Navier–Stokes equations, while the subgrid scale ones are modelled using the Smagorinsky closure model [24]. Dhotre et al. [25] stated that the Smagorinsky model could obtain results similar to the higherlevel Smagorinsky model with reasonable model constants. Therefore, the constant CSGS model in the Smagorinsky model was selected to be 0.094, according to the literature [26].
3 Numerical procedures
3.1 Structure of Venturitype bubble generator
The venturetype bubble generator used in this study is a rectangular structure, which is a quasithreedimensional channel. Figure 1 shows the schematic and geometry of the venturetype bubble generator, divided into a converging section, throat, and diverging section with different lengths and angles. The diverging section has sufficient width to improve the liquid phase shear force and promote the shear breakup process of bubbles. Considering practicality, the angles of the converging section and the diverging section of the bubble generator are 45° and 15°, respectively.
Fig. 1 Schematic drawing of the Venturitype bubble generator. 
3.2 Grid independence and comparison with experimental results
The grid is structured and divided into three parts ABC, as shown in Figure 2. A previous study found that the interaction between the liquid and gas is most pronounced at the entrance location of the evanescent region [27]. Hence, it is possible to densify the grid in zone B. The grid of zone A used a rougher grid. Based on the above requirements, four grids were studied independently, i.e., with 80,000, 100,000, 120,000, and 130,000 cell grids, respectively. In this case, the inlet flow rate was set to 0.162 m^{3}/h, and the bubble generation was in the middle of the converging section with an initial bubble size of 1 mm. The timeaveraged velocity distribution at x = 8 mm is obtained as the grid size increases to 130,000, as shown in Figure 3. The difference in timeaveraged velocity between the cell sizes of 120,000 and 130,000 is minimal. Therefore, a grid size of 120,000 is a relatively economical solution, satisfying y+ ≈ 1 for the throat and adjacent zones.
Fig. 2 Meshing schematic of the calculation domain. 
Fig. 3 Results for the mesh independence study timeaveraged velocity (Umean magnitude) distribution at x = 8 mm. 
3.3 Boundary conditions and bubble patching method
The fluids for the numerical simulations are water and air. Due to the low velocity of bubbles in the venturi channel (v < 100 m/s) and negligible gas density. Both water and air are set as incompressible fluids. The physical parameters of air and water are shown in Table 1.
The motion behaviour of a bubble in a venturetype bubble generator under different flow conditions was investigated by changing the inlet velocity. The fluid inlet was set to velocity inlet, and the outlet was set to pressure outlet. The numerical simulation working conditions are shown in Table 2, where Re represents the inlet Reynolds number, and Ca represents the capillary number.
Pressure−velocity coupling was solved by the semiimplicit pressurelinked equation (SIMPLE) via a relationship between velocity and pressure corrections. The PIMPLE algorithm, which uses a combination of PISO and SIMPLE (SemiImplicit Method for PressureLinked Equation), is able to handle transient calculations better under significant timestep conditions. The basic idea of the PIMPLE algorithm is to consider the iterative calculations within each time step as a steady state and utilize the SIMPLE algorithm to calculate the steadystate flow within a single time step. The standard PISO algorithm is used to compute the last step and finalize the iteration when the steadystate computation reaches a predetermined number of iterations [15]. The calculation step was adjusted to 2 × 10^{−6} s. Every 1000 steps of calculation were automatically saved. The gravity direction was set to the negative Y direction.
During the numerical simulation, gravity and the presence of residual curves closely monitor the stability of the flow field. When the flow field is stable, we suspend the simulation to add a spherical bubble to the flow field. The simulation is continued afterwards. The diameter of the patched bubble is ϕ. This is done as follows: (1) a spherical region with a specific location and size is marked for adaptation through the region; (2) the gas phase is patched to initialize the specific region during the initialization period, and the α_{G} is set to 1.
Main parameters for simulation calculation.
Numerical simulation of working conditions.
3.4 Experimental validation
Figure 4 shows a macroscopic comparison between the numerical results in the present study and previous experimental results of Mo et al. [28] based on the same venturi channel. The results demonstrate the change in velocity magnitude along the flow direction of a typical bubble as it enters the throat. The experimental and numerical simulation results show a consistent trend. In the throat section, the bubble velocity remained stable along the flow direction and gradually slowed down at the entrance of the diverging section. The bubble experienced a sharp deceleration due to the rapid splitting and bursting, which indicated that the bubble was experiencing tremendous resistance. The decrease in flow velocity led to increased static pressure in the diverging section along the flow direction, preventing the bubble from moving downstream and slowing down the bubble along the flow direction. In summary, the simulation results are consistent with the experimental results.
Fig. 4 Experimental and numerical simulation results of the velocity of bubble motion along the flow direction. 
4 Results and discussion
4.1 Effect of the vortex on single bubble motion fragmentation
To investigate the effect of vortex motion on the fragmentation characteristics of individual bubbles under operating conditions with Re = 27,950. After the flow field was stable, a single bubble with a diameter of ϕ1 mm was injected into the patch at the position of space coordinate (0.02, 0, 0), which meant the single bubble was injected from the water inlet. The transient diagram and streamline diagram postprocessing results of bubble movement, deformation and fragmentation are shown in Figure 5. The complete transient diagram of the bubble movement and breakup between the two red dotted lines. All the transient diagrams are combined to form the evolutionary process of the bubble in vortex generation, breaking, collapse, and block coalescence. Through observation, it can be found that there are three prominent areas in the process: steadystate disturbed breakup, unsteadystate disturbed breakup and concomitant vortex disturbance. Vorticity (Ω) is usually used to measure the intensity and direction of the vortex, which is defined as:
$$\text{\Omega}=\frac{\partial vx}{\partial z}\frac{\partial vz}{\partial x},$$(9)
where vx and vz are radial velocity and axial velocity, respectively.
The vorticity variation curve on the axial centerline of the Venturitype bubble generator is shown in Figure 7. The fluid has almost no vortex intensity before entering the throat and produces fluctuations for the first time after entering the throat. Vorticity decreases rapidly and remains stable after entering the steady state disturbed breakup region. The vorticity increases sharply and fluctuates significantly when it enters the unsteady state disturbed break area, begins to fall back, and continues to decrease until it reaches the region concomitant vortex disturbance. The vorticity increases and fluctuates sharply once it enters the unsteady state disturbed breakup region, and it begins to descend until it reaches the vortex disturbance region. Although Li et al. [15] mentioned the mechanism of bubble fragmentation in the first two regions, the influence of vorticity on bubble disturbance fragmentation was not analyzed.
The bubble undergoes four processes: microdeformation, stretching and deformation, breakup, collapse and chunking, as shown in Figure 5. This gasphase transport process is further divided into 18 frames with a time interval of 0.6 ms in Figure 6. The gas–liquid phase interface is extracted using the isovolumetric method, with iso = 0.75 set, and the flow field cloud is a transient velocity distribution. The image identified by the blue border in Figure 6 captures the motionbreaking process of the bubble in the converging section. The bubble moves smoothly initially and undergoes small deformations. The bubble begins to stretch and deform under the impact of the fluid, then gradually evolves into a bullet shape, and finally, the tail begins to concave. The tail of the bubble breaks through the head of the bubble and splits into two bubbles. The flow field is relatively stable until the bubble enters the throat. The image identified by the purple border in Figure 6 captures the motionbreaking process and flow field characteristics of the bubble in the throat. The bubble's velocity increases significantly after entering the throat. The bubbles pass quickly through the throat in two parts. The bubble does not break up as the fluid velocity increases significantly in the narrow interval of the throat. The image identified by the red border in Figure 6 captures the motionbreaking process of the bubbles in the diverging section. The bubbles undergo a sharp breakup and collapse and form more small bubbles after entering the diverging section. The flow velocity is high in the middle region of the diffusion section. Vortices are formed near the wall due to the sidewall effect. Bubbles break up under the shear of the vortex. Combined with the vorticity change curve, it can be seen that the phenomenon of bubble fragmentation and disintegration is most significant in the unsteady state disturbed breakup region. However, bubbles can coalesce in the vortex disturbance region since bubbles experience a long period of turbulent collision and energy dissipation from convergence to divergence. Bubbles are more likely to break and collapse due to the increased instability of the bubble surface. The bubble breakup mechanism is a typical steadystate disturbed breaking mode. In this process, the velocity difference and forward motion state of the captured bubble surface indicate that viscous shear plays a dominant role in the fragmentation of bubbles, which are split into two main parts under the action of viscous shear force. Therefore, the mechanism of the steadystate perturbation breakup mode is due to interfacial instability caused by KelvinHelmholtz instability perturbation [29,30].
The motion characteristics of the bubbles are incredibly complex, and the bubbles are rapidly deformed and broken under the surface tension in the process. The velocity contour shows that the bubble surface has uneven velocity distribution and an unsmooth shape, which indicates that the acceleration effect acts in the different directions of the bubble. This bubble interface fluctuation can be described as the instability of the density interface under average acceleration, namely RayleighTaylor instability [31]. It can be seen that the captured bubble gradually enters the region with larger vortex pulsation from the corresponding position of the streamline diagram. The bubble undergoes significant stretching deformation and moves forward with an inclination angle at the end of the unsteadystate disturbed breakup region.
It is important to note that the bubbles captured in the upper part are giant. Viscous shear still plays an essential role in bubble deformation because of a tremendous velocity gradient at the bubble surface. Meanwhile, the bubble edges and lower ends are constantly shedding tiny bubbles, a typical erosive breakup shear process. It can be found that RT instability induces bubble deformation, viscous shear enhances defects on the bubble surface, and turbulent pulsations in the flow field cause vorticity changes, which ultimately lead to bubble breakage. The vorticity in the vicinity of the bubble gradually increases as the bubble continues to move downstream. However, the vortex moves diagonally upwards at x = 0.095 m due to the turbulent kinetic energy of the liquid not being sufficient to maintain the vortex's wallhugging motion after a certain distance from the jet has occurred. The broken bubbles at the upper and lower parts gather to form two core bubbles under the effect of the vortex when the static pressure in the channel rises. Meanwhile, the volume of the single bubble decreases significantly in the third stage, so the corrosive bubble fragmentation disappears. The core bubbles mainly accompanied the vortex for rotational stretching deformation without further breakup.
To investigate the effect of the axial position of the initial bubble on the bubble breakup characteristics under operating conditions with Re = 27,950. A single bubble with a diameter of ϕ1mm was injected into the patch at the spatial coordinate (0.035, 0, 0) after the flow field was stable, which indicated that the throat injected the single bubble. The bubble split into the upper and lower parts of the throat, as shown in Figure 8. The bubble began to deform more extensively after entering the bifurcation section. In addition, the reduction of water flow velocity and rebound of static pressure caused the bubble to break gradually. Precisely, the bubble undergoes concave deformation and tensile deformation before entering the throat when the bubble is injected into the converging section. Whether in the converging section or the throat, the bubble appeared to have significant deformation and tended to split into upper and lower parts. However, the degree of fragmentation in the subsequent flow field of the bubble at the converging section is more potent than that of the bubble at the throat. Since the bubble surface instability is enhanced after a more extended turbulent collision and energy dissipation in the converging section, it is easier to break up after entering the diverging section.
The relationship between the trajectory of a single bubble and the trajectory of the centre of the vortex is further explored by analyzing the trajectory of the bubble and the vortex. The upper and lower vortex centre coordinates of the Venturitype bubble generator were extracted at different times. The centre coordinates of the bubble movement and the centre coordinates of the vortex in the upper and lower parts were compared and analyzed, as shown in Figure 9. It could be seen that the upper part of the vortex and the bubble have been maintained against the wall motion, and the motion trajectory could be linearly fitted as a straight line. In contrast, the lower part is detached, and the motion trajectory is fitted as a parabola. The four trajectory lines were fitted separately to obtain the results shown in Table 3. It can be seen that all the curve fit coefficients of determination (R^{2}) are close to 1. In the process, the motion trajectory of the bubble centre and the vortex centre are highly coincident with the deformation, breakup, collapse and aggregation of the bubble vortex.
Fig. 5 Single bubble deformation and fragmentation transient diagram and streamline diagram in the Venturitype bubble generator. 
Fig. 6 Single bubble transport process with 18 frames. 
Fig. 7 The curve of vorticity along the flow direction on the axial centerline. 
Fig. 8 Transient diagram of throat injection single bubble motion fragmentation. 
Fig. 9 The trajectory of bubble centre and vortex centre movement. 
Results of fitting the trajectory line of bubble and vortex centre motion.
4.2 Velocity and turbulent kinetic energy distributions
According to the work of Baylar [32], turbulent kinetic energy can be given as
$$k={\left(\frac{\epsilon l}{{C}_{\mu}^{3/4}}\right)}^{2/3},$$(10)
where k is the turbulent kinetic energy, ε is the turbulent dissipation rate, l is the turbulence length scale, C_{μ} is the turbulent constant, and d_{h} is the diameter of the throat.
Two crosssectional locations in each of the three regions in the diverging section were selected to extract their radial velocity and radial turbulence parameters for characterization, as marked by the blue dashed lines in Figure 5 (x =0.05 m, 0.06 m, 0.07 m, 0.085 m, 0.105 m, and 0.12m). The results of the curves are shown in Figures 10 and 11. At x = 0.05 m and 0.06 m positions, the maximum velocity and minimum turbulent kinetic energy are distributed in the central region (–0.004 m ≤ y < 0.004 m). The velocity in the central region does not change with the radial position, and the shear effect weakens. The fluctuations and turbulence between the layers are very weak. In contrast, the velocity gradient is more prominent, and the shear stress and turbulence are more substantial in the sidewall region, where the bubble fragmentation is more intense. It can be seen that the turbulence parameter in the middle of the diverging section is much lower than in the sidewall region, which results in a weaker fragmentation of the bubble in the middle of the diverging section. Meanwhile, the bubble size of the subbubbles produced in the intermediate region is more significant. The radial velocity and turbulent kinetic energy change gradually as the flow progresses into the unsteady disturbance breakup region. The radial velocity and turbulent kinetic energy centre deviate from the crosssection at x = 0.085 m, and the maximum radial velocity is still concentrated in the central region. However, the turbulent kinetic energy in the middle region is more prominent than in the side wall region. The radial velocity gradient and turbulent kinetic energy along the flow direction at the interface between the centre and the side wall region decrease gradually. The deviation trend is also maintained, and the interface gradient between radial velocity and turbulent kinetic energy is further reduced in the region vortex disturbance. The main reason is that the vortex cannot keep moving forward and tends to move upward after the jet reaches a certain distance, which leads to the deviation between the velocity centre and the distribution range of turbulent kinetic energy. The growth deviation of the vortex also causes the turbulent kinetic energy to increase gradually in the central region, resulting in the gradient of radial velocity and turbulent kinetic energy at the boundary between the side wall and the central position decreasing progressively. The growth deviation of the vortex also causes the turbulent kinetic energy to increase gradually in the central region, resulting in the gradient of the radial velocity and turbulent kinetic energy at the boundary between the side wall and the central position decreasing progressively.
Fig. 10 The velocity profiles are at different radial positions in different sections. 
Fig. 11 The turbulent kinetic energy is at different radial positions in different sections. 
4.3 Effect of liquid flow rate
Different patterns of bubble fragmentation are clearly observed in Figure 12. The liquid Reynolds number covers a range of 27,950–55,899. The acceleration, rotation and tensile deformation processes occur in three flow conditions after the bubbles enter the throat from the converging section for Re = 27950. The bubble enters the diverging section and begins to decelerate. Accordingly, the bubble rotation also experiences two similar stages: fast and slow. In both processes, the bubble is stretched and lengthened continuously, mainly due to the formation and growth of the vortex. The movement of the vortex causes the bubble to break, and the larger the flow velocity of the liquid, the more pronounced the bubble breaks and collapses in the diverging section.
As the liquid Reynolds number increases, new patterns of breakup appear in the venturetype bubble generator. As soon as the mother bubble enters the diverging section, the bubble undergoes violent deformation in the process of movement and breakup up into a mass of tiny bubbles rapidly. The larger the Reynolds number bubble, the smaller the bubble size of the subbubbles. This breakup mode can be named “static erosion breakup”, proposed by Wang [33]. Meanwhile, the increase in Reynolds number leads to the increase in turbulence intensity. Static erosion breakup and dynamic erosion breakup are more evident under the effect of turbulent fluctuation and collision [17]. Static erosion breakup and dynamic erosion breakup are surrounded by blue and red lines, respectively, as shown in Figure 12. The mother bubble is broken into upper and lower parts by fluid erosion in a static erosion breakup mode. The mother bubble moves in the radial direction and deforms severely in the dynamic erosive breakup mode. The parent bubble breaks into a cluster of tiny bubbles during rotation and stretching.
The fluid velocity and turbulent energy at the position of section (x = 0.085 m) are distributed along the radial direction, as shown in Figures 13 and 14. The radial velocity and turbulent kinetic energy at the crosssectional position (x = 0.085 m) become progressively more extensive overall. However, the distribution form does not change as the Reynolds number increases. The more significant the Reynolds number, the more drastic the change of turbulent kinetic energy. The middle region still forms a low turbulent region, similar to the radial velocity distribution. The number of subbubbles produced in the side wall region is more significant, and the bubble size is smaller with the flow rate increase, while the bubble size of subbubbles in the centre region is less affected. The difference in velocity between the sidewall region and the centre region leads to the formation of vortexes in the sidewall region. This is the reason for the intense turbulence in the sidewall region. The interaction between vortexes and the main jet flow generates a stagnation effect and forms stagnation zones. Bubbles have a higher probability of experiencing vortices and stagnation zones in steadystate disturbed breakup regions. The vortex near the wall has a strong interaction with the bubble, increasing the deformation rate of the bubble. The bubble breakup is more favoured in the motion state. The stagnation zone prevents the radial movement of the bubble and causes the bubble to be sheared by the strong flow, resulting in a static erosive breakup. When the bubbles are not affected by the stagnation zone, the inward driving force in the central lowpressure region pushes the bubbles to move in the radial direction, which is the main reason for the formation of dynamic erosion breakup.
Fig. 12 Transient diagram of bubble motion fragmentation under different liquid flow rates. 
Fig. 13 Radial velocity distribution at x = 0.085 m section for different flow conditions. 
Fig. 14 Radial turbulence energy distribution at x = 0.085 m crosssection for different flow conditions. 
4.4 Effect of diverging section tension angle
Four models with diverging angles (β) of 15°, 20°, 25° and 30° were constructed to compare the effect of different diverging angles on single bubble fragmentation. The comparison results are shown in Figure 15. The bubble fragmentation becomes more apparent, and the bubble size of subbubbles generated by bubble fragmentation decreases gradually in the diverging section with the gradual increase of diverging angle. The side wall effect weakens the bubble fragmentation, resulting in the decrease of bubble size distribution uniformity of bubbles. The tiny bubbles moving with the core bubble are no longer visible, and the accompanying subbubbles have already collapsed along the flow direction in the region vortex disturbance when the divergence angle is 25° and 30°.
The fluid velocity and turbulence parameters along the radial distribution in the middle diverging section (x = 0.085 m) with different diverging angles are shown in Figures 16 and 17. As can be seen from the figure, the radial velocity at the section position shows a decreasing trend with the increase of diverging angle. The radial velocity is no longer unimodal but a fluctuating trend from the side wall to the central area. The greater the angle of the diverging section, the more serious the loss of kinetic energy and the faster hydrostatic pressure rises in the direction of flow. The distribution of radial turbulent kinetic energy parameters is similar to that of radial velocity. The fluctuations become more and more evident with the increase of diverging angles. The greater the turbulent kinetic energy is, the more significant the gradient of change from the edge wall to the central region is. On the contrary, the more pronounced the fluctuation is. The centre area of turbulent kinetic energy also tends to move to the lower part with the increase of diverging angle. The centre of the bubble tends more towards the lower part as the angle of divergence increases, as seen in Figure 14. Moreover, the dip angle of the bubble moving upward in the region concomitant vortex disturbance is more significant. In this condition, the bubbles are more susceptible to the vortex.
Fig. 15 Transient diagram of single bubble motion fragmentation in Venturitype bubble generator with different diverging angles. 
Fig. 16 Radial velocity at position x = 0.085 m for different diverging angle conditions. 
Fig. 17 Radial turbulence energy at position x = 0.085 m for different diverging angle conditions. 
5 Conclusions
In this paper, the VOF model and the LES large vortex simulation were used to simulate the fragmentation behaviour of a single bubble in a venturitube bubble generator. The instantaneous motion state of the bubble and the factors influencing bubble motion and fragmentation were analyzed. Based on this research, the following conclusions are obtained:
Under the given working conditions, there are three fragmentation mechanisms and motion states of a single bubble: the steadystate disturbed breakup, the unsteadystate disturbed breakup, and the concomitant vortex disturbance region. Bubble motion tracks are different in all three regions. The fitted data revealed that the bubble trajectory was highly similar to that of the vortex centre. Vortex motion is the critical factor for bubble motion and fragmentation.
The turbulent kinetic energy in the central region is smaller than in the sidewall region in the steadystate disturbed breakup region. The difference in velocity between the sidewall region and the centre region leads to the formation of vortexes in the sidewall region. This is the reason for the intense turbulence in the sidewall region. The interaction between vortexes and the main jet flow generates a stagnation effect and forms stagnation zones. The stagnation zone prevents the radial movement of the bubble and causes the bubble to be sheared by the strong flow, resulting in a static erosive breakup. When the bubbles are not affected by the stagnation zone, the inward driving force in the central lowpressure region pushes the bubbles to move in the radial direction, which is the main reason for the formation of dynamic erosion breakup. Differences in the radial position distribution of turbulent kinetic energy directly lead to uneven bubble size distribution. The fragmentation of the bubble injected by the converging section is more dramatic than that injected by the throat in the same axial direction.
As the liquid Reynolds number increases, static erosion breakup and dynamic erosion breakup are more pronounced. The size of the subbubbles formed by the mother bubble is much smaller. The central region of the fluid's turbulent kinetic energy tends to move lower in the evanescent section. The more intense the bubble breakup, the more subbubbles are produced in the sidewall region. The boundary layer appears to play the role of separation to strengthen the fragmentation and refinement of bubbles with the increase of diverging angles. The bubble size of the subbubbles also decreases gradually. However, the side wall effect weakens bubble fragmentation, reducing the uniformity of bubble size distribution. Subbubbles are more likely to break in the vortex disturbance region with a diverging angle between 25° and 30°.
The feasibility of a highly efficient bubble generator is demonstrated at a theoretical level by investigating the movement of single bubbles within a Venturitype bubble generator. An academic reference is provided for acquiring highly efficient microbubbles in the petrochemical industry.
Acknowledgements
Acknowledgement of all individuals involved in project of Open Fund of Shandong Key Laboratory for Oilfield Produced Water Treatment and Environmental Pollution Control.
Funding
This research was funded by [Open Fund of Shandong Key Laboratory for Oilfield Produced Water Treatment and Environmental Pollution Control] grant number [No. ZC06070007].
Conflicts of interest
The authors have nothing to disclose.
Data availability statement
The manuscript has relevant research data which can be disclosed upon request.
Author contribution statement
Conceptualization, Junliang Chen and Mao Lei; Methodology, Junliang Chen; Software, Mao Lei; Validation, Junliang Chen, Mao Lei and Shaobo Lu; Formal Analysis, Xiaolong Xiao; Investigation, Mingxiu Yao; Resources, Junliang Chen; Data Curation, Junliang Chen and Mao Lei; Writing − Original Draft Preparation, Junliang Chen; Writing − Review & Editing, Junliang Chen; Visualization, Shaobo Lu; Supervision, Xiaolong Xiao and Qiang Li; Project Administration, Mingxiu Yao and Qiang Li; Funding Acquisition, Qiang Li.
References
 A. Gordiychuk, M. Svanera, S. Benini, P. Poesio, Size distribution and Sauter mean diameter of micro bubbles for a Venturi type bubble generator, Exp. Therm. Fluid Sci. 70, 51–60 (2016) [CrossRef] [Google Scholar]
 A. Mosavi, S. Shamshirband, E. Salwana, K.W. Chau, J.H.M. Tah, Prediction of multiinputs bubble column reactor using a novel hybrid model of computational fluid dynamics and machine learning, Eng. Appl. Comp. Fluid 13, 482–492 (2019) [Google Scholar]
 S. Shamshirband, M. Babanezhad, A. Mosavi, N. Nabipour, E. Hajnal, L. Nadai, K.W. Chau, Prediction of flow characteristics in the bubble column reactor by the artificial pheromonebased communication of biological ants, Eng. Appl. Comp. Fluid. 14, 367–378 (2020) [Google Scholar]
 L. Fang, W. Li, Q. Li, Z. Wang, Numerical investigation of the cavity shedding mechanism in a Venturi reactor, Int. J. Heat Mass Transf. 156, 119835 (2020) [CrossRef] [Google Scholar]
 V.M. Tikhomirov, On the breakage of drops in a turbulent flow, in Selected Works of Kolmogorov A.N. (1991) [Google Scholar]
 Y. Nomura, S.I. Uesawa, A. Kaneko, Y. Abe, Study on bubble breakup mechanism in a venturi tube, in: ASMEJSMEKSME 2011 Joint Fluids Engineering Conference ASME (2011), pp. 2533–2540 [CrossRef] [Google Scholar]
 L. Zhao, Z. Mo, L. Sun, G. Xie, H. Liu, M. Du, J. Tang, A visualized study of the motion of individual bubbles in a Venturitype bubble generator, Prog. Nucl. Energ. 97, 74–89 (2017) [CrossRef] [Google Scholar]
 L. Zhao, L. Sun, Z. Mo, J. Tang, L. Hu, J. Bao, An investigation on bubble motion in liquid flowing through a rectangular Venturi channel, Exp. Therm. Fluid Sci. 97, 48–58 (2018) [CrossRef] [Google Scholar]
 S.I. Uesawa, A. Kaneko, Y. Nomura, Y. Abe, Study on bubble breakup behavior in a Venturi tube, Multiphase Sci. Technol. 24, 257–277 (2012) [CrossRef] [Google Scholar]
 A. Fujiwara, K. Okamoto, K. Hashiguchi, J. Peixinho, S. Takagi, Y. Matsumoto, Bubble breakup phenomena in a venturi tube, in 2007 Proceedings of the 5th Joint ASME/JSME Fluids Engineering Summer Conference ASME (2007), pp. 553–560 [CrossRef] [Google Scholar]
 Y. Song, D. Wang, J. Yin, J. Li, K. Cai, Y. Qian, W. Liu, H. Li, Study of bubble breakup mechanism in a venturi bubble generator applied in TMSR, in 17th International Topical Meeting on Nuclear Reactor Thermal Hydraulics ACM (2017) [Google Scholar]
 G.D. Ding, J.Q. Chen, Z.L. Li, X.L. Cai, Numerical simulation on the motion and breakup characteristics of a single bubble in a Venturi channel, Ind. Eng. Chem. Res. 60, 14613–14624 (2021) [CrossRef] [Google Scholar]
 X. Ju, L. Sun, W. Tang, H. Yun, C. Yan, Analysis of the operating characteristics of a Venturitype bubble generator for MSR, Nucl. Technol. 37, 120605–1–120605–6 (2014) [Google Scholar]
 G. Ding, J. Chen, C. Wang, C. Shang, M. Liu, X. Cai, Y. Ji, Structural design and numerical simulation of axialswirling type microbubble generator, Guocheng Gongcheng Xuebao/Chin. J. Process Eng. 18, 934–941 (2018) [Google Scholar]
 L. Zhao, L. Sun, J. Tang, L. Hu, Z. Mo, H. Liu, G. Xie, Investigation on bubble motion in liquid flowing through a rectangular Venturi channel, in 17th International Topical Meeting on Nuclear Reactor Thermal Hydraulics, ACM (2017) [Google Scholar]
 Q. Li, D. Ming, M. Lei, X. Guo, J. Liu, H. Zhu, L. Fang, Z. Wang, Numerical investigation on the coupled mechanisms of bubble breakup in a Venturitype bubble generator, Eng. Appl. Comp. Fluid. 16, 229–247 (2022) [Google Scholar]
 H.Y. Bie, Y.X. Li, L.C. Xue, Y. Wang, G. Liu, Z.R. Hao, W.Z. An, A visualized investigation of bubble breakup in a swirlVenturi bubble generator, Aiche J. 69, 14 (2023) [Google Scholar]
 S.J. Osher, R.P. Fedkiw, Level Set Methods and Dynamic Implicit Surfaces, SpringerVerlag (2022) [Google Scholar]
 H. Jasak, Error analysis and estimation for the finite volume method with applications to fluid flows, Imperial College London 8, A385 (1996) [Google Scholar]
 C.W. Hirt, B.D. Nichols, Volume of fluid (VOF) method for the dynamics of free boundaries, J. Comput. Phys. 39, 201–225 (1981) [Google Scholar]
 S.S. Deshpande, L. Anumolu, M.F. Trujillo, Evaluating the performance of the twophase flow solver interFoam, Comput. Sci. Discov. 5, 014016 (2013) [Google Scholar]
 D.J.E. Harvie, M.R. Davidson, M. Rudman, An analysis of parasitic current generation in volume of fluid simulations, Appl. Math. Model. 30, 1056–1066 (2006) [CrossRef] [Google Scholar]
 L.L. Wang, Large eddy simulation theory and its application, J. Hohai Univ. Nat. Sci. 32, 261–265 (2004) [Google Scholar]
 J. Smagorinsky, General circulation experiments with the primitive equations, J. Monthly Weather Rev. 99–164 (1963) [CrossRef] [Google Scholar]
 M.T. Dhotre, B. Niceno, B.L. Smith, Large eddy simulation of a bubble column using dynamic subgrid scale model, Chem. Eng. J. 136, 337–348 (2008) [CrossRef] [Google Scholar]
 D.K. Liliy, Proceedings of the IBM scientific computing symposium on environmental Sciencer, IBM Form no 195, 320–1951 (1967) [Google Scholar]
 J. Huang, L. Sun, M. Du, Z. Liang, Z. Mo, J. Tang, G. Xie, An investigation on the performance of a microscale Venturi bubble generator, Chem. Eng. J. 386, 120980 (2020) [Google Scholar]
 Z.Y. Mo, M. Du, Z.Y. Shao, G. Xie, H.T. Liu, L.C. Sun, Investigation of bubble transportation in Venturitype bubble generator, At. Energy Sci. Technol. 50, 1034–1039 (2016) [Google Scholar]
 Helmholtz, XLIII. On discontinuous movements of fluids, Philos. Mag. J. Sci. London, Edinburgh Dublin 36, 337–346 (1868) [CrossRef] [Google Scholar]
 B. Thomson, W. Kelvin, Hydrokinetic solutions and observations, in Baltimore Lectures on Molecular Dynamics and the Wave Theory of Light (2012) [Google Scholar]
 Rayleigh, Investigation of the character of the equilibrium of an incompressible heavy fluid of variable density, Proc. Lond. Math. Soc. 114, 170–177 ( 1882) [CrossRef] [Google Scholar]
 A. Baylar, M.C. Aydin, M. Unsal, F. Ozkan, Numerical modeling of venturi flows for determining air injection rates using fluent V6. 2, Math. Comput. Appl. 14, 97–108 (2009) [MathSciNet] [Google Scholar]
 X.Y. Wang, Y. Shuai, H.M. Zhang, Bubble breakup in a swirlventuri microbubble generator, Chem. Eng. J. 403, 126397 (2021) [CrossRef] [Google Scholar]
Cite this article as: J. Chen, M. Lei, S. Lu, X. Xiao, M. Yao, Q. Li, Numerical simulation of single bubble motion fragmentation mechanism in Venturitype bubble generator, Mechanics & Industry 25, 21 (2024)
All Tables
All Figures
Fig. 1 Schematic drawing of the Venturitype bubble generator. 

In the text 
Fig. 2 Meshing schematic of the calculation domain. 

In the text 
Fig. 3 Results for the mesh independence study timeaveraged velocity (Umean magnitude) distribution at x = 8 mm. 

In the text 
Fig. 4 Experimental and numerical simulation results of the velocity of bubble motion along the flow direction. 

In the text 
Fig. 5 Single bubble deformation and fragmentation transient diagram and streamline diagram in the Venturitype bubble generator. 

In the text 
Fig. 6 Single bubble transport process with 18 frames. 

In the text 
Fig. 7 The curve of vorticity along the flow direction on the axial centerline. 

In the text 
Fig. 8 Transient diagram of throat injection single bubble motion fragmentation. 

In the text 
Fig. 9 The trajectory of bubble centre and vortex centre movement. 

In the text 
Fig. 10 The velocity profiles are at different radial positions in different sections. 

In the text 
Fig. 11 The turbulent kinetic energy is at different radial positions in different sections. 

In the text 
Fig. 12 Transient diagram of bubble motion fragmentation under different liquid flow rates. 

In the text 
Fig. 13 Radial velocity distribution at x = 0.085 m section for different flow conditions. 

In the text 
Fig. 14 Radial turbulence energy distribution at x = 0.085 m crosssection for different flow conditions. 

In the text 
Fig. 15 Transient diagram of single bubble motion fragmentation in Venturitype bubble generator with different diverging angles. 

In the text 
Fig. 16 Radial velocity at position x = 0.085 m for different diverging angle conditions. 

In the text 
Fig. 17 Radial turbulence energy at position x = 0.085 m for different diverging angle conditions. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.