A high-precision numerical simulation method for strong shock compressible multiphase fluids

Abstract

The present application discloses a high-precision numerical simulation method for strong impact compressible multiphase fluid, which relates to the field of compressible multiphase flow. The method determines the flow field state quantity of the material interface based on the multi-medium Riemann problem according to the state of the material interface at the current time step, defines the virtual fluid state according to the flow field state quantity of the material interface determined by the solution, divides the entire flow field calculation domain into several single-medium flow calculation areas, and solves the flow field state prediction results of each single-medium flow calculation area at the next time step respectively, and then solves and updates the state of the material interface at the next time step based on the particle level set algorithm, and continues to solve at the next time step until all time steps are solved. This method can further improve the prediction accuracy of the material interface state, has engineering applicability and strong robustness, and has important engineering application significance.

Description

Technical Field

[0001] The present application relates to the field of compressible multiphase flow, and in particular to a high-precision numerical simulation method for strongly impacted compressible multiphase fluids. Background Art

[0002] For compressible multiphase flows formed under extreme conditions such as underwater explosions, they include the movement and interaction of material interfaces and strong shock waves, and have the characteristics of large density ratio, large pressure ratio, multi-scale space, multi-scale time, and high nonlinearity. The physical phenomena and processes are extremely complex. With the continuous improvement of research and application needs, high-precision numerical simulation of such strong impact compressible multiphase fluids has gradually become an indispensable research method and a very challenging topic in the fields of computational mathematics and computational fluid dynamics. For the numerical simulation of such strong impact compressible multiphase fluids, most of the simulation technologies currently used in the industry have problems such as accuracy and robustness. Summary of the invention

[0003] In view of the above problems and technical requirements, the applicant has proposed a high-precision numerical simulation method for strong impact compressible multiphase fluids. The technical solution of this application is as follows:

[0004] A high-precision numerical simulation method for strong impact compressible multiphase fluid, the high-precision numerical simulation method comprising:

[0005] Divide the computational grid in the flow field computational domain of the strong impact compressible multiphase fluid, and initialize the multiphase flow field distribution and the state of the material interface between the multiphase flow fields;

[0006] According to the state of the material interface at the current time step, the flow field state quantity of the material interface is determined based on the solution of the multi-medium Riemann problem;

[0007] The virtual fluid state is defined according to the flow field state quantity of the material interface determined by the solution, the entire flow field calculation domain is divided into several single medium flow calculation areas, and the flow field state prediction results of each single medium flow calculation area at the next time step are solved respectively;

[0008] According to the prediction results of the flow field state of each single-medium flow calculation area in the next time step, the state of the material interface in the next time step is updated based on the particle level set algorithm, and the step of determining the flow field state quantity of the material interface based on the multi-medium Riemann problem according to the state of the material interface in the current time step is executed again in the next time step until all time steps are solved.

[0009] The beneficial technical effects of this application are:

[0010] The present application discloses a high-precision numerical simulation method for strong impact compressible multiphase fluids. The method combines the modified virtual fluid method and the particle level set method, develops a method for constructing local multi-media Riemann problems more reasonably and accurately in multidimensional situations, and develops a more general and stable multi-media Riemann problem calculation method, further improving the prediction accuracy of the material interface state, thereby realizing a highly robust and high-precision simulation method with engineering applicability that can be applied to strong impact compressible multiphase fluids of underwater explosions. The method has important engineering application significance and can be easily extended to unstructured finite volume grids. BRIEF DESCRIPTION OF THE DRAWINGS

[0011] Figure 1 It is a method flow chart of a high-precision numerical simulation method according to an embodiment of the present application.

[0012] Figure 2 This is an example of using ray tracing to determine the spatial position P on both sides of the material interface. L and spatial position P R Schematic diagram of .

[0013] Figure 3 It is a flow chart of a method for determining the flow field state quantity of a material interface according to an embodiment of the present application.

