Numerical investigation of partial cavitation in a Venturi tube by Eulerian-Lagrangian multiscale modelling

The cavitating flow in a Venturi tube has a complex flow structure. In this paper, the partial cavitation in an axisymmetric Venturi tube, dominated by bubbly shock, is investigated by utilizing Eulerian-Lagrangian multiscale modelling. The unsteady partial cavitation is simulated by Large Eddy Simulation coupled with the Volume of Fluid and the Discrete Bubble Model. The numerical results show a well agreement with the high-speed photography. Through a comprehensive analysis of the development and transportation of macro cavities and micro bubbles, a detailed explanation of the cavity shedding process triggered by bubbly shock is provided. The results highlight the precise capability of the multi-scale method in elucidating the intricate flow field induced by partial cavitation. The findings may pave the way for the further investigations on the underlying mechanisms of partial cavitation, fostering a deeper understanding of bubbly shock dominated cavitating flow.

The existing experimental phenomena generally show that there are two mechanisms leading to cavity shedding: the re-entrant jet mechanism and the bubbly shock mechanism.The former usually occurs in the development of attached cavity.The adverse pressure gradient existed in the cavity closure region drives fluids to move upstream under the attached cavity [8].Brunhart et al. [9] reliably realized the shedding mechanism of the re-entrant jet in the axisymmetric convergent divergent Venturi nozzle.Trummler et al. [10] conducted further research on the re-entrant jet driving mechanism.
Bubbly shock, also known as condensation shock, usually occurs under the condition of flow with small cavitation numbers, that is, the attached cavity extends further downstream.The detached cavity forms cloud cavitation and creates a series of pressure waves after collapse.Due to the extremely small local sound velocity in the attached cavity, the pressure waves can become shock waves.This

Experimental method
The Venturi experimental system is illustrated in figure 1 (a).This system drives water circulate in the pipe through a centrifugal pump.The material of Venturi tube is polymethyl methacrylate (PMMA).The camera used in the visualization experiment is Phantom VEO 1310, combined with Tokina's atx-i 100 mm F2.8 FF lens. Figure 1    where is the velocity vector field, is time, is pressure, is the mixed density, is the gravitational acceleration, is the stress tensor.

Governing equations
The VOF method is used to accurately reproduce the interface between water and vapor.The transport equation of the vapor volume fraction is as follows: (3) where is the mass source term, which indicates the rate of condensation and evaporation.
The mass transfer in the VOF is determined by the simplified Rayleigh-Plesset equation.The mass transfer rate without the second-order term and surface tension is expressed as follow: (4) where is the saturated pressure, is the bubble radius: (5)

Discrete bubble model
In order to achieve multi-scale simulation, small-scale bubbles are simulated according to the discrete bubble model developed by Li et al. [15,16].Newton's second law is used to track each discrete bubble. ( The term on the right represents the force acting on it.The drag force is as follows: is the relative Reynolds number which is expressed as： The additional force, comprising of the virtual mass force and the pressure gradient force ( and ), is expressed as: (9) (10)

Simulation setup
The three-dimensional modeling of the axisymmetric convergent divergent Venturi is conducted by using SolidWorks software.Figure 2 (a) is a schematic diagram of the computation domain, and figure 2 (b) presents the size of each section of the Venturi tube.The Venturi length is 128 mm, the throat diameter is 4 mm, the throat length is 2 mm, and the diameter of the upstream and downstream sections is 25 mm.The grid of convergent divergent tube section is shown in figure 2 (c).Hexahedral grid is established, and the grid of cavitation area is encrypted.The height of the first layer grid is specified as 2×10 -5 m to ensure that the Y + of the wall is close to 1.In order to ensure the accuracy of the simulation calculation, the independence of the grid is verified in this paper.Based on the Venturi tube calculation model, four kinds of numerical grids are set, which are 5.9×10 5 , 1.5×10 6 , 2.3×10 6 , 4.0×10 6 , respectively.The pressure of the 375 mm section in front of the Venturi tube with different grids are extracted, as shown in figure 3. From the pressure distribution along the axis, it can be seen that with the increase of the number of grids, the pressure distribution deviation becomes smaller.When the number of grids is greater than 1.5 × 10 6 , the deviation is within the acceptable range.The number of grids of 2.3×10 6 is used for simulation.All walls are specified as no-slip boundary.The inlet velocity and the outlet pressure are set to 0.56 m/s, and 96,954 Pa, respectively.According to the formula , the cavitation number is 0.39, and the numerical simulation and experimental comparison are carried out under this operating condition.The saturated vapor pressure is set to 3500 Pa.The total calculation time is 0.0425 s, and the time step is set to 2.5×10 −6 s.The convergence criterion of all equations is 1×10 −4 , and the maximum number of iterations per time step is 25.For the discrete bubble modeling, a vapor core with an initial diameter of 2×10 -6 m is randomly added to the area.The number density of 1×10 5 is used to distribute the nuclei in the computational domain.Figure 4 shows that the Y + distribution of the wall is about 6 in the cavity closure area and about 1 in the downstream region.

