Numerical investigation on the coupled mechanisms of bubble breakup in a venturi-type bubble generator

ABSTRACT The bubble breakup mechanism in a venturi-type bubble generator is a hot research topic, and the bubble breakup mechanisms induced by the axial regional flow field characteristics has long been studied and debated. A coupled volume of fluid and large eddy simulation numerical simulation method was conducted within the OpenFOAM® framework to track the two-phase interfaces and transient turbulence characteristics of the flow field in a venturi-type bubble generator. A turbulence inlet database was established to ensure the reliability of the numerical solution and the economy of the calculation. The results indicated that the bubble breakup induced by the axial regional flow field characteristics of the divergent section can be divided into a steady-state disturbed breakup pattern and an unsteady-state disturbed breakup pattern macroscopically, the latter playing a dominant role in generating fine bubbles. Furthermore, a series of typical bubble breakup processes were captured via a Lagrangian approach that showed that interfacial instability played an important role at each stage of bubble transport: Kelvin–Helmholtz instability caused bubbles to detach from the gas film during the steady-state disturbed breakup pattern, and Rayleigh–Taylor instability induced the deformation of detached bubbles in the unsteady-state disturbed region. Additionally, viscous shear force strengthened the interfacial instability in the binary breakup process, and turbulent fluctuations and collisions contributed to the shattering of defective bubbles. Moreover, a shearing-off process tended to occur in regions with small vorticity fluctuations, therefore playing a small role in creating large clusters of bubbles.


Introduction
Multiphase flow is an important method of heat and mass transfer, and the formation of microbubbles is a typical physical phenomenon in multiphase flow. Owing to their characteristics of large specific surface area, strong adsorption, and long residence time in water (Takahashi et al., 2003;Zhang et al., 2016), microbubbles occur widely in engineering fields Zhang et al., 2016;Zhao et al., 2017). It is recorded that Bashforth and Adams (1883) began the study of microbubble formation as early as 1883 and, since then, a variety of bubble formation methods has been developed. As a typical microbubble generator, the venturi-type bubble generator, based on the hydrodynamic method, has attracted extensive attention in industry and research fields because of its simplicity and remarkably low energy consumption, e.g. process enhancement, water purification (Fujiwara et al., 2003;Yao et al., 2016), mineral flotation (Paramar & Majumder, 2016), molten salt reactors  and water oxygenation (Andinet et al., 2016). CONTACT Qiang Li liqiangsydx@163.com In-depth understanding of the breakup mechanism during bubble transport is of great importance for optimizing the design of venturi bubble generators. At present, the general view is that bubble transport in a venturi-type bubble generator is a complex unsteady flow phenomenon (Mosavi et al., 2019;Shamshirband et al., 2020), including many complex flow characteristics, such as turbulence, two-phase flow and even cavitation (Fang et al., 2020). Early theories (Kolmogorov, 1949) about bubble breakup have proposed that bubble breakup in a turbulent field is the result of collision between turbulent vortices of a size similar to that of the bubbles. Liao and Lucas (2009) summarized mechanisms of bubble breakup in turbulence and expressed the opinion that the breakup of a bubble is determined by the hydrodynamic conditions in the surrounding liquid and its own characteristics. In Liao and Lucas's review, breakup mechanisms can be classified into four main categories: (1) turbulent fluctuation and collision (Coulaloglou & Tavlarides, 1977;Hui & Wei, 2007); (2) viscous shear stress (Grace, 1982); (3) a shearing-off process (Fu & Ishii, 2003); and (4) interfacial instability (Wang et al. 2005). It is noted that the studies mentioned above only focused on dealing with conventional channels. However, bubbly flow characteristics in a venturi generator is much more complex than that in a conventional channel owing to its divergent section design. Thang and Davis (1979) conducted experiments on the flow of bubbles in venturi channels, and the results indicate a significant difference between bubble breakup in a venturi divergent section and in a flat channel. Based on research in turbulent dispersions, more researchers have analyzed the breakup mechanism in venturi divergent sections through visualization experiments. Fujiwara et al. (2003Fujiwara et al. ( , 2007 and Nomura et al. (2011) measured the wall pressure in a venturi divergent section by an experimental method on a water-air system, and proposed that the fragmentation of bubbles was caused by shock waves. Zhao et al. (2017Zhao et al. ( , 2018) used a highspeed camera to study the process of bubble breakup in a venturi-type bubble generator, and they found an obvious stagnation region in the venturi divergent section that impeded the movement of the bubbles and enhanced interaction with the water. Obviously, the bubbles in this region were broken by strong shear in a flow field with a huge velocity gradient. Uesawa et al. (2012) believed that it was the pressure difference on the surface of the bubbles that led to the fragmentation of bubbles in a venturi-type bubble generator. Based on the force analysis of individual bubbles, Song et al. (2019) found that the breakup of bubbles is the combined action of the destructive stress (turbulence stress, viscous stress) and the allowable stress. Additionally, researches on bubble breakup with different liquid Reynolds numbers (Yin et al., 2015) and different venturi structural parameters (Lee et al., 2019) have also been reported in the literature.
As a practical and effective experimental approach, high-speed photography is able to capture the liquid-gas two-phase interface morphology and other visible characteristics of venturi tubes. However, in terms of flow field information analysis, the measurement of fluctuating velocity and pressure in turbulence presents great challenges experimentally. Apart from experiments, numerical simulation has become an efficient tool for investigating the complex multiphase flow characteristics of venturi channels so as to gain a basic understanding of the bubble breakup mechanism. Numerous simulation works have been reported in the past two decades. In the early stages, Wang and Chen (2002) simulated bubbly flows numerically and found that the relative motion between the two phases can be predicted well using a twofluid model. Rzehak and Krepper (2013) used Population Balance Modelling (PBM) coupled with the Reynolds-Averaged Navier-Stokes (RANS) method to study the gas-liquid two-phase flow in venturi structural channels. Zhao et al. (2019) investigated the effect of divergence angle on bubble transport based on the RANS model. However, the above studies on bubble transport failed to take into account the phase interface directly, and the RANS method is not suitable for the microscopic mechanism of bubble breakup. Considering that high grid resolution is required to capture the sharp gas-liquid interfaces and appropriate modeling cost, a coupled Volume Of Fluid (VOF) (Hirt & Nichols, 1981) and Large Eddy Simulation (LES) numerical simulation method has been conducted and is reported in this paper.
In the current work, a solver named interFoam was adopted from the open-source library OpenFOAM ® 6.0, which is based on the finite volume method for solving incompressible Navier-Stokes equations and isothermal immiscible fluids using a VOF method. Firstly, a turbulence inlet database was established to ensure the reliability of the numerical solution and the economy of the calculation. Then, the regional characteristics of the flow field in the divergent section were analyzed to study the conversion of the breakup mechanisms in the process of bubble transport, and three factors affecting the conversion were discussed. Finally, a series of typical bubble breakup processes were captured with the Lagrangian approach, and different breakup mechanisms were evaluated during the process.

