1 Introduction

Bubbles can be seen everywhere in daily life and industrial production fields. Their typical phenomena include the cavitation and cavitation erosion in marine propellers, underwater weapons, liquid rocket pumps, microbubbles, etc. The study of bubble dynamics is one of the hotspots in the research fields of aeronautical and astronautical engineering, naval architecture and ocean engineering, and environmental engineering.

Bubbles are roughly divided into two types. One is gas bubbles, in which the gas is insoluble. The other one is vapor bubbles, in which the gas is generated by liquid-phase vaporization^[1]. The formation process of vapor bubbles in a liquid is referred to as bubble nucleation, and the phase change process dominated by pressure is referred to as cavitation, which is caused by a local drop in the pressure of the water stream. In a vortex commonly seen in the flow caused by the rotation of the fluid and the centrifugal force associated with the rotation, the pressure at the center is lower than that at the outer part. This makes it easy for cavitation to occur in the center of the vortex.

The process, by which microbubbles grow into macroscopic cavitation bubbles while moving along the suction side of a foil, is shown in Fig. 1^[2]. Cavitation damage often results in serious damage to the propeller blades, turbomachinery, and hydraulic equipment structure of a ship. The cavitation damage to the impeller blades of a pump is shown in Fig. 2^[3].

To discover the mechanism of cavitation erosion, the bubble dynamics near the boundary have been extensively studied in the last 50 years with research tools including the high-speed photography technology and the acoustic measurement. The destructive effects of cavitation can be roughly summarized as follows: the shock waves generated by bubble collapse and the high-speed jets pointing at the solid boundary.

The research on cavitation phenomena has a history of more than a hundred years. The numerical simulation of cavitation bubbles has developed rapidly in the last several years. Blake and Gibson^[4] studied the growth and collapse of vapor bubbles close to the horizontal free surface. Cerone and Blake^[5] studied the streamline, trace, and pressure nephogram characteristics of cavitation bubbles near a boundary based on the assumption that the liquid was incompressible and cohesionless and that the effects of the surface tension and gravity could be ignored. Chisum^[6] studied the dynamic behavior of underwater explosion bubbles. Chan et al.^[7] carried out a phenomenological study on the interactions between bubbles and structures. Pearson et al.^[8] studied the formation and evolution of different types of cavitation jets of bubbles near rigid boundaries. Afterwards, using some new methods such as the finite volume method (FVM), Chahine^[9] successfully introduced nonlinear compressible effects into the calculation of the dynamics of non-spherical bubbles. Zhang et al.^[10] studied the collapse of a central bubble affected by the surrounding bubbles under a high-pressure environment and with the interactions of multiple bubbles, i.e., the number of bubbles was greater than six. Koch et al.^[11] proposed a numerical method based on the FVM and the volume of fluid (VOF) method^[12], which could capture the shock waves caused by bubble collapse.

This study mainly focuses on the basic mechanism of bubble generation by cavitation. The objectives are to simulate the bubble interface of complex geometric shapes by using a compressible two-phase flow solver based on the bubble dynamic theory and the open source software OpenFOAM^[13] and to obtain a deep understanding of the basic characteristics of bubbles through a parametric study with the variations of the number of bubbles, the pressure, the radius and in the absence or presence of a free surface.

2 Governing equations and mathematical model

In this paper, a two-phase flow solver based on the FVM and the VOF method is used to simulate the bubble interfaces of complex geometric shapes. The FVM is the basis of fluid dynamics, and its basic process is to divide the entire computational domain into several subdomains, integrate the basic equations on each subdomain, transform the volume integrals into the surface integrals by means of Gauss's theorem, and then transform the differential equations into algebraic equations by using an interpolation scheme. The VOF method used in the OpenFOAM captures the interfaces by introducing the void fraction. In this study, the standard solver compressible inter foam (CIF) is modified for the simulation of bubble expansion and collapse.

2.1 VOF method

The VOF method is a numerical method for processing two-phase flow, which generally includes a gas phase and a liquid phase. The basic idea of the VOF method is to use a function to represent the fluid type, which is referred to as the "volume function" F. F and (1-F) represent the volume fractions of the gas and liquid in the computational domain, respectively. For a certain computational unit, (ⅰ) when F=0, it is completely filled with gas; (ⅱ) when F=1, it is completely filled with liquid; (ⅲ) when 0 < F < 1, it has a gas-liquid interface.

2.2 Compressible two-phase flow governing equations

In the two-phase flow solver used in this study, the continuity equation, the compressible momentum conservation equation, and the state equation of each phase are as follows:

(1)

(2)

(3)

where i=1 represents the liquid phase, and i=2 represents the gas phase. ρ_i is the density, and α_i is the phase fraction. U is the velocity. μ is the kinetic viscosity. p is the pressure. g is the acceleration of gravity. σ is the surface tension coefficient. is the outward normal vector of the interface between the two phases. ρ_oi is the initial density, and ψ_i represents the compressibility.

2.3 Algorithms and validation of the compressible two-phase flow solver

The PIMPLE algorithm^[14] in the FVM is used to solve the equations. It has a second-order spatial precision and a first-order time precision. The governing equations are discretised as follows: a first-order Euler implicit scheme is used in the time derivative, the Gaussian linear schemes are used in the pressure and velocity gradients, a Gaussian linear correction scheme is used in the Laplace operator, a first-order Gaussian upwind scheme is used in the convection term, a first-order non-orthogonal correction scheme is used in the surface normal gradient term, and a central difference scheme is used in the interface difference. The grid is a hexahedral structured grid.

The discretised coupling of the governing equations for the pressure-velocity coupling in the transient flow is obtained with the PISO^[15] algorithm as follows:

(ⅰ) Momentum prediction

The pressure field of the previous time step is used to solve the momentum equation and to obtain the predicted velocity field.

(ⅱ) Pressure solution

The predicted velocity field is used to solve the pressure equation and to obtain a first estimate of the new pressure field.

(ⅲ) Explicit velocity correction

The new pressure field is used to solve the momentum equation again and to correct the velocity field.

(ⅳ) Repeat the above steps until the allowable error is reached.

The difference between the implementation of the PISO algorithm in the CIF solver and that of the general PISO algorithm lies in the calculation of the α equation. When the momentum prediction is carried out, the new density and the surface density satisfying the mass continuity equation are needed, which can be obtained by using the old volume flux of the previous step so as to solve the α equation. Since the time step required by the compressed differential scheme^[16] is small enough, it is the same as that in other variable-density flow, and the change in the density field generated by the new volume flux obtained by the PISO algorithm is negligible. The solution process is as follows:

(Ⅰ) Initialize the variables.

(Ⅱ) Calculate the Courant number and adjust the time step.

(Ⅲ) Solve the α equation according to the volume flux of the old time level.

(Ⅳ) Obtain new estimates of the viscosity, density, and surface density according to the α value and the constitutive relationship.

(Ⅴ) Carry out momentum prediction according to the new values and continue the PISO algorithm.

(Ⅵ) If the final time has not been reached, enter the next time step, and return to Step Ⅱ.

All numerical calculations in this paper are carried out on the high-performance parallel cluster at School of Aeronautics and Astronautics, Dalian University of Technology. The operating system is CentOS6.4, and the OpenFOAM version is 2.1.1. The CIF solver used in this study is validated by means of the laser bubble collapse experimental data in Ref. [17]. The experimental data and the results of the CIF solver used in this study are plotted in Fig. 3. It can be found that the curves are coincident, verifying the effectiveness of the CIF solver used in this paper.

3 Numerical simulation of bubbles in a closed volume in the absence of a free surface

The computational domain of this study is a cube with a side length of 0.1 m. The boundary conditions are that the dynamic pressure on each of the six surfaces of the computational domain is 100 kPa and the velocity is zero. The acceleration of gravity is g=9.8 m·s^-2, and the initial velocity of the liquid is

3.1 Collapse of a single bubble

The initial conditions of a single bubble are as follows:

(C₁) The minimum grid size is 0.4 mm.

(C₂) The center of the bubble is located at (0.05 m, 0.05 m, 0.05 m).

(C₃) The initial radius of the bubble R₀ is 0.005 m.

(C₄) The pressure inside the bubble is 50 kPa.

(C₅) The pressure of the liquid is 100 kPa.

Figure 4 illustrates the three-dimensional (3D) physical process of the collapse and rebound of a single bubble under a bubble-liquid pressure ratio of 1:10.

The ratio t_Ray under the condition of a bubble-liquid pressure ratio of 1:2 in physical time is taken as the abscissa, and the ratio of the instantaneous radius to the initial radius is taken as the ordinate. The curves of the bubble radius under different initial pressure ratios over time are plotted in Fig. 5. It can be found that when the initial pressure of the bubble decreases, the collapse period decreases, and the volume compression ratio increases.

