Application of higher-order FV-WENO scheme to the interaction between shock wave and bubble

The high-order finite volume-WENO (Weighted Essentially Non-Oscillatory) scheme combines the finite volume method with the WENO method, allowing for high-order accuracy and accurate simulation of complex physical phenomena. It has advantages in handling shock waves and bubble interactions. The idea of this method is to discretize the physical equations into a set of conservation equations and use the WENO method to calculate the numerical fluxes on each finite volume cell. This approach is effective in dealing with problems such as shock wave propagation, bubble deformation, and evolution. In fluid dynamics simulation and research, the high-order FV-WENO scheme has significant application prospects. It can provide accurate numerical solutions and simulate complex physical phenomena, making it widely applicable in scientific research and engineering. In this study, we simulated the interactions between shock waves and single or double bubbles, obtaining the complete process of bubble collapse with clear bubble interfaces and strong program stability.


Introduction
The Weighted Essentially Non-Oscillatory (WENO) method is a modification of the Essentially Non-Oscillatory (ENO) method.The ENO method was initially proposed by Harten, Engquist, Osher, and others in 1987 [1] .It involves constructing interpolation polynomials directly on a given template of m nodes to obtain reconstructed values at the grid midpoints.However, directly using these reconstructed values is not ideal.The ENO method improves the accuracy of the interpolation method and achieves high resolution and non-oscillatory behavior by using sequentially expanding node templates and selecting templates based on the principle of minimizing the absolute values of the order differences.However, many intermediate calculation results are discarded in the implementation of the ENO scheme.To overcome the limitations of the ENO scheme, Liu [2] , Shu [3] , and others proposed the WENO scheme.The WENO scheme optimally combines the respective templates based on their smoothness to improve the overall accuracy of the discretization scheme and achieve high resolution.The WENO scheme can be divided into Finite Volume WENO (FV-WENO) [4] and Finite Difference WENO (FD-WENO) formats.The FV-WENO format requires the combination of WENO reconstruction with an exact or approximate Riemann solver, while the FD-WENO format requires the splitting of the convective flux function.These two types of WENO schemes have significant differences in their applications, with FV-WENO suitable for unstructured grids and FD-WENO limited to uniform structured grids.

( ) ( )
The main steps of the 5th-order WENO reconstruction and HLLC approximate Riemann solver are as follows: ① Compute the integral averages , i j q of the two-dimensional grid.
② Perform a 5th-order WENO reconstruction along the y-direction to obtain high-order approximations of the integral averages at the points located above and below , +1/ 2 i j q ± .This reconstruction is based on two templates, one biased to the left and the other biased to the right.
③ Similarly, perform a 5th-order WENO reconstruction along the x-direction to obtain high-order approximations of the integral averages +1/ 2, i j q ± at the points located to the left and right.This reconstruction is also based on two templates, one biased to the left and the other biased to the right.④ Use the reconstructed values obtained in the previous steps as input for the HLLC approximate Riemann solver to calculate the numerical flux.

Control equations
For the motion of multiphase fluid, assuming each phase is inviscid and compressible, neglecting surface tension and phase transitions between different phases, the governing equations of the fluid can be unifiedly represented as follows [11] : Where ρ represents the density of the fluid medium, p represents the pressure in the flow field, u, v represent the components of fluid velocity in x and y directions.E is the total energy per unit volume, ( ) The stiffened gas equation of state has been widely used to simulate multiphase flow problems involving shock wave propagation.The interface position of the multiphase flow is captured by different stiffened parameters within the fluid cells.The stiffened gas equation of state is defined as follows [9] : ( ) γ π = are two phase stiffness parameters.For ideal gases, =1.4,=0 γ π ; for water media, =5.5, =4.92 8 e Pa γ π . Since there is convective motion at the interface of multiphase flow, the fluid stiffness parameters a and b must satisfy the following convective equation [9][10] 0 The equation ( 4) is combined with equations (5) (7) to form the governing equations for compressible multiphase flow.Equation ( 4) is a conservation equation that ensures the conservation of important physical quantities such as density, momentum, and energy in the flow field.Equation (7), on the other hand, is a non-conservation equation primarily used to determine the position of the interface between the two phases.