Validation of numerical simulation results
The multi-scale modeling framework is used to simulate the cavitation state of cloud cavitation.Figure 5 is the comparison between the transient simulation results of cavity evolution and the experimental observation results under the typical period σ = 0.39.The numerical results of the transition of the partial cavitation including vortex shedding and collapse.Since the cloud cavitation, composed of numerous tiny bubbles, collapses in the downstream area, the conditions such as temperature and pressure change drastically, therefore, it is difficult to ensure that the predicted cavity is fully consistent with the experiment in space and time.Nonetheless, the numerical simulation still showed a reasonable agreement with the high-speed photography.
Figure 6 demonstrates that the vortex cloud shedding process is reproduced from the perspective of axisymmetric cross section.From a microscopic perspective, it can be observed that cavitation bubbles are formed when the local pressure falls below the saturation pressure.Subsequently, as the bubbles flow into high-pressure regions it undergoes a collapse process.Because the VOF is employed to accurately capture the macroscopic cavity and the VOF-DBM conversion is considered, there are always discrete bubbles in pure water and vapor.The multi-scale model considers large cavities and microbubbles in two different equations respectively, which can obtain more detailed characteristics and fully reveal cavitation phenomena such as large interfaces and microbubbles on an order-of-magnitude scale.

Morphological characteristics of the cavity
Figure 7 shows the cavity pattern in the axisymmetric Venturi under different cavitation numbers.Three typical moments under three representative cavitation numbers are selected, and the time intervals of these conditions are subsequent different.A large-scale cloud cavitation phenomenon begins to appear, and the cavity length begins to change dramatically.After the cloud cavitation disappears, it will grow again until the next shedding occurs.As the cavitation number decreases, there is an observable increase in the length of the cavity.At σ = 0.30, it can be observed that the maximum length of the cavity remains relatively consistent, and the shedding process of cavitation clouds exhibits remarkable similarity when compared to the case at σ = 0.39.When σ = 0.39, the cavity continues to develop, the cavitation cloud changes more violently, and the large-scale cavity in the downstream collapses, forming pressure waves.They propagate upstream and become bubbly shocks in the attached cavity region.This paper mainly aims at the phenomenon of cavity shedding induced by bubbly shock wave, and deeply studies the shedding mechanism by numerical simulation.