[0014] Figure 4 It is a local schematic diagram of the material interface obtained by continuous reconstruction when using a grid-centered calculation grid.

[0015] Figure 5 It is a local schematic diagram of the material interface obtained by continuous reconstruction when using a grid-type calculation grid.

[0016] Figure 6 It is a schematic diagram of the flow field calculation domain in an example.

[0017] Figure 7 yes Figure 6 In the example, the high-precision numerical simulation method of this application is used to obtain t=1.8×10 -3 Schematic diagram of the flow field pressure and Y-direction velocity contour diagram at s. DETAILED DESCRIPTION

[0018] The specific implementation of the present application is further described below in conjunction with the accompanying drawings.

[0019] This application discloses a high-precision numerical simulation method for strong impact compressible multiphase fluids. Please refer to Figure 1 As shown in the flowchart, the high-precision numerical simulation method includes:

[0020] Step 1: divide the computational grid in the flow field computational domain of the strong impact compressible multiphase fluid, and initialize the multiphase flow field distribution and the state of the material interface between the multiphase flow fields, wherein the state of the material interface includes the shape and position of the material interface.

[0021] In this step, the numerical model of the strong impact compressible multiphase fluid is first determined. For the strong impact compressible multiphase fluid such as underwater explosion, the influence of factors such as medium viscosity, surface tension and gravity can be ignored when establishing the numerical model. Assuming that the media are immiscible and there is no mass conversion and heat transfer, the control equation of the strong impact compressible multiphase fluid can be expressed in the Cartesian coordinate system as follows:

[0022]

[0023] Among them, the expressions of U and F(U), G(U), and H(U) are:

[0024]

[0025] In the above formula, Ω is the flow field calculation domain, ρ is the medium density, p is the medium pressure, u, v, w are the velocity components of the medium along the x, y, and z directions of the Cartesian coordinate system, respectively, and E is the total energy per unit volume of the medium. e is the internal energy per unit mass.

[0026] In order to make the equation group closed, the state equation of the strong impact compressible multiphase fluid of underwater explosion can be expressed in a unified form as ρe=f(ρ)p+g(ρ), where f(ρ) and g(ρ) have different expressions according to different fluid media. The expressions under various fluid media are shown in the following table:

[0027]

[0028] Among them, γ, c 0 , 0 , B, A, ω, R 1 and R 2 are constants related to the physical properties of the fluid medium

[0029] After the numerical model of the strong impact compressible multiphase fluid is determined, the initialization of the multiphase flow field distribution and the state of the material interface can be completed according to the numerical model based on the actual working conditions.

[0030] Step 2: Determine the flow field state quantity of the material interface based on the state of the material interface at the current time step and the solution of the multi-medium Riemann problem. This includes:

[0031] (1) The ray tracing method is used to determine the normal direction of the material interface at the current time step, and the IThe initial state U of the left side of the fluid domain separated by the material interface at L And the initial state U of the right fluid domain R The spatial position of the left initial state U L and the initial state U on the right R The spatial position of the left initial state U is determined by interpolation. L and the initial state U on the right R The flow field state quantity.

[0032] Please refer to Figure 2 , determine the initial state U on the left L and the initial state U on the right R The flow field state quantity method includes: first, using the ray tracing method to determine the projection point A of the grid point A of the computational grid in the real fluid domain near the material interface along the normal direction of the material interface on the material interface. p , then for the projection point A p The predetermined width is extended along the normal direction of the material interface to the fluid domains on both sides of the material interface. The predetermined width is generally selected as 1.5 to 2 times the grid size. Then the spatial position P in the fluid domain on the left side of the material interface can be obtained. L And the spatial position P in the fluid domain on the right R .

[0033] Since the spatial position P L and spatial position P R Usually not on the grid point, so the spatial position P L and spatial position P R The flow field state quantity needs to be obtained by interpolating the flow field state quantity of the known grid points in the surrounding fluid medium domain. This application adopts an inverse distance weighted interpolation method, which is not limited to the topological form of the calculation grid, that is:

[0034] Using spatial position P L The flow field state quantity of the grid points in the surrounding fluid domain is calculated according to Perform inverse distance weighted interpolation to obtain the spatial position P L The flow field state quantity q L ,q i is the spatial position P L The flow field state quantity of any i-th grid point around, d i is the spatial position P L and the Euclidean distance between any i-th grid point around it. Using the spatial position P R The flow field state quantity of the grid points in the surrounding fluid domain is calculated according to Perform inverse distance weighted interpolation to obtain the spatial position P R The flow field state quantity q R ,qj is the spatial position P R The flow field state quantity of any j-th grid point around, d j is the spatial position P L and the Euclidean distance between any j-th grid point around it.

[0035] In the above interpolation process, the grid points around the spatial position are grid points whose Euclidean distances to the spatial position do not exceed a set value, and the number of grid points around the spatial position is determined according to actual conditions.

[0036] (2) According to the determined initial state U on the left L and the initial state U on the right R The flow field state quantity is constructed along the normal direction of the material interface at the current time step to construct the local one-dimensional compressible flow multi-medium Riemann problem as follows: In this local one-dimensional compressible flow multi-medium Riemann problem, t represents time, x represents the coordinate along the x direction of the Cartesian coordinate system, and x I represents the position on the material interface, x <x I Indicates that the position is located in the fluid domain on the left side of the material interface, x>x I It indicates that the position is located in the fluid domain on the right side of the material interface. ρ is the medium density, u is the flow velocity along the normal direction of the material interface, that is, the x-direction of the Cartesian coordinate system, p is the medium pressure, E is the total energy per unit volume of the medium and e is the internal energy per unit mass.

[0037] This step develops a method to construct the local one-dimensional compressible flow multi-medium Riemann problem more reasonably and accurately in multi-dimensional conditions, and this method is not limited to a certain grid topology form.

[0038] (3) Then, a universal and stable multi-medium Riemann problem calculation method is used to accurately solve the local one-dimensional compressible flow multi-medium Riemann problem constructed above, and the flow field state quantity of the material interface is obtained, including the pressure p of the material interface. * and normal velocity u * .

[0039] For underwater explosion strong shock compressible multiphase fluid, both sides of the material interface of the multi-medium Riemann problem are fluids (gas or liquid). Its wave system structure includes the left nonlinear wave corresponding to the left fluid domain, the right nonlinear wave corresponding to the right fluid domain, and the linear degenerate contact discontinuity wave corresponding to the middle material interface. The entire flow field calculation domain is divided into the left initial state U L , the initial state U of the right side of the fluid domain R The key to solving the above multi-medium Riemann problem is to solve the pressure p of the material interface.* , the pressure satisfies the following nonlinear equation:

[0040] Φ(p)≡f L (p,U L )+f R (p,U R )+u R -u L =0;

[0041] in, k is a variable parameter and can be L or R, p is a variable parameter, u L is the initial flow velocity of the fluid domain on the left, u R is the initial flow velocity of the fluid domain on the right, p L represents the initial medium pressure of the left fluid domain, p R represents the initial medium pressure of the fluid domain on the right, c is the speed of sound, ρ represents the density of the medium represented by the variable parameter p, ρ L represents the initial medium density of the fluid domain on the left.

[0042] Based on this, this embodiment provides a strategy combining internal and external iterations to achieve accurate solution of the multi-medium Riemann problem. The solution method includes: Figure 3 :

[0043] (a) Initialize the iteration step m = 0 and give the pressure p of the material interface at the 0th iteration step 0 .

[0044] (b) Pressure p based on the mth iteration step m , left initial state U L The flow field state quantity and the initial state U on the right R The flow field state quantity determines the nonlinear wave type of the left and right fluid domains of the material interface. The nonlinear wave type is shock wave or rarefaction wave, including the following four cases:

[0045] When satisfied , it is determined that the nonlinear wave type of the left fluid domain is shock wave, and the nonlinear wave type of the right fluid domain is shock wave.

[0046] When it is determined that p L <p R and , it is determined that the nonlinear wave type of the left fluid domain is shock wave, and the nonlinear wave type of the right fluid domain is rarefaction wave.

[0047] When it is determined that p R <p L and , it is determined that the nonlinear wave type of the left fluid domain is rarefaction wave, and the nonlinear wave type of the right fluid domain is shock wave.

[0048] When it is determined that p R <p L and , it is determined that the nonlinear wave type of the left fluid domain is sparse wave, and the nonlinear wave type of the right fluid domain is sparse wave.

[0049] in, The various parameters are explained above.

[0050] (3) According to the solution method corresponding to the nonlinear wave type of the fluid domain on the left, the f of the mth iteration step is obtained L (p m ,U L ), and the solution method corresponding to the nonlinear wave type of the fluid domain on the right side is used to obtain the f of the mth iteration step. R (p m ,U R ). For the fluid domain on the left side of the material interface and the fluid domain on the right side, solve f according to the solution method corresponding to the nonlinear wave type of the fluid domain. k (p m ,U k ) methods include:

[0051] When the nonlinear wave type in the fluid domain is a shock wave, the nonlinear algebraic equation of the shock wave branch is established based on the RH relationship of the internal energy e and combined with the state equation:

[0052]

[0053] The pressure p is obtained by iterative solution using the bisection method m The corresponding density ρ m , and then substitute Get f k (p m ,U k ).

[0054] When the nonlinear wave type in the fluid domain is a rarefaction wave, according to the pressure p m and the corresponding density ρ m The isentropic relationship combined with the above get The adaptive time step Runge-Kutta method is used to solve f k (p m ,U k ).

[0055] (4) According to f L (p m ,U L ) and f R (p m ,U R) Update the pressure p of the m+1th iteration step m+1 . Including determining the pressure p at the m+1th iteration step m+1 for:

[0056]

[0057] (5) When the pressure p in the m+1th iteration step is m+1 If the tolerance TOL is not met, set m = m + 1 and execute the pressure p based on the mth iteration step again. m , left initial state U L The flow field state quantity and the initial state U on the right R The steps of determining the nonlinear wave types of the left fluid domain and the right fluid domain of the material interface by using the flow field state quantity;

[0058] (6) When the pressure p in the m+1th iteration step is m+1 When the tolerance TOL is met, the pressure p of the m+1th iteration step is m+1 As the pressure p of the material interface * ;

[0059] (7) Pressure p based on material interface * according to Calculate the normal velocity u of the material interface * .

[0060] Step 3: Define the virtual fluid state according to the flow field state quantity of the material interface determined by solving, and divide the entire flow field calculation domain into several single medium flow calculation areas, including:

[0061] The constructed local one-dimensional compressible flow multi-medium Riemann problem Decoupling is a single medium Riemann problem corresponding to each single medium flow calculation area:

[0062]

[0063] Among them, the virtual fluid state in the single-medium Riemann problem and It is defined by the PGFM method based on the flow field state quantity of the material interface solved.

[0064] Step 4, respectively solve the flow field state prediction results of each single medium flow calculation area at the next time step. After decoupling the local one-dimensional compressible flow multi-medium Riemann problem into the single medium Riemann problem corresponding to each single medium flow calculation area, for each single medium flow calculation area, the corresponding single medium high-precision solver can be used for time advancement solution to obtain the flow field state prediction results of each single medium flow calculation area at the next time step. By combining the flow field state prediction results of each single medium flow calculation area, the flow field state prediction results of the entire flow field calculation domain at the next time step can be obtained.

[0065] Step 5: Based on the prediction results of the flow field state of each single medium flow calculation area at the next time step, the state of the material interface at the next time step is updated based on the particle level set algorithm. For strong impact compressible multiphase fluids, this application adopts a conservation form of Level-Set equation Among them, ρ is the medium density, u is the velocity, φ contains the Level-Set function value of each grid point and represents the signed distance function from the grid point to the material interface, and the zero isosurface φ = 0 is the material interface.