Phase equation
In the VOF method, a variable α is introduced to represent the phase fraction of fluid in the grid cell and the phase equation can be expressed as where U is the mixture velocity. As a special passive scalar, the phase fraction α has obvious numerical characteristics. For a vapor-liquid two-phase flow, α L and α G are strictly constrained by α L + α G = 1, and the subscripts L and G represent the liquid and vapor phases, respectively. Clearly, a value of α between zero and one in some grid cells will inevitably lead to an unclear phase interface. To make the interface sharp, an artificial compression is added to the VOF phase equation, which is solved based on the Multi-dimensional Universal Limiter with Explicit Solution (MULES) algorithm of the Flux Corrected Transport (FCT) method (Boris & Book, 1997). Therefore, the phase equation (Equation 1) is rewritten as where the third term on the left-hand side is the artificial compressible term and c represents the controllable compression factor. As the value of c increases, the interface compression effect becomes more obvious. Considering the limitations of the venturi model size and calculation conditions in this work, the value of c is set to be one.

Governing equations
The VOF method assumes that all phases share a set of continuity and momentum equations, so the continuity and momentum equations considering gravity and source terms are given as where ρ is the mixture density, μ is the viscosity and the last term in Equation (4) (i.e. F σ ) is the added momentum source term due to the addition of the surface tension model. The source term of surface tension is the surface integral of the liquid volume, which is very difficult to solve directly. Therefore, the Continuous Surface Force (CSF) model (Brackbill et al., 1992) is used to convert this integral into solving a continuous volume force acting on the gas-liquid mixing zone. F σ can be simplified as F σ = σ κ∇α and κ is given by κ = −∇ · (∇α/|∇α|), where σ is the surface tension coefficient. As for the pressure term, it can be rewritten as p rgh = p − ρg · h. The modified pressure term makes the definition of boundary conditions simpler, where h represents the position vector of the body center of the grid cell. Thus, the final form of the momentum equation can be expressed as