The ratio t_Ray under the condition of a bubble-liquid pressure ratio of 1:2 in physical time is taken as the abscissa, and the ratio of the instantaneous pressure to the initial pressure is taken as the ordinate. The curves of the pressure at the center of the bubble under different initial pressure ratios over time are plotted in Fig. 6. It can be found that when the initial pressure of the bubble decreases, the peak value of the pressure at the center of the bubble increases.

3.2 Simulation of the interactions between two bubbles 3.2.1 Breakdown caused by the bubble pressure ratio

In this section, we aim to find the bubble spacing and the bubble pressure ratio which can cause the occurrence of the breakdown phenomenon while keep the radii of the two bubbles constant. First, set up two bubbles with the radii of 5 mm, the minimum grid side length of 0.4 mm, and the velocity field of U=(0, 0, 0) within the domain and at the boundary. The bubble spacing and the bubble pressure ratio are adjusted by adjusting the initial positions and the initial pressures of the two bubbles. The initial conditions are found through a series of numerical experiments with the initial conditions as follows:

(c₁) The center of the high-pressure bubble is located at (0.042 5 m, 0.05 m, 0.05 m).

(c₂) The center of the low-pressure bubble is located at (0.057 5 m, 0.05 m, 0.05 m).

(c₃) The pressure of the high-pressure bubble is 400 kPa.

(c₄) The pressure of the low-pressure bubble is 50 kPa.

(c₅) The pressure of the liquid is 100 kPa.

Under such initial conditions, the high-pressure bubble will break down the low-pressure bubble.

The 3D physical process of the low-pressure bubble being broken down by the high-pressure bubble is shown in Fig. 7.

The radii of the two bubbles are extracted and plotted in Fig. 8. It can be found that under the condition of a large or small bubble pressure (higher or lower than the pressure of the liquid), when oscillation progresses, the volume of the low-pressure bubble first decreases and then increases, while the volume of the high-pressure bubble first increases and then decreases. When the volume of the high-pressure bubble reaches its maximum value, the low-pressure bubble is broken down, and a jet is generated. It can also be found in Fig. 8 that the oscillation of the radius curve of the low-pressure bubble is not a simple harmonic.

The diagram of the vorticity at the moment that the low-pressure bubble is broken down is shown in Fig. 9. The vorticity is found to be greater in the areas where the low-pressure bubble is broken down as well as where the high-pressure bubble is close to the low-pressure bubble. The white transparent part in the figure is the gas-liquid interface.

It can be seen from Fig. 10 that the pressure and the velocity are not uniform at the front and back sides of the low-pressure bubble. Since one side near the high-pressure bubble has a high pressure and a high velocity while the other side has a low pressure and a low velocity, a breakdown phenomenon is generated. In Fig. 10(a), the thick black line is the contour of the bubble, and the colors represent the velocity component in the x-direction.

3.3 Breakdown caused by the difference in the radii

In this section, our goal is to find the bubble spacing and the bubble radius ratio related to the occurrence of the breakdown phenomenon. We simulate the collapse process of bubbles with a pressure less than that of the liquid. First, we set up two bubbles with an initial pressure of 50 kPa, a minimum grid side length of 0.4 mm, and a velocity field of U=(0, 0, 0) within the domain and at the boundary, and adjust the bubble spacing and the bubble radius ratio by adjusting the initial positions and the initial radii of the two bubbles. The initial conditions are found through a series of numerical experiments, and are summarized as follows:

(D₁) The center of the large bubble is located at (0.044 5 m, 0.05 m, 0.05 m).

(D₂) The center of the small bubble is located at (0.054 5 m, 0.05 m, 0.05 m).

(D₃) The radius of the large bubble is 7.5 m.

(D₄) The radius of the small bubble is 2.5 m.

(D₅) The pressure of the liquid is 700 kPa.

Under such initial conditions, the phenomenon of large bubbles breaking down small bubbles occurs.

Figure 11 gives the collapse process of bubbles with a 1:3 radius ratio in the absence of a free surface.

It can be seen in Fig. 12 that the collapse period of the small bubble is smaller than that of the large one, and both bubbles depress another bubble in the volume expansion phase. Moreover, the large bubble has a greater effect on the small one, and it can break the small one down. When the small bubble is broken down in the fourth collapse period, it has its minimum volume. At that moment, its radius curve will no longer be a simple harmonic since its topology changes from a simple connection to a double one, which makes it possible to determine that the bubble is broken down directly when the equivalent radius curve is no longer a simple harmonic oscillation.