[0066] The solution process for the above Level-Set equation includes:

[0067] (1) First, the Level-Set equation The finite volume method is used to discretize the space into dΩ is the volume of the computational grid, S f To calculate the area of ​​the mesh surface, n f Calculate the unit external normal vector of the mesh surface.

[0068] (2) Using the third-order TVD Runge-Kutta method for time marching solution Get the Level-Set function value φ of each grid point.

[0069] In addition, due to the inherent dissipative effect of the numerical method, after a few time steps of time advancement, φ will no longer satisfy the signed distance property. Therefore, it is necessary to reinitialize it at each time step by solving the following initial value problem: x represents the spatial position and t represents the time step.

[0070] (3) Particle level set algorithm Based on the traditional level set algorithm, two types of Lagrangian particles with volume but no mass are introduced near the material interface to correct the moving material interface. The particles in the region of φ>0 are positive particles, and the particles in the region of φ<0 are negative particles. The particle migration equation is u(x p ) is x pThe velocity of the particle at is obtained by interpolating the velocities of the grid points around the particle. The particle migration equation is solved by time marching using the third-order TVD Runge-Kutta method. Get the position of all particles, and use the escaped particles to correct the Level-Set function value φ of each grid point. Define the position x of the grid point at any position x near the particle and the center of the particle p The distance is a local Level-Set function φ p (x) = s p (r p -|xx p |), r p is the radius of the particle, s p is the particle sign of the particle, the particle sign of the positive particle is positive, and the particle sign of the negative particle is negative. For each grid point, there may be multiple escaped particles around the grid point, so the Level-Set function value φ of the grid point is compared with the calculated correction value φ of each escaped positive particle. p (x) and compare the Level-Set function value φ of the grid point + Corrected to φ + =max(φ p (x),φ + ). The same is true for the escaped negative particles. After calculating the correction value φ of each escaped negative particle p (x), the Level-Set function value φ of the grid point - Corrected to φ - =min(φ p (x),φ - ). Using φ + and φ - The minimum value in is corrected to reconstruct the material interface.

[0071] In the actual calculation process, particles are usually attracted to a position very close to the material interface in the smooth area of ​​the material interface, while in the area where the material interface is greatly distorted, particles may escape from the material interface far away. In addition, sometimes the number of particles in some areas may be reduced due to the distortion of the material interface, or the particles in some areas may be too concentrated. Therefore, it is also necessary to reset the particles at each time step to ensure sufficient number of particles to correct the interface and reduce unnecessary particle tracking.

[0072] (4) Use the Level-Set function value φ of each grid point to continuously reconstruct the state of the material interface at the next time step. After the Level-Set function value φ of each grid point is determined by numerical methods, it is also necessary to continuously reconstruct the material interface based on the specific computational grid used. There are mainly two cases:

[0073] When the computational grid of the flow field computational domain is a grid-centred computational grid, make sure that when the Level-Set function values ​​φ of the grid points of two computational grids are opposite in sign, the shared grid edge or face of the two computational grids coincides with the material interface of the next time step. Please refer to Figure 4 .

[0074] When the computational grid of the flow field computational domain is a grid-type computational grid, the material interface passing through each computational grid is reconstructed according to the slicing characteristics of the Level-Set function value φ. For example, taking the computational grid as a two-dimensional triangular grid-type computational grid as an example, please combine Figure 5 , the coordinates x of the two end points of the material interface passed through in the calculation grid a and x b They are Among them, the coordinates of node A of the calculation grid are x A 、Level-Set function value is φ A , calculate the coordinates of node B of the grid to be x B 、Level-Set function value is φ B , the coordinate of the node C of the calculation grid is x C 、Level-Set function value is φ C .