Velocity-pressure coupling and implementation
For a long time, the Pressure Implicit with the Splitting of Operators (PISO) algorithm has been widely used as a common discretization algorithm in transient simulation calculations. However, the linearization of the convection term in the governing equation causes a lag in the flux at each step of the PISO algorithm (Jasak, 1996). Therefore, the PISO algorithm must adopt a small time step when conducting transient simulation, which means that it will consume a lot of time in the whole computation procedure. In this paper, the PIMPLE algorithm (Pham & Choi, 2021), i.e. PISO coupled with the Semi-Implicit Method for Pressure-Linked Equations (SIM-PLE), is adopted, which can better deal with transient computation for large time steps. The discretization of the PIMPLE algorithm will be introduced briefly, and the solution process of the governing equations in this paper is shown in Figure 1. Firstly, the predicted velocity U r can be obtained by semi-discretization of Equation (5) with the expression where U 0 and P 0 rgh are the defined initial velocity and pressure, respectively. H(U) is the H operator of the semi-discretized momentum equation and A p is the coefficient matrix of velocity for each grid in the momentum equation. Next, calculating the divergence of Equation (6) and coupling the continuity equation (Equation 3), the pressure Poisson equation of P n rgh in PISO to be corrected at the next time iteration can be solved as where n represents the time step in PISO, and P n rgh is obtained (here n = 1), and U n to be updated is given by When the PISO algorithm converges, output U n and P n rgh to the next loop of the PIMPLE algorithm. The convergence criteria for the PIMPLE loop are based on a fixed number of iterations, which is set as two in this paper. When PIMPLE converges, the solution for a time step is complete. The time step Δt in the solution process is automatically controlled by the CFL number condition, which is the maximum allowable CFL number a numerical scheme can use. The maximum number CFL max was set to be 0.5 in this work, and the CFL number condition becomes where |v| is the modulus of the velocity in a grid cell and Δx is the length of the grid in the direction of velocity.

Turbulence model
As mentioned above, in a venturi-type bubble generator, the gas phase is transformed into a multi-scale bubble group through a complex crushing process that takes place in a matter of milliseconds. 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 (Wang, 2004;Zhang, 2012;Zhou, 2014). In this work, LES obtained the large-scale turbulent structure by directly solving the Navier-Stokes equations, while the sub-grid scale ones are modeled using the Smagorinsky closure model (Smagorinsky, 1963). Dhotre et al. (2008) stated that the Smagorinsky model could obtain results similar to the higher level Smagorinsky model with reasonable model constants. Therefore, the model constant C SGS in the Smagorinsky model is selected to be 0.094 in accordance with the literature (Liliy, 1967).

Structure of the venturi-type bubble generator
The venturi-type bubble generator used in this study is a venturi channel with a rectangular structure improved by our laboratory. Figure 2 shows the schematic drawing and geometrical dimensions of the venturi-type bubble generator. This venturi channel has a long throat of varying width that is designed to enhance the shear of the liquid phase against the gas phase and encourage more bubbles to break up at the throat for the first time. Considering the practicability, the diameter of the injection hole is 1 mm and the angles of convergence and divergence are 22.5°and 7.5°, respectively. Its physical properties are shown in Table 1. It needs to be emphasized that the venturi channel will also be a simple cavitation reactor owing to a sufficient drop of local pressure in the throat. In order to prevent cavitation, the liquid flow rate Q L is limited to 0.324 m 3 /h, corresponding to a liquid Reynolds number Re L at the inlet of 5918.