Next, we simulate the expansion process of bubbles with a pressure greater than that of the liquid. Except for an initial bubble pressure of 200 kPa, other initial conditions are exactly the same as those for the collapse simulation in the first half of this section. It can be seen from Fig. 13 that the small bubble is broken down at the first moment of its minimum volume. The large bubble and the small bubble have basically the same pulsation period, but a situation is completely different when the two bubbles collapse simultaneously. For two bubbles with the same pressure, when the initial bubble pressure is lower than the liquid pressure, the bubbles will have different sizes and thus different pulsation periods. Moreover, the larger the bubble radius is, the longer the period is. When the initial bubble pressure is higher than the liquid pressure, even if the sizes of the two bubbles are different, their pulsation periods will be the same.

3.4 Interactions among four bubbles

In this section, four bubbles are initially set symmetrically around the geometric center of the computational domain, and the morphological changes of the bubbles are studied by adjusting the bubble spacing and the pressure inside the bubbles. The initial radii of the bubbles are 5 mm, the minimum grid side length is 0.5 mm, the velocity field within the domain and at the boundary is U =(0, 0, 0), and the pressure of the liquid is 100 kPa. The initial conditions, which make the bubbles fuse and then separate, are found through a series of numerical experiments as follows:

(d₁) The center of Bubble 1 is located at (0.057 5 m, 0.05 m, 0.057 5 m).

(d₂) The center of Bubble 2 is located at (0.057 5 m, 0.05 m, 0.042 5 m).

(d₃) The center of Bubble 3 is located at (0.042 5 m, 0.05 m, 0.042 5 m).

(d₄) The center of Bubble 4 is located at (0.042 5 m, 0.05 m, 0.057 5 m).

(d₅) The bubble pressure is 600 kPa.

As shown in Fig. 14, in the bubble pulsation process, due to the interactions among the bubbles, the pressure of each bubble is not uniformly distributed on the sides near the center and far away from the center, and the average impulse of the bubble points towards the center, thereby causing the four bubbles to gradually draw closer and be broken down.

When the bubbles have been broken down but are not yet separated, since they are compressed, due to the conservation of momentum, the water flowing from the horizontal direction and meeting at the symmetry center will move along the upward and downward directions, which explains the upward and downward separation phenomena of the bubbles (see Fig. 15).

As shown in Fig. 16, after the bubbles are separated upward and downward, they will form two annular bubbles and have a different topology from that of the circular bubbles. When the annular bubbles are compressed, a local high pressure is generated at the symmetry center, making the liquid at this location flow along the upward and downward directions. The difference in the topology of the bubbles leads to different interactions among them, i.e., mutual attraction for the circular topology while mutual repulsion for the annular topology.

4 Numerical simulation of bubbles in a closed volume in the presence of a free surface

The effects of the free surface on the evolution and morphology of bubbles are investigated in this section. The geometric size of the 3D computational domain is Δx×Δy×Δz=0.1 m× 0.2 m× 0.1 m.

The boundary conditions are as follows: the pressure on each of the six surfaces of the computational domain is 100 kPa, and ∇α=0.

The initial conditions are as follows:

(e₁) The liquid is within the region of 0 m < y < 0.1 m, and the gas is within the region of 0.1 m < y < 0.2 m.

(e₂) The center of the bubble is located at (0.05 m, 0.09 m, 0.057 5 m) or (0.05 m, 0.075 m, 0.057 5 m).

(e₃) The initial radius of the bubble is R=5 mm.

(e₄) The bubble pressure is 400 kPa.

(e₅) The velocity field within the domain and at the boundary is U=(0, 0, 0).

(e₆) The pressure on each of the six boundaries of the computational domain is p_b=100 kPa, and the pressure of the liquid is p_w=100 kPa.

(e₇) The acceleration of gravity is g =9.8 m·s^-2.

4.1 Effects of the distance between a single bubble and the free surface on the interactions

By adjusting the distance of the bubble from the free surface, two working conditions are set up, in which the distance from the bubble center to the free surface is 0.025 m and 0.01 m, respectively. The simulation results are shown in Figs. 17 and 18, respectively. It is found that when the distances between the bubble and the free surface are large, the bubble has no effect on the free surface; when the bubble is closer to the free surface, the free surface forms a vertical upward water column. However, under both working conditions, the bubble is affected by the free surface to a certain extent, and the upper wall surface is depressed downward, thus forming a jet.

The curves of the radii of the bubbles under the aforementioned working conditions are shown in Fig. 19. It can be seen that the pulsation of the bubble near the free surface has a small amplitude and disappears quickly, demonstrating the effects of the free surface on a single expanded bubble.