Interaction between single bubble and shock wave
Richtmyer-Meshkov instability is a typical shock-driven bubble problem, where the interaction between the shock wave and the two-phase interface leads to the generation of vorticity, causing deformation and instability of the interface.Studying such problems helps to verify the accuracy of interface capturing in multiphase flow.Haas and Sturtevat (1987) [12] conducted an experimental study on the interaction between shock waves and cylindrical helium bubbles in a shock tube, which is a well-known Richtmyer-Meshkov instability problem.Now, we will simulate this experiment using a two-phase flow model.The calculation model is shown in figure 1.The initial parameters are as follows: ( )  The calculation domain is that the diameter of the cylindrical Helium bubble is 50 mm and the shock Mach number is 1.22.The initial conditions are obtained according to the normal shock wave relation, and the experimental results are compared with the numerical simulation results in figure 2. It can be seen from the diagram that the simulation results are close to the experimental results in several typical stages from the initial deformation of the bubble to the collapse of the bubble, and the capture accuracy of the two-phase interface is higher.

Interaction between double bubbles and shock wave
In order to investigate the interaction mechanism between shock wave and bubble group, the interaction process between shock wave and double bubble is studied on the basis of single bubble acted by shock wave in section 4.1.The initial parameters are the same as in the previous section, and the calculation model is shown in Fig. 3.

Parallel implementation of FV-WENO format
For the more complex processes involving shock waves and multiphase flow interfaces, especially in the numerical simulation of three-dimensional models, parallel processing is essential to adapt to highperformance computing.The programming is done using Fortran language in a message passing parallel programming environment (MPI).Let's take the example of the shock wave interacting with a single bubble problem shown in Section 3. First, the entire domain is divided into uniformly spaced sub-domains based on the number of processors used.In this case, for computational convenience, only the grid in the Y direction is divided.The computational tasks on each sub-domain are assigned to different processors.The main tasks include assigning initial conditions to each processor before the computation, and exchanging data between processors after the computation.The results obtained from the parallel program are identical to those obtained from the serial program, proving the correctness of the parallel transformation.
Table 1 lists the results of the calculations performed on a personal PC for the problem discussed in Section 4. The CPU used is an Intel(R) Core (TM) 2 Quad with a frequency of 2.33GHz and 3.25GB of memory.The table compares the average time for a single time step calculation, the speedup, and the parallel efficiency for different numbers of threads.As seen from Table 1, even for small-scale problems, the parallel efficiency reaches over 80%.This demonstrates that the parallel performance of the high-order FV-WENO scheme is highly satisfactory.
In summary, for the simulation of complex processes and three-dimensional models involving shock waves and multiphase flow interfaces, parallel processing is crucial for achieving highperformance computing.The use of MPI and Fortran programming language allows for efficient parallelization, leading to significant speedup and improved computational efficiency.

Conclusion
The high-order finite volume-WENO (Weighted Essentially Non-Oscillatory) scheme is a powerful numerical method that combines the finite volume method with the WENO method.This scheme allows for high-order accuracy and accurate simulation of complex physical phenomena, making it particularly suitable for handling shock waves and bubble interactions.The main idea behind this method is to discretize the physical equations into a set of conservation equations and utilize the WENO method to calculate the numerical fluxes on each finite volume cell.By doing so, we can effectively deal with problems such as shock wave propagation, bubble deformation, and evolution.In the field of fluid dynamics simulation and research, the high-order FV-WENO scheme holds great promise.It not only provides accurate numerical solutions but also enables the simulation of complex physical phenomena.As a result, it finds wide applications in scientific research and engineering.In this study, we focused on simulating the interactions between shock waves and single or double bubbles.Through this simulation, we were able to capture the complete process of bubble collapse with clear bubble interfaces and strong program stability.
the specific internal energy of the fluid per unit mass.The parameter n is related to the model being used: when n = 1, it represents a two-dimensional planar model; when n = 2, it represents a two-dimensional axisymmetric model, where the velocity v represents the radial velocity component.When n = 1 and n = 0, it degenerates into a one-dimensional planar model; when n = 2 and n = 0, it degenerates into a one-dimensional spherically symmetric model.

4 Figure 1 .
Figure 1.Schematic diagram of the interaction between shock wave and Helium bubble in air.

Figure 2 .
Figure 2. Cylindrical Helium bubble process under shock wave action (density cloud)

Figure 3 .Figure 4 .
Figure 3. Schematic diagram of the interaction between shock wave and double Helium bubble in air.The calculation domain is [ ] [ ] 0, 0.325 0, 0.178 × , in which the diameter of the cylindrical Helium bubble is the same as the shock wave Mach number above.Figure 4 shows the density cloud diagram of the flow field of t = 32, 245, 427, and 983us starting from the shock wave at the time of the shock wave.

Table 1 .
List of program parallel computing data