Inlet conditions used for LES
For an unsteady turbulence model such as LES, the flow behavior within the domain is largely determined by the behavior at the inlet. Tabor and Baba-Ahmadi (2010) generalized that inlet conditions for LES should be stochastically varying, 'look' like turbulence, and be easy to implement and to adjust to new inlet conditions, etc. However, an accurate and effective method for generating turbulent inlet conditions does not exist. Previous research in the literature can be classified into two categories: precursor simulation methods (Lund et al., 1998) and synthetic methods (Bechara et al., 2012;Jarrin et al., 2006;Kraichnan, 1970;Wu, 2017). The method given in the present work belongs in the synthetic method category, and it produces a spatiotemporal-variant field by summing a set of pseudo-random numbers and spatiotemporal-invariant mean fields as where x p is the patch value, x ref refers to the spatiotemporal-invariant patch scalar, n is the time level, s represents the fluctuation scale and R is a pseudo-random number. ζ is the fraction of new random component added in a previous time value and its default value is 0.1 in this work. c corresponds to an automatically calculated heuristic correction term introduced to compensate for the energy level losses due to ζ . Apparently, in order to implement the inlet conditions correctly, the coefficients in Equation (10) need to be constrained. Therefore, a section of reducer pipe connected to the venturi-type bubble generator is introduced to develop the initial turbulence, and the liquid flow in the reducer tube is calculated based on the simpleFoam solver provided in OpenFOAM, which is a steady-state turbulence solver based on RANS. The schematic is shown in Figure 3. When a statistically steady state is obtained, the average velocity at the outlet of the reducer pipe will be used as the initial reference patch value x ref 0 , the extremum and average values of the outlet velocity components will be recorded and stored in a turbulence database. Therefore, the components of fluctuation scale s in the x-, y-and z-directions can be expressed as [|u| max − |u|, |u| − |u| min ] max ,[|v| max − |v|, |v| − |v| min ] max and [|w| max − |w|, |w| − |w| min ] max . The databases of numerical simulations involved in this paper are shown in Table 2.

Grid independence and comparison with experimental results
The mesh of the venturi-type bubble generator was drawn with a hexahedral grid using the mesh generating software Gambit and the mesh was divided into four parts as shown in Figure 4. Previous research (Huang et al., 2020;Zhao et al., 2018) has found that the interaction between the liquid and gas phases was most obvious at the entrance position of the divergent region, so it was necessary to encrypt the grid of region C in Figure 4. In addition, in order to generate turbulence better, the grid of region A in Figure 4 was also coarsened. Based on the above requirements, a mesh independence study was carried out on four meshes, i.e. with 0.82, 1.04, 1.32 and 1.61 million cells, respectively. A liquid flow rate of 0.162 m3/h and a void fraction β (the ratio of the gas volumetric flow rate to the total volumetric flow rate) of 3% were adopted for the tested case. The sensitivity of the spatial discretization was examined in terms of the gas phase distribution and the velocity distribution. As shown by Figure 5(a), with the number of grids increasing to 1.32 million, a steady gas phase distribution can be observed and the bubbles' shape are clear enough for this study. After temporal averaging, the time-averaged velocity distribution at x = 4.5 mm (along the blue line in Figure 5(a)) could be obtained and the result is presented in Figure 5(b). It seems that, in either the near wall regions or the center regions, the difference in Umean magnitude between the cell numbers of 1.32 and 1.61 million is very small. Hence, a mesh size of 1.32 million is a relatively economic scheme and this also satisfies y + ≈ 1 for the surface of the throat section and its  adjacent regions. A typical case of the gas phase transport process running 250,000 time steps under 20 processors costs nearly 140 h in total. Figure 6 shows a macroscopic comparison between the numerical results in the present study and previous experimental results of Zhao et al. (2019) based on the same venturi channel. Gray processing for the gas phase is adopted in the numerical results and bubble shape is extracted from an isovolume with α L between 0.5 and 0.82. It is apparent from Figure 6 that the gas phase distribution and bubble shape are quite similar to the images captured by Zhao et al. In the direction of the flow, dramatic changes in bubble transport and morphology can be clearly displayed in the area shown. Specifically, under the three operating conditions, the gas phase in the entrance of the expansion section starts to come off the wall in the form of stretched bubbles. Entering the diverging section, the stretched bubbles begin to be broken up further under the strong interaction between the gas and liquid, and the coalescence of bubbles also takes place at the back of the diverging section. In summary, the simulation results are in good agreement with the experimental results and the numerical method proposed in this paper is reliable and applicable for simulating the venturi-type bubble generator. Figure 7 shows a typical time-dependent evolution of gas phase transport in the venturi channel and the velocity field of the central slice. As observed, the first strong interaction of the gas-liquid phase occurs at the throat, which has not been the focus of previous studies. In fact, this phenomenon will have a huge impact on bubble transport in the divergent section.