Flow characteristics of cloud cavitation
The quasi-periodic development and shedding processes of the partial cavitation along the wall of the diffusion section can be clearly observed by numerical simulation (figure 8).Once the attached cavity grows to its maximum length, it begins to detach.The vapor-liquid interface becomes sharp, and large shedding clouds move with the mainstream towards the downstream high-pressure zone.As shown in figure 9, brighter regions with high density variation can be seen, especially at T = t0 + 4 ms, indicating that the bubbly shock existed in the attached cavity greatly compresses the vapor-liquid mixture after the collapse of the cloud cavity.Figure 10 (a) illustrates the temporal evolution of the vapor distribution of the middle plane of the Venturi at σ = 0.39, where the blue line is the original data and the orange is the smoothed data which are obtained by the Savitzky-Golay method and make it easier to identify different cavitation stages.The power spectral density (PSD) of the variation of vapor area of the middle plane shows a maximum energy frequency of 117.36 Hz (figure 10 (b)).This observation suggests that the shedding process of the cavity is rather rapid and intense.The one-dimensional Rankine-Hugoniot jump condition describes the variation of the physical quantities before and after a shock wave, which can be used to calculate shock wave properties if the density and pressure are known, e.g., the propagation velocity [12].This condition is mathematically expressed by equation ( 11): (11) where s represents the wave propagation speed, and u represents the axial velocity.The density, axial velocity, and pressure values presented in table 1 are taken at the same position before and after the impact at T = t0+3 ms and T = t0+3.25 ms (figure 11).The calculated propagation velocities of the bubbly shock are 4.29 m/s and 4.38 m/s, exhibiting a remarkable proximity to the velocity of 4.51 m/s depicted in figure 12. Detailed description of this can be found in the works by Brunhart [9] and Budich [17].Table 1.Flow properties before and after the shock wave for the 1D Rankine-Hugoniot verification, obtained from the location shown in figure 10.

Conclusion
In this paper, a new multi-scale model, combining an interface model and a coupled Lagrangian discrete bubble tracking model (VOF-DBM), is adopted for the partial cavitation in an axisymmetric Venturi tube.The cavitation characteristics of both large cavities and microbubbles are characterized.
The simulation results exhibit a commendable concordance with the experimental observations, validating the correction of the multi-scale model.When σ = 0.39, the change of cavitation volume has well quasi-periodic behaver, and the frequency is approximately 117 Hz.According to the change of cavity volume, the occurrence of different stages of cloud cavitation can be clearly distinguished.
The bubbly shock can be well captured by the multi-scale model.At σ = 0.39, the cloud cavitation is dominant, and the sharp cavity trailing edge can be observed.The multi-scale model can well describe generation, growth and collapse of bubbles in the Venturi tube, the generation of large cavities, the formation of bubble clouds by bubbly shock, and the recovery of new cavities.
(b) illustrates the high-speed photography.The water and PMMA are demonstrated as black, and the cavitation clouds are white, which can be used to determine the shape and position of the cavity (figure 1 (c)).
In this work, LES is utilized as the turbulence model.Using the VOF method, water and vapor are considered impermeable to each other.The governing equations are the compressible Navier-Stokes equation of the VOF and the transport equation of the vapor mass fraction.The CFD simulations were run in ANSYS Fluent.The equations for solving the conservation of mass and momentum are as follows:

Figure 2 .
Figure 2. (a) Illustration of the physical domain and boundary conditions, (b) Venturi geometry, and (c) grid.

Figure 3 .
Figure 3. Absolute pressure distribution along the axis of the Venturi tube.

Figure 4 .
Figure 4. Distribution of Y + of the Venturi wall.

Figure 5 .
Figure 5. Validation of the numerical simulation by comparing the experimental cavity morphology at σ = 0.39.

Figure 8 .
Figure 8. Variation of the vapor volume fraction distribution of the Venturi middle plane at σ = 0.39.

Figure 9 .
Figure 9. Variation of the axial density gradient distribution of the Venturi middle plane at σ = 0.39.Figure10(a) illustrates the temporal evolution of the vapor distribution of the middle plane of the Venturi at σ = 0.39, where the blue line is the original data and the orange is the smoothed data which are obtained by the Savitzky-Golay method and make it easier to identify different cavitation stages.The power spectral density (PSD) of the variation of vapor area of the middle plane shows a maximum energy frequency of 117.36 Hz (figure10 (b)).This observation suggests that the shedding process of the cavity is rather rapid and intense.

Figure 10 .
Figure 10.(a) Variation of vapor area of the Venturi middle plane, and (b) the corresponding power spectral density.

DensityFigure 11 .
Figure 11.Data collection position of the flow properties before and after the shock wave for the 1D Rankine-Hugoniot verification.

Figure 12 .
Figure 12.Variation of the shock wave position in 1 ms.