4.2 Effects of the bubble pressure on the interactions

Based on the results for the bubble-liquid pressure ratio of 4:1 in the previous section, the calculation for a bubble-liquid pressure ratio of 8:1 is carried out in this section. Except that the initial pressure is changed to 800 kPa, other initial conditions are the same as those in the numerical example in which the distance between the bubble and free surface is 0.01 m. It can be found from Fig. 20 that under the working conditions of a bubble-liquid pressure ratio of 8:1, the water column on the free surface is high. The secondary water column is shown in Fig. 20(b).

The curves of the bubble radii over time under the working conditions of bubble-liquid pressure ratios of 4:1 and 8:1 are shown in Fig. 21. It can be found that the larger the bubble-liquid pressure ratio is, the larger the pulsation amplitude is, and the longer the pulsation phenomenon lasts.

4.3 Interactions between two bubbles and the free surface

In this section, two bubbles are placed at the same locations as those in the working conditions of the closed volume so as to study the effects of the presence of a free surface on the interactions of the bubbles. The initial conditions are as follows:

(E₁) The center of the high-pressure bubble is located at (0.042 5 m, 0.05 m, 0.05 m).

(E₂) The center of the low-pressure bubble is located at (0.057 5 m, 0.05 m, 0.05 m).

(E₃) The pressure of the high-pressure bubble is 1 000 kPa.

(E₄) The pressure of the low-pressure bubble is 50 kPa.

(E₅) The bubble radius is 5 mm.

As shown in Fig. 22, under the conditions of the presence of a free surface as well as a large pressure (higher than the pressure of the liquid) and a small pressure (lower than the pressure of the liquid) in the two bubbles, the low-pressure bubble is broken down first by the water flow from the direction of the wall (see Fig. 22(c)), after which the high-pressure bubble is broken down by the water flow from the direction of the free surface (see Fig. 22(d)). Finally, the low-pressure bubble fuses with the high-pressure bubble (see Fig. 22(f)).

4.4 Interactions among four bubbles

The initial radii of the bubbles are set to be 5 mm, the same as the radii of the four bubbles in the numerical example under the working condition of a closed volume. The initial conditions found through a series of numerical experiments are as follows:

(f₁) The center of Bubble 1 is located at (0.057 5 m, 0.08 m, 0.057 5 m).

(f₂) The center of Bubble 2 is located at (0.057 5 m, 0.08 m, 0.042 5 m).

(f₃) The center of Bubble 3 is located at (0.042 5 m, 0.08 m, 0.042 5 m).

(f₄) The center of Bubble 4 is located at (0.042 5 m, 0.08 m, 0.057 5 m).

(f₅) The bubble pressure is 600 kPa.

Under such conditions, the morphological evolution of the four bubbles in the presence of a free surface is similar to that of the four bubbles in a closed volume.

As in Fig. 23, under the condition that the bubbles are near the free surface, due to the expansion of the bubbles, a very weak water column is excited on the free surface above them. At the same moment, the bubbles fuse with each other in the expansion process, and the fused bubble is decomposed into two bubbles, i.e., a big bubble and a small bubble, under the pressure of the water column. Under the expansion action of the large bubble, the small bubble excites a finer and higher secondary water column in the center of the water column.

5 Conclusions

Bubbles and cavitation phenomena are hotspot topics in the research fields of aeronautical and astronautical engineering, naval architecture and ocean engineering, and environmental engineering. This study mainly focuses on the basic mechanism of bubble generation by cavitation. Under the working conditions with and without a free surface, we carry out numerical simulations on the interactions of a single bubble and multiple bubbles. Under the working condition with a free surface, for adjacent bubbles under equal pressures, the bubble morphology and dynamic characteristics are mainly affected by the relative sizes of the bubbles and the pressure of the liquid, while for adjacent bubbles under unequal pressures, the bubble morphology and dynamic characteristics are mainly affected by the difference between the pressures of the bubbles. Under the working condition with a free surface, for adjacent bubbles under equal or unequal pressures, the bubble morphology and dynamic characteristics are mainly affected by the pressure of the liquid between the upper surface of the bubbles and the free surface.

Although many interesting phenomena related to bubbles are obtained by adjusting the parameters, the mechanisms underlying such phenomena are not yet elucidated. An important physical phenomenon in the bubble collapse process is the shock wave, which has a strong destructive effect on the structure. Therefore, the use of a solver that can capture the shock wave to study the bubble problem will be a focus of our subsequent research in the future.