Patterns of bubble transport
In the venturi throat, after gas was injected from millimeter-scale injection holes in the vertical direction, it began to form a film near the wall and then moved along the axial direction, effectively disturbed the flow of the liquid as it did so. It can clearly be seen from the velocity streamline diagram at T = 5 ms in Figure 7 that there was a large velocity difference between the mainstream and the gas film. As the contact area increased, the mixture flow in the throat developed into a typical shear flow and the shear velocity perturbed the interface through viscous action, which eventually led to Kelvin-Helmholtz instability (Hermann 2009;Lord 2009). It can clearly be seen in Figure 7 that K-H waves existed in the throat and presented axisymmetric characteristics at T = 10 ms and T = 15 ms. As the mixture flow escaped the constraint of the throat walls, serious deformation of the gas film was caused by the K-H instability at the entrance of the divergent section. Subsequently, the gas film began to move along the side wall of the divergent section, and the movement of the gas phase was stagnant at 20.2 ms. In the meantime, the stretched bubbles formed at the throat began to break away continuously from the gas film and split into two or more daughter bubbles and flowed with the liquid downstream. At a more macroscopic level, the mainstream surface distortion increased gradually, and the detached bubbles demonstrated repeated coalescence and fragmentation in the direction of the flow. It can be considered that the transport of bubbles was completely developed.  During the gas phase moving process, the flow field in the divergent section showed obvious axial regional characteristics, which may be the reason for different breakup mechanisms of the detached bubbles and the stretched bubbles attached to the gas film. A more detailed flow field description is revealed in the red frame at 20.2 ms in Figure 7: the dynamical vortices downstream of the gas film stagnation region hindered the bubbles from moving forward, giving the bubble flow near the wall in the first half of the divergent section distinct visual characteristics of continuity and steady state. This phenomenon has also been observed in visualized studies by Huang et al. (2018), Lee et al. (2019) and Zhao et al. (2019). Moreover, the mainstream shows obvious K-H instability, which indicates that gas-liquid shear flow is the main flow pattern in the first half of the divergent section.
Three macroscopic characteristics of the flow field are discussed below, including the velocity components distribution, the vorticity distribution and the Reynolds stress distribution. Figure 8 provides velocity components and vorticity distributions along the central axis in the x-direction at 20.2 ms. The axial velocity component presents a fluctuating downward trend after entering the diverging section, while the pulsation evolves to be more obviously disordered behind the stagnation region. Additionally, the vorticity and the components of velocity in the other directions also evolve more intensively behind the stagnation region, producing here fine bubbles with minimized diameter as expected.
Another more intuitive manifestation of the turbulence characteristics is shown in Figure 9. For statistically steady turbulent flow, a prime-squared mean of the velocity field gives the Reynolds stress field as where U denotes the instantaneous velocity field and is decomposed into a time-averaged mean component (i.e. U) and a fluctuating component (i.e.U ). Interestingly, the Reynolds stress also has obvious axial regional characteristics in the divergent section.
Integrating the above phenomena, it is concluded that the existence of the vortex has a significant impact on the bubble flow transport and breakup mechanisms macroscopically. The position where the irregular pulsation of the macroscopic characteristics occurs is also the transition point of the change of the bubble breakup mechanism in the venturi channel, which is the position of vortex blocking of the downstream extension of the gas film. For convenience, the dynamical vortex is called the blocking-vortex in this work. Additionally, bubble breakup occurring in the region before this transition point is signified as a steady-state disturbed breakup pattern. After this point, the pattern is denoted as unsteadystate disturbed breakup. In this section, the two patterns are represented by pink and blue shading in Figures 8  and 9. In fact, the location of the blocking-vortex in the venturi channel has a certain controllability, which will be elaborated in the following sections.