[0075] Step 6, then in the next time step, execute again the step of determining the flow field state quantity of the material interface based on the state of the material interface in the current time step based on the solution of the multi-medium Riemann problem, that is, enter the next time step and repeat the process of steps 2-5 until all time steps are solved.

[0076] In one example, the evolution of a strong shock wave from an underwater explosion on the air-water interface was simulated, where the explosion was generated by a high-pressure bubble underwater. Figure 6 , the entire flow field calculation domain is X×Y∈[-5m,5m]×[-5m,5m], and the initial position of the free surface between air and water is y=2.5m. The high-pressure bubble with a radius of R1=1.0m is located at the origin (0.0,0.0). The parameters and initial conditions of the state equations of the three fluid media are shown in the following table. The non-reflecting boundary conditions are used on all sides of the calculation:

[0077] Fluid medium γ B(Pa) u(m / s) v(m / s) <![CDATA[ρ(kg / m 3 )]]> p(Pa) water 7.15 <![CDATA[3.309×10 8 ]]> 0.0 0.0 <![CDATA[1.0×10 3 ]]> <![CDATA[1.0×10 5 ]]> Explosion bubble 2.0 0.0 0.0 0.0 <![CDATA[1.27×10 3 ]]> <![CDATA[8.29×10 8 ]]> Air 1.4 0.0 0.0 0.0 1.0 <![CDATA[1.0×10 5 ]]>

[0078] Once the explosion begins, a very strong underwater shock wave is generated below the free surface and propagates radially outward with gradually decreasing intensity. After reaching the free surface, a reflected expansion wave and a weaker transmitted shock wave entering the air medium are formed. The reflected rarefaction wave further degenerates into a transmitted rarefaction wave and a reflected compression wave when it contacts the expanding bubble. The numerical simulation using the method of the present application obtains t = 1.8 × 10-3 The flow field pressure at s is Figure 7 As shown in (a), t = 1.8 × 10 -3 The Y-axis velocity contour diagram of s is shown in Figure 7 As shown in (b), from Figure 7 It can be seen that the method of the present application clearly captures all the main and local features of the underwater explosion flow field.

[0079] The above is only a preferred embodiment of the present application, and the present application is not limited to the above embodiments. It is understood that other improvements and changes directly derived or associated by those skilled in the art without departing from the spirit and concept of the present application should be considered to be included in the protection scope of the present application.

Claims

