© J. Chen et al., Published by EDP Sciences 2024

Abbreviations

VOF: Volume of fluid

OpenFOAM: Open Field Operation and Manipulation

CFD: Computational fluid dynamics

LES: Large eddy simulation

R² : Curve fitting determination coefficient

ϕ: Bubble size

β: Diverging angle

U: Mixed phase velocity

α: Phase fraction

αL: Liquid phase fraction

αG: Gas phase fraction

c: Controllable compression factor

Fσ: Surface tension model momentum source term

σ: Surface tension coefficient

Uc: Compression velocity

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 Venturi-type bubble generator have been extensively studied. However, less research has been done on the bubble breakup process and mechanism of Venturi-type bubble generators. The general view is that bubble transport in a Venturi-type bubble generator is a complex unsteady flow phenomenon [2,3], including many complex flow characteristics, such as turbulence, two-phase 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 vortex-to-bubble 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, high-speed photography technology has helped study the Venturi-type 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 scaled-up design and industrial application of Venturi-type 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 bubble-breaking 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 Venturi-type 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 swirl-venturi bubble generator (SVBG) using a high-speed 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 Venturi-type 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 two-dimensional 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 Venturi-type bubble generator were discussed and analyzed. The above research provides a theoretical reference for improving the foaming performance of Venturi-type 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

(1)

where U is the velocity of the mixed phase in the grid cell, the phase fraction α has distinct numerical characteristics in a gas-liquid two-phase 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 (Multi-dimensional 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:

(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

(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.

(4)

2.2 Governing equations

For the gas–liquid two-fluid model, the multinomial fluid continuity and momentum equations considering gravity and source terms can be expressed as

(5)

(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

(7)

where σ is the surface tension coefficient, the final expression of the momentum equation in the interFoam solver takes the form.

(8)

2.3 Turbulence model

The gas phase is transformed into a multi-scale cluster of bubbles through a complex crushing process that occurs in only a few milliseconds in a venture-type bubble generator. It is crucial to select a turbulence model that can capture the large-scale structural evolution and transient turbulence characteristics of the flow field. LES is applied in many related studies [23]. In this work, LES obtained the large-scale turbulent structure by directly solving the Navier–Stokes equations, while the sub-grid 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 higher-level 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 Venturi-type bubble generator

The venture-type bubble generator used in this study is a rectangular structure, which is a quasi-three-dimensional channel. Figure 1 shows the schematic and geometry of the venture-type 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 Venturi-type 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³/h, and the bubble generation was in the middle of the converging section with an initial bubble size of 1 mm. The time-averaged velocity distribution at x = 8 mm is obtained as the grid size increases to 130,000, as shown in Figure 3. The difference in time-averaged 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 time-averaged 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 venture-type 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 semi-implicit pressure-linked equation (SIMPLE) via a relationship between velocity and pressure corrections. The PIMPLE algorithm, which uses a combination of PISO and SIMPLE (Semi-Implicit Method for Pressure-Linked Equation), is able to handle transient calculations better under significant time-step 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 steady-state flow within a single time step. The standard PISO algorithm is used to compute the last step and finalize the iteration when the steady-state computation reaches a predetermined number of iterations [15]. The calculation step was adjusted to 2 × 10⁻⁶ 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.

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 post-processing 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: steady-state disturbed breakup, unsteady-state disturbed breakup and concomitant vortex disturbance. Vorticity (Ω) is usually used to measure the intensity and direction of the vortex, which is defined as:

(9)

where vx and vz are radial velocity and axial velocity, respectively.

The vorticity variation curve on the axial centerline of the Venturi-type 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: micro-deformation, stretching and deformation, breakup, collapse and chunking, as shown in Figure 5. This gas-phase 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 motion-breaking 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 motion-breaking 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 motion-breaking 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 steady-state 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 steady-state perturbation breakup mode is due to interfacial instability caused by Kelvin-Helmholtz 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 Rayleigh-Taylor 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 unsteady-state 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 R-T 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 wall-hugging 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 Venturi-type 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²) 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 Venturi-type 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.

4.2 Velocity and turbulent kinetic energy distributions

According to the work of Baylar [32], turbulent kinetic energy can be given as

(10)

(11)

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 cross-sectional 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 sub-bubbles 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 cross-section 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 venture-type 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 sub-bubbles. 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 cross-sectional 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 sub-bubbles 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 sub-bubbles 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 steady-state 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 low-pressure 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 cross-section 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 sub-bubbles 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 sub-bubbles 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 Venturi-type 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 venturi-tube 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 steady-state disturbed breakup, the unsteady-state 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 steady-state 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 low-pressure 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 sub-bubbles 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 sub-bubbles 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 sub-bubbles 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 Venturi-type 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. ZC0607-0007].

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

All Tables

All Figures

	Fig. 1 Schematic drawing of the Venturi-type bubble generator.
In the text

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

	Fig. 3 Results for the mesh independence study time-averaged 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 Venturi-type 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 cross-section for different flow conditions.
In the text

	Fig. 15 Transient diagram of single bubble motion fragmentation in Venturi-type 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