Effect of the inlet disturbance on the gas film
The disturbance is a necessary condition for turbulence to occur. Regarding the case with a minimum liquid flow  rate of 0.162 m3/h in this work, the average velocity of the main flow stream at the throat is 12.6 m/s, corresponding to the average Reynolds number of 2.2 × 10 4 , which is far beyond the scope of a laminar flow. As mentioned above, for the venturi channel studied in this paper, the generation of turbulent inlet conditions depends on a synthetic method; thus, the influence of the disturbance on the blocking-vortex can be studied by artificially changing the modulus s * of the fluctuation scale. Here, the results for three fluctuation scales (s * , 0.5s * and 0.25s * ) are chosen as the sample to set the operating conditions for tests, and the tested case still has the operating conditions demonstrated in Section 4.1.
As is discussed in Section 4.1, the dynamical blockingvortex downstream of the stagnation region blocks the gas film extension; therefore, the length of the gas film can be used to predict the position of the vortex in a venturi tube. Figure 10 illustrates grayscale images of the gas film, which is completely developed by extracting from an iso-volume after temporal averaging, and the velocity streamline of the whole flow field is also shown after adjusting the transparency. It shows that, as the fluctuation scale modulus decreases, the gas film characteristic lengths increases. This means that the position of the blocking-vortex is moving downstream, and further leads to the postponement of the unsteady-state disturbed breakup pattern. In addition, considering the gas phase distribution diagram in Figure 10, the distribution of the gas phase with different fluctuation scales share similar interface characteristics. In more detail, this continuous gas film structure can be divided into the wallattached gas film and the mainstream-attached gas film, depending on the object to which the film is attached. The reasons for such structural differences can be seen from the velocity streamline diagram, where there are some non-negligible small vortices ahead of a blocking-vortex forcing the gas phase to leave the wall, resulting in the formation of a mainstream-attached gas film.
Here, the gas film length L 0 is defined, and L 1 , L 2 correspond to the length of the wall-attached gas film and  the mainstream-attached gas film, respectively. Figure 11 presents a more concise relationship between the above lengths. The results express that, under three inlet disturbances, there is little difference between the wallattached gas films, but the mainstream-attached gas film is the main reason for the difference in the length of the gas film. Given the discussions above, it is reasonable to consider that enhancing the inlet disturbance is an effective way to narrow the field of the steady-state disturbance.  Figure 12 demonstrates the temporal average distribution of the gas film under different operating conditions in the venturi channel. As presented in Figure 12(a), under two different liquid flow rates, the length of the gas film increases when increasing the gas volume fraction. This indicates that the gas volume fraction also has a significant effect on the location of the blocking-vortex. In the limited space of the divergent section, the steady-state disturbed breakup pattern occupies a greater region with an increase of gas volume fraction, which further macroscopically explains why the venturi bubble generator has difficulty generating fine bubbles under large gas volume fractions (such as β = 7%).

Effect of the operating conditions on the gas film
The distribution of the gas film with two different gas volume fractions is shown in Figure 12(b). At a fixed gas volume fraction, the length of the gas film is not greatly affected by the liquid flow rate. This suggests that, when changing the working environment of a specific venturitype bubble generator, in order to achieve the desired effect (such as the desired bubbles size), the flow rate of the main stream and gas can be predicted, which is helpful for engineering applications to some extent.
Of course, the geometrical size of the venturi channel also has a significant impact on the blocking-vortex (Zhao et al., 2019). Further studies that take these variables into account will need to be undertaken.