1. A high-precision numerical simulation method for strong impact compressible multiphase fluid, characterized in that: The high-precision numerical simulation method comprises: Divide the computational grid in the flow field computational domain of the strong impact compressible multiphase fluid, and initialize the multiphase flow field distribution and the state of the material interface between the multiphase flow fields; According to the state of the material interface at the current time step, the flow field state quantity of the material interface is determined based on the solution of the multi-medium Riemann problem; The virtual fluid state is defined according to the flow field state quantity of the material interface determined by the solution, the entire flow field calculation domain is divided into several single medium flow calculation areas, and the flow field state prediction results of each single medium flow calculation area at the next time step are solved respectively; According to the flow field state prediction results of each single-medium flow calculation area at the next time step, the state of the material interface at the next time step is updated based on the particle level set algorithm, and the step of determining the flow field state quantity of the material interface based on the multi-medium Riemann problem according to the state of the material interface at the current time step is performed again at the next time step until all time steps are solved; The method of determining the flow field state quantity of the material interface based on the state of the material interface at the current time step based on the solution of the multi-medium Riemann problem includes: using a ray tracing method to determine the normal direction of the material interface at the current time step, and determining the flow field state quantity of the material interface based on the solution of the multi-medium Riemann problem. I The initial state U of the left side of the fluid domain separated by the material interface at L And the initial state U of the right side of the fluid domain R The spatial position of the L and the initial state U on the right R The spatial position of the left initial state U is determined by interpolation. L and the initial state U on the right R The flow field state quantity; according to the determined initial state U on the left L and the initial state U on the right R The flow field state quantity is constructed along the normal direction of the material interface at the current time step to construct the local one-dimensional compressible flow multi-medium Riemann problem as follows: And solve it to get the flow field state quantity of the material interface including the pressure p of the material interface * and normal velocity u * ; where t represents time, x represents the coordinate along the x direction of the Cartesian coordinate system, and x I represents the position on the material interface, x <x I Indicates that the position is located in the fluid domain on the left side of the material interface, x>x I It indicates the fluid domain located on the right side of the material interface; ρ is the medium density, u is the flow velocity in the normal direction along the material interface, p is the medium pressure, E is the total energy per unit volume of the medium and e is the internal energy per unit mass; The method for solving the flow field state quantity of the material interface includes: initializing the iteration step m = 0, and giving the pressure p0 of the 0th iteration step of the material interface; based on the pressure p m , left initial state U L The flow field state quantity and the initial state U on the right R The flow field state quantity determines the nonlinear wave type of the left and right fluid domains of the material interface, and the nonlinear wave type is shock wave or rarefaction wave; according to the solution method corresponding to the nonlinear wave type of the left fluid domain, the f of the mth iteration step is obtained L (p m ,U L ), and the solution method corresponding to the nonlinear wave type of the fluid domain on the right side is used to obtain the f of the mth iteration step. R (p m ,U R ); according to f L (p m ,U L ) and f R (p m ,U R ) Update the pressure p of the m+1th iteration step m+1 ; When the pressure p of the m+1th iteration step m+1 If the tolerance TOL is not met, set m=m+1 and execute the pressure p based on the mth iteration step again. m , left initial state U L The flow field state quantity and the initial state U on the right R The flow field state quantity determines the nonlinear wave type of the left and right fluid domains of the material interface; when the pressure p of the m+1th iteration step m+1 When the tolerance TOL is met, the pressure p of the m+1th iteration step is m+1 As the pressure p of the material interface * ; Pressure p based on material interface * according to Calculate the normal velocity u of the material interface * ;u L is the initial flow velocity of the fluid domain on the left side of the material interface, u R is the initial flow velocity of the fluid domain on the left side of the material interface; The method for determining the nonlinear wave types of the left fluid domain and the right fluid domain of the material interface includes: when determining the nonlinear wave types satisfying When the nonlinear wave type of the left fluid domain is determined to be shock wave, the nonlinear wave type of the right fluid domain is determined to be shock wave; when it is determined to satisfy p L <p R And f L (p L ,U L )+f R (p L ,U R ) L -u R <f L (p R ,U L )+f R (p R ,U R ), the nonlinear wave type of the left fluid domain is determined to be a shock wave, and the nonlinear wave type of the right fluid domain is a rarefaction wave; when it is determined that p R <p L And f L (p R ,U L )+f R (p R ,U R ) L -u R <f L (p L ,U L )+f R (p L ,U R ), the nonlinear wave type of the left fluid domain is determined to be a rarefaction wave, and the nonlinear wave type of the right fluid domain is a shock wave; when it is determined that p R <p L and , the nonlinear wave type of the left fluid domain is determined to be a rarefaction wave, and the nonlinear wave type of the right fluid domain is determined to be a rarefaction wave; where, k is a variable parameter and can be L or R, p is a variable parameter, u L is the initial flow velocity of the fluid domain on the left, u R is the initial flow velocity of the fluid domain on the right, p L represents the initial medium pressure of the left fluid domain, p R represents the initial medium pressure of the fluid domain on the right, c is the speed of sound, ρ represents the density of the medium represented by the variable parameter p, ρ L represents the initial medium density of the fluid domain on the left.​​ 2. The high-precision numerical simulation method according to claim 1, characterized in that: Determine the initial state U on the left L and the initial state U on the right R The flow field state quantity methods include: Using a ray tracing method to determine the projection points of the grid points of the computational grid in the real fluid domain near the material interface along the normal direction of the material interface on the material interface; The projection point is extended along the normal direction of the material interface to the fluid domains on both sides of the material interface by a predetermined width to obtain a spatial position P located in the fluid domain on the left side of the material interface. L And the spatial position P in the fluid domain on the right R ; Using spatial position P L The flow field state quantity of the grid points in the surrounding fluid domain is calculated according to Perform inverse distance weighted interpolation to obtain the spatial position P L The flow field state quantity q L ,q i is the spatial position P L The flow field state quantity of any i-th grid point around, d i is the spatial position P L and the Euclidean distance between any i-th grid point around it; Using spatial position P R The flow field state quantity of the grid points in the surrounding fluid domain is calculated according to Perform inverse distance weighted interpolation to obtain the spatial position P R The flow field state quantity q R ,q j is the spatial position P R The flow field state quantity of any j-th grid point around, d j is the spatial position P L and the Euclidean distance between any j-th grid point around it.