Bubble morphology and breakup mechanisms in the divergent section
In order to obtain more detailed dynamic characteristics of bubbles, a series of typical bubble breakup processes are captured via a Lagrangian approach in this section, which will be helpful for better understanding of the bubble breakup mechanism. In the case tested here, the working condition of the processes has a liquid rate of 0.275 m3/h and a gas volume fraction of 3%. Figure 13(a) shows the starting and ending points of captured bubbles and the trace line of the bubbles in a three-dimensional view. Figure 13(b) shows the vorticity distribution of the central section at the final moment in the two-dimensional view. Figure 14 further divides the bubble transport process into 28 frames with an interval of 0.2 ms, and the starting time is denoted as t 0 with a value of 0.4 s. Each frame has a side size of 100 pixels to ensure a uniform bubble size scale. The shape of bubbles is extracted from an iso-volume with α L = 0.66, and the surface of bubbles is colored according to their transient velocity.
The process of a captured bubble falling off from the gas film can clearly be seen in the red frame in Figure 14, which belongs to the steady-state disturbed breakup pattern covered in detail in Section 4.1. When the time is t 0 + 1.2 ms, the bubble breakup turns into an unsteadystate disturbed breakup pattern, corresponding to a horizontal position of the bubble of 0.038 m in Figure 13. This exactly corresponds to the position of a blockingvortex under the corresponding working conditions in Figure 12(a), further proving that it is feasible to reflect the position of blocking-vortex by the gas film length after temporal averaging. In the steady-state disturbed breakup pattern region, the velocity distribution on the surface of the ligamentous bubble and the sloping forward movement state indicate that the viscous shear force plays a dominant role in the detached bubbles' breakup. Hence, the mechanism of the bubble breakup can be attributed to the interface instability caused by K-H disturbance.
As the detached bubbles move forward, the bubbles first undergo a second division as shown in the purple frames in Figure 14. In fact, the characteristics of the bubble movement during this process are extremely complex. After the bubbles are detached from the gas film, they will quickly change from ligamentous to elliptical under the action of surface tension. In addition, Figure 14 shows the uneven velocity distribution on the bubble surface during the second division process, which means there are accelerations in different directions on the bubble. Therefore, such bubble interface fluctuations can be described  as instability of the density interface under the action of normal acceleration, namely Rayleigh-Taylor instability (Lewis, 1950;Rayleigh, 1882).
At t 0 + 2.2 ms, it can be seen from the corresponding position in Figure 13 that the bubble enters an area where the vorticity changes dramatically. After 0.2 ms, an obvious deformation occurs, and the bubble is elongated and moves forward at a fixed inclined angle. The bubble size is still too large, and there is still a large velocity gradient at the interface, so that the viscous shear force still plays an important role in the bubble deformation. In addition, it should be noted that there are constantly tiny bubbles falling off at the edge of the bubble, which is a typical shearing-off process. With the gradual disorder of vorticity near the bubble, the bubble finally breaks up (to be more precise, collapses) at t 0 + 3.2 ms into two sub-bubbles. Based on the above results, it can be concluded that R-T instability induces the deformation of large bubbles, and the viscous shear force strengthens bubble deformation and gives the bubbles a neck structure. Subsequently, turbulent pulsation shatters the neck structure and eventually decomposes the large bubbles.
The blue frames in Figure 14 show the third stage of bubble breakup, which is a multiple breakup process. This process is similar to the binary breakup, the difference being that the shearing-off phenomenon is rarely seen. It can be found that such an aggressive shearing-off tends to occur in the region with small vorticity fluctuations, and bubbles of large size are more likely to occur. The bubble breakup mechanisms reported in the current work are summarized in Figure 15.

Conclusions
Numerical simulation within the OpenFOAM framework was conducted to investigate the coupling mechanisms of bubble breakup in a venturi-type bubble generator. In order to track the two-phase interfaces and transient turbulence characteristics of the flow field, the volume of fluid method and large eddy simulation were combined in this research. Furthermore, a turbulence inlet database was established to ensure the reliability of the numerical solution and the economy of the calculation. Based on the simulation results under a variety of operating conditions, the following conclusions can be drawn.
(1) The flow field in the divergent section has obvious axial regional characteristics, which leads to the difference in bubble breakup patterns. The breakup patterns can be divided into steady-state disturbed breakup and unsteady-state disturbed breakup macroscopically. In a venturi-type bubble generator, the unsteady-state disturbed breakup pattern plays a dominant role in the generation of micro bubbles.
(2) The location of the blocking-vortex marks the transition point of bubble flow transport and breakup mechanisms, which have a certain controllability.
The forward movement of a blocking-vortex can enlarge the unsteady-state disturbed breakup region, and further improve the efficiency of generating microbubbles.
(3) Enhancing the inlet disturbance of a venturi channel is an effective way to advance the blocking-vortex. In the limited space of the divergent section, the steadystate disturbed breakup pattern occupies a greater region with the increase of gas volume fraction. However, at the same volume fraction of the gas, the liquid flow rate has little effect on the position of the blocking-vortex. (4) Interfacial instability plays an important role at each stage of bubble transport. Kelvin-Helmholtz instability causes bubbles to detach from the gas film during the steady-state disturbed breakup pattern, and Rayleigh-Taylor instability induces the deformation of detached bubbles in the unsteady-state disturbed region. (5) In the unsteady-state disturbed region, the viscous shear force strengthens the interfacial instability. Turbulent fluctuation and collision contribute to the shattering of defective bubbles. Additionally, the shearing-off process tends to occur in the region of small vorticity fluctuations, while it plays an insignificant role in creating large clusters of bubbles.
The present study provides a series of numerical investigations based on the VOF method for the whole area and enhances understandings of the mechanisms of bubble breakup in a venturi-type bubble generator. As only one geometry was tested in the present work, performing systematic tests (different structural parameters or more realistic operating conditions) is encouraged in the future to obtain more detailed and widely applicable flow mechanisms. On the other hand, the refinement of mesh or the optimization of numerical methods to obtain more precise output data are other issues needing to be addressed in future work.