3. The high-precision numerical simulation method according to claim 1, characterized in that: According to f L (p m ,U L ) and f R (p m ,U R ) Update the pressure p of the m+1th iteration step m+1 include:

4. The high-precision numerical simulation method according to claim 1, characterized in that: For the fluid domain on the left side of the material interface and the fluid domain on the right side, solve f according to the solution method corresponding to the nonlinear wave type of the fluid domain. k (p m ,U k ) methods include, k = L or k = R: When the nonlinear wave type in the fluid domain is a shock wave, the binary method is used to iterate and solve the problem. Get the pressure p m The corresponding density ρ m , substitute Get f k (p m ,U k ); When the nonlinear wave type in the fluid domain is a rarefaction wave, according to the pressure p m and the corresponding density ρ m The isentropic relationship combined with get The adaptive time step Runge-Kutta method is used to solve f k (p m ,U k ); Among them, p L represents the initial medium pressure of the left fluid domain, p R represents the initial medium pressure of the fluid domain on the right, c is the speed of sound, ρ L represents the initial medium density of the fluid domain on the left, ρ R represents the initial medium density of the fluid domain on the right.

5. The high-precision numerical simulation method according to claim 1, characterized in that: Methods for dividing the entire flow field calculation domain into several single-medium flow calculation regions include: The constructed local one-dimensional compressible flow multi-medium Riemann problem Decoupling is a single medium Riemann problem corresponding to each single medium flow calculation area: Among them, the virtual fluid state U in the single-medium Riemann problem is R g and U L g It is defined by the PGFM method based on the flow field state quantity of the material interface solved.

6. The high-precision numerical simulation method according to claim 1, characterized in that: Methods for solving and updating the state of the material interface at the next time step based on the particle level set algorithm include: Level-Set Equation The finite volume method is used to discretize the space into φ contains the Level-Set function value of each grid point and represents the signed distance function from the grid point to the material interface, ρ is the medium density, dΩ is the volume of the calculation grid, S f To calculate the area of ​​the mesh surface, n f is the unit external normal vector of the mesh surface, and u is the velocity; Time marching solution using the third-order TVD Runge-Kutta method Get the Level-Set function value φ of each grid point; The particle transport equation is solved by time marching using the third-order TVD Runge-Kutta method Get the positions of all particles and use the escaped particles to correct the Level-Set function value φ, u(x p ) is x p The velocity of the particle at The Level-Set function value φ of each grid point is used for continuous reconstruction to obtain the state of the material interface in the next time step.

7. The high-precision numerical simulation method according to claim 6, characterized in that: Methods for using the Level-Set function value φ of each grid point to continuously reconstruct the state of the material interface at the next time step include: When the computational grid of the flow field computational domain is a grid-centred computational grid, determining that when the Level-Set function values ​​φ of the grid points of two computational grids are opposite in sign, the common grid edge or face of the two computational grids coincides with the material interface of the next time step; When the computational grid of the flow field computational domain is a grid-type computational grid, the material interface passing through each computational grid is reconstructed according to the fragmentation characteristics of the Level-Set function value φ.

Images