Underwater explosion multiphase flow forecasting method based on level-set function

Abstract

The invention discloses an underwater explosion multiphase flow forecasting method based on a level-set function, and relates to the field of explosion numerical simulation, the method uses an interface distance function to capture a two-phase flow interface and determine a fluid coverage area, then uses a third-order TVD-Runge-Kutta method to solve an Euler control equation in the fluid coverage area to obtain a multiphase flow forecasting result, and finally, the multiphase flow forecasting result is obtained. And an updating formula of the interface distance function is given, and the time dispersion of the interface distance function and the time dispersion of the Euler control equation adopt different strategies, so that the overall solving calculation efficiency can be improved on the basis of achieving relatively good forecasting accuracy, thereby having relatively good performance in forecasting accuracy, precision, efficiency and stability. And a solid theoretical support can be provided for deep research of underwater explosion dynamics. The method is easy to popularize to accurate capture simulation of a complex interface topology motion process, so that an excellent solution is provided for engineering application.

Description

Technical Field

[0001] The present application relates to the field of explosion numerical simulation, and in particular to a method for predicting multiphase flow of underwater explosions based on a level-set function. Background Art

[0002] Underwater explosions have a wide range of applications in both the military field and in civil defense processes such as channel excavation and underwater protection. How to accurately predict the shock wave load and bubble movement process during underwater explosions has always been a difficult problem in the field of underwater near-free surface explosions.

[0003] There are many different models for numerical simulation of underwater explosion scenes, among which the most widely used are Lagrange method, ALE method and Front tracking method. These methods have their own characteristics, but they all have different degrees of shortcomings and limitations: Lagrange method is particularly suitable for calculating problems with small deformation and can provide clear material interfaces, but for large deformation fluid problems, the grid needs to be re-divided, and the grid itself may intersect and cause calculation interruption. When the deformation is large, the ALE method can be considered, but for multiple large deformation interfaces or when the interface undergoes topological changes, the grid movement is difficult to operate, which may cause errors and destroy the conservation of the numerical format. The Front tracking method usually uses the classic Euler method to calculate in areas other than the vicinity of the interface, and uses a specially designed method to calculate near the interface. Its disadvantage is that it is difficult to handle changes in the interface topological structure, and more computing resources are required as the number of interfaces increases. It is difficult for these existing numerical prediction models to achieve good performance in both prediction accuracy and prediction efficiency, and the prediction effect is not ideal. Summary of the invention

[0004] In view of the above problems and technical requirements, this application proposes a method for underwater explosion multiphase flow prediction based on level-set function. The technical solution of this application is as follows:

[0005] A method for predicting multiphase flow of underwater explosion based on level-set function, the method comprising:

[0006] At any nth time step, according to any grid G ​​along the i-th column in the x direction and the j-th column in the y direction in the computational domain, i,j The interface distance function φ at the nth time step i,j (x,n) is solved to obtain the motion interface Γ(n) and the fluid coverage area Ω at the nth time step f (n); where n, i, j are integer parameters, and x is the grid G i,j Any point within

[0007] The fluid coverage area Ω at the nth time step f (n) Using the third-order TVD-Runge-Kutta method to solve the Euler fluid governing equations, we can obtain the multiphase flow prediction results at the nth time step;

[0008] According to the grid G i,j The interface distance function φ at the nth time step i,j (x,n), using the second-order Runge-Kutta method to solve Get the grid G i,j The interface distance function φ at the n+1th time step i,j (x,n+1) and enter the forecast of the n+1th time step, where Δτ is the time step, sign() is the sign function, φ0=φ i,j (x,1).

[0009] Its further technical solution is to obtain the grid G i,j The interface distance function φ at the n+1th time step i,j (x,n+1) includes:

[0010] according to The interface distance function φ is calculated i,j (x,n+1), where and φ i,j Any parameter φ in (x,n) i,j have:

[0011]

[0012] in,

[0013]

[0014] and,

[0015]

[0016] Where Δx is the grid G i,j The width of the grid, Δτ = Δx / 2, ε = αΔx, α is a parameter.

[0017] Its further technical solution is:

[0018] in,

[0019]

[0020] For any parameters a and b, the function and

[0021] A further technical solution is to solve the motion interface Γ(n) and the fluid coverage area Ω at the nth time step. f (n) Includes:

[0022] Determine any grid G i,j The inner i,j The point x where (x,n)=0 is located on the motion interface, and the point x that makes φ i,j The point x where (x,n)>0 is located in the fluid coverage area. The motion interface Γ(n) of the nth time step is obtained by combining all the points x in all grids that are located on the motion interface. The fluid coverage area Ω of the nth time step is obtained by combining all the points x in all grids that are located in the fluid coverage area. f (n).

[0023] A further technical solution is to obtain the multiphase flow prediction result at the nth time step including:

[0024] The governing equations for Euler fluids On any grid Ω k Integrate and introduce the divergence formula, F = (F(U), G(U)) on the grid Ω i There are continuous first-order partial derivatives, grid Ω k The boundary of is composed of block-smooth directed closed curves. The unit external normal vector of any point on the boundary of the directed closed curve is Then there exists a semi-discrete finite volume format:

[0025]

[0026] Among them, U k is the grid Ω k The conserved variable U at any point x in the computational domain is ρ is the fluid density, u is the fluid velocity at point x along the x direction, v is the fluid velocity at point x along the y direction, p is the fluid pressure, E is the total energy, and t represents time;

[0027] The grid Ω k boundary Points on Split into the integral sum of each line segment on the directed closed curve, and simplify the semi-discrete finite volume format to obtain the semi-discrete equation:

[0028]

[0029] Among them, K kr is the grid Ω k The border The projection of the numerical flux on the rth line segment in the directed closed curve along its unit external normal vector, grid Ω k The border The directed closed curve of includes L line segments in total;

[0030] The semi-discrete equation is solved by the third-order TVD-Runge-Kutta method to obtain the fluid coverage area Ω at the nth time step. f The conservation variable U of (n) is used to obtain the multiphase flow prediction result at the n+1th time step.

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

[0032] The present application discloses a method for underwater explosion multiphase flow prediction based on level-set function. The method uses interface distance function to capture two-phase flow interface based on level-set method, solves Euler control equation in fluid coverage area using third-order TVD-Runge-Kutta method to predict multiphase flow, and gives time discretization and update formula of level-set function. Time discretization of level-set function and time discretization of Euler control equation adopt different strategies, which can improve the overall solution calculation efficiency on the basis of achieving better prediction accuracy, thus having better performance in prediction accuracy, precision, efficiency and stability, and can provide solid theoretical support for in-depth research on underwater explosion dynamics. The method can be easily extended to the accurate capture simulation of complex interface topological motion process, thus providing an excellent solution for engineering application. BRIEF DESCRIPTION OF THE DRAWINGS

[0033] Figure 1 It is a method flow chart of a method for underwater explosion multiphase flow prediction according to an embodiment of the present application.

[0034] Figure 2 It is a schematic diagram of the computational domain and the initial state of the shock wave in a test example of a double Mach reflection strong discontinuity problem.

[0035] Figure 3 yes Figure 2 In the example, the pressure wave front distribution predicted by the underwater explosion multiphase flow prediction method of the present application is compared with the reference simulation result of Titarev.

[0036] Figure 4 The present invention is a cloud diagram of the evolution of fluid density, fluid pressure and motion interface obtained by using the underwater explosion multiphase flow prediction method of the present invention to simulate an underwater near-free surface explosion evolution process.

[0037] Figure 5 The invention discloses an evolution cloud diagram of fluid density, fluid pressure and motion interface obtained by simulating and predicting the evolution process of a bubble water jet of an underwater near-wall explosion using the underwater explosion multiphase flow prediction method of the present invention. DETAILED DESCRIPTION

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

[0039] This application discloses a method for predicting underwater explosion multiphase flow based on level-set function. Please refer to Figure 1 The flowchart shown in the figure, the underwater explosion multiphase flow prediction method includes:

[0040] First, the computational domain is meshed. The computational domain in the underwater explosion scenario is a rectangular computational domain of a predetermined size. The coordinate system xy is established with the long side of the computational domain as the x direction and the wide side as the y direction. The computational domain is divided along the x direction and the y direction to obtain several rectangular grids. The grid along the i-th column in the x direction and the j-th column in the y direction in the computational domain is defined as grid G. i,j . Where i and j are both integer parameters.

[0041] This application uses the level-set function to predict multiphase flow and define any grid G i,j The interface distance function at the nth time step is φ i,j (x,n), where x is the grid G i,j The initial value of the integer parameter n is 1 at any point in the range.

[0042] When n=1, the interface distance function of each grid at the nth time step is initialized as φ i,j (x,n), and then according to any grid G i,j The interface distance function φ at the nth time step i,j (x,n) is solved to obtain the motion interface Γ(n) and the fluid coverage area Ω at the nth time step f (n) including:

[0043] Determine any grid G i,j The inner i,j The point x where (x,n)=0 is located on the motion interface, and the point x that makes φ i,j The point x where (x,n)>0 is located within the fluid coverage area, and the point x such that φ i,j The point x where (x,n)<0 is located in the gas coverage area. The motion interface Γ(n) of the nth time step is obtained by integrating all the points x in all grids that are located on the motion interface, and the fluid coverage area Ω of the nth time step is obtained by integrating all the points x in all grids that are located in the fluid coverage area. f (n), the gas coverage area Ω of the nth time step is obtained by integrating all points x in all grids within the gas coverage area. s (n), that is, the computational domain of the nth time step can be defined as:

[0044]

[0045] In addition, when determining the fluid coverage area Ω at the nth time step f (n), the fluid coverage area Ω at the nth time step f (n) The third-order TVD-Runge-Kutta method is used to solve the Euler fluid governing equation to obtain the multiphase flow prediction result at the nth time step.

[0046] The governing equations for Euler fluids are:

[0047]

[0048] Among them, the conserved variable at any point x in the computational domain is ρ is the fluid density, u is the fluid velocity at point x along the x direction, v is the fluid velocity at point x along the y direction, p is the fluid pressure, E is the total energy, and t represents time.

[0049] The above Euler fluid control equation is applied to any grid Ω k Integrate and introduce the divergence formula, F = (F(U), G(U)) on the grid Ω i There are continuous first-order partial derivatives, grid Ω k The boundary of is composed of block-smooth directed closed curves. The unit external normal vector of any point on the boundary of the directed closed curve is Then there exists a semi-discrete finite volume format:

[0050]

[0051] Among them, U k is the grid Ω k The conserved variable U at any point x in the

[0052] The grid Ω k boundary Points on The semi-discrete equation is obtained by splitting it into the integral sum of each line segment on the directed closed curve and simplifying the above semi-discrete finite volume format:

[0053]

[0054] Among them, K kr is the grid Ω k The border The projection of the numerical flux on the rth line segment in the directed closed curve along its unit external normal vector, grid Ω k The border The directed closed curve of includes L line segments. Numerical flux K kr Gaussian integral can be used, and the numerical flux K krThe construction method is the key to the spatial discretization of the control equation, and it is also the difference between different solution formats and their computational efficiency and accuracy.

[0055] Using the third-order TVD-Runge-Kutta method to solve the above semi-discrete equation, we can get the fluid coverage area Ω at the nth time step. f The conservation variable U of (n) is used to obtain the multiphase flow prediction result at the n+1th time step. When solving, the above semi-discrete equation can be rewritten as:

[0056]

[0057] The third-order TVD-Runge-Kutta method is as follows:

[0058]

[0059] Among them, Δt is the time step. The key to the above discretization method is to solve the numerical flux K in the normal direction of the grid boundary. kr For three-dimensional problems, it is necessary to solve the numerical flux in the normal direction of the three-dimensional unit surface. Therefore, this method is not limited to two-dimensional models and can be directly extended to three-dimensional models.

[0060] After completing the prediction of the nth time step, it is necessary to obtain the interface distance function of each grid in the next time step and enter the numerical prediction of the next time step. Theoretically, the interface distance function of each grid in the next time step can be calculated based on the unit external normal vector and curvature at each point on the moving interface Γ(n) of the nth time step. However, the interface distance function of the grid near the moving interface Γ(n) calculated in this way in the next time step does not strictly meet the properties defined by the signed distance function, which affects the multiphase flow prediction accuracy and interface capture accuracy. Therefore, in order to maintain the properties of the interface distance function of the next time step while not changing the interface position and sign defined by the interface distance function of the next time step, this application updates the arbitrary grid G ​​by solving the following function update formula i,j The interface distance function φ at the n+1th time step i,j (x,n+1):

[0061]

[0062] Where Δτ is the time step, which is independent of the time step Δt for solving the Euler fluid control equation and can be determined according to the grid size. For example, Δτ = 0.1Δx, where Δx is the grid G. i,j Width. φ0=φ i,j (x,1), sign() is the sign function, when φ0 is greater than 0, sign(φ0) = 1, otherwise sign(φ0) = -1.

[0063] When solving the above formula, the second-order Runge-Kutta method is used. Different time discretization accuracy is used for solving the Euler control equation. This can reduce the efficiency of computing resources and improve computing efficiency while ensuring the accuracy of the solution. The time discretization of the second-order Runge-Kutta method is as follows:

[0064]

[0065] in, is an intermediate variable in the solution process. and φ i,j Any parameter φ in (x,n) i,j have:

[0066]

[0067] Among them,

[0068]

[0069] So we can sort it out:

[0070]

[0071] in,

[0072]

[0073] In the above formula, Δx is the grid G i,j The width of the prism. ε=αΔx, α is a parameter, usually α=2.5.

[0074]

[0075]

[0076] in,

[0077]

[0078] In the above formula,

[0079]

[0080] The same calculation can be done. For any parameters a and b, the function φ i-1,j Represents the grid G i,j The adjacent grid G ​​on the left along the x direction i-1,j The interface distance function at the nth time step, φ i+1,j Represents the grid G i,j The adjacent grid G ​​on the right along the x directioni+1,j The interface distance function at the nth time step, φ i,j-1 Represents the grid G i,j The adjacent grid G ​​along the y direction i,j-1 The interface distance function at the nth time step, φ i,j+1 Represents the grid G i,j The adjacent grid G ​​below along the y direction i,j+1 Interface distance function at the nth time step

[0081] Each grid G ​​can be updated through the above formula i,j The interface distance function φ at the n+1th time step i,j (x,n+1) then enters the calculation of the next time step and repeats the above process for numerical prediction.

[0082] In a test example of a double Mach reflection strong discontinuity problem, please refer to Figure 2 , the computational domain is taken as [0,4]×[0,1], the number of grids in the computational domain is 960*240, the reflection wall is located at the bottom of the computational domain, an oblique strong shock wave with a Mach number of 10 is placed at x=1 / 6, y=0 and forms an angle of 60° with the x-axis, the bottom wall before x=1 / 6 adopts the accurate shock wave post-wave condition, and the other wall surfaces adopt the reflection boundary condition, such as Figure 1 For a detailed discussion of this model, please refer to the famous paper by Woodward and Colella. The initial conditions of the numerical simulation are obtained based on the normal shock wave relationship, and the relationship between the shock wave before and after is as follows:

[0083]

[0084] The underwater explosion multiphase flow prediction method of the present application is used for prediction and compared with the reference simulation results of Titarev. At T = 0.2s, the shock wave propagates to the lower right corner area in the calculation domain. The pressure wave front distribution predicted by the underwater explosion multiphase flow prediction method of the present application is as follows: Figure 3 As shown in (a) of the figure, the pressure wave front distribution shown by Titarev's reference simulation results is as follows Figure 3 As shown in (b), compared with Figure 3 It can be seen from (a) and (b) that the prediction results of the pressure wave front distribution of the present application are very close to the reference simulation results of Titarev, indicating that the present application has good reliability and accuracy in dealing with the propagation and reflection of shock waves.

[0085] In a simulation example of the evolution process of an underwater near-free-surface explosion, the computational domain is taken as [0,9]×[0,7], the free surface is located at y=3.5, the center of the charge is located at (4.5,3.0), and the radius is 0.1. The upper boundary and the left and right boundaries of the computational domain are all transmission boundaries, and the lower boundary is a solid wall sliding boundary, which is used to simulate near-free-surface shallow water explosions. The reflection of the pool bottom is considered as total reflection in the calculation. The underwater explosion multiphase flow prediction method of the present application is used for prediction, and the evolution cloud diagrams of the fluid density, fluid pressure, and moving interface at 0.56s, 1.14s, 2.23s, and 6.39s after the explosion are obtained as shown in the following figure. Figure 4 As shown by Figure 4 It can be seen that part of the shock wave generated at the beginning acts on the bottom of the pool and then reflects back to the gas-liquid motion interface, generating a compression wave inside the bubble and a rarefaction wave at the same time, while the other part is transmitted to the interface and reflects back to the rarefaction wave near the free surface, which in turn acts on the gas-liquid motion interface, transmits the rarefaction wave into the bubble, and reflects back into the water as a compression wave. The simulation of underwater explosions near the water surface has been a hot topic in the field of interface capture research in recent years. It can be seen from this example that the underwater explosion multiphase flow prediction method of this application can be well used for the prediction of such problems.

[0086] In a simulation example of the evolution process of an underwater near-wall explosion bubble water jet, a solid wall boundary is set at the left boundary at x=0 of the computational domain, a symmetric boundary is set at the lower boundary, and a non-reflecting boundary is set at the upper and right boundaries. There is a bubble of known radius at a certain distance from the solid wall boundary. When the early shock wave propagation stage and the bubble motion process are treated separately, the same two-dimensional axisymmetric model is used. Both stages start with an initial high-pressure, high-density bubble, but the bubble radius is different. The initial condition of this example is set to The underwater explosion multiphase flow prediction method of the present application is used to predict the evolution cloud diagram of the fluid density, fluid pressure and motion interface at multiple different typical moments after the explosion. Figure 5 As shown, it can be seen that the underwater explosion multiphase flow prediction method of the present application can also be well used for the prediction of such problems.

[0087] 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 method for predicting multiphase flow of underwater explosion based on level-set function, characterized in that: The underwater explosion multiphase flow prediction method comprises: At any nth time step, according to any grid G ​​along the i-th column in the x direction and the j-th column in the y direction in the computational domain, i,j The interface distance function φ at the nth time step i,j (x,n) is solved to obtain the motion interface Γ(n) and the fluid coverage area Ω at the nth time step f (n); where n, i, j are integer parameters, and x is the grid G i,j Any point within The fluid coverage area Ω at the nth time step f (n) Using the third-order TVD-Runge-Kutta method to solve the Euler fluid governing equations, we can obtain the multiphase flow prediction results at the nth time step; According to the grid G i,j The interface distance function φ at the nth time step i,j (x,n), using the second-order Runge-Kutta method to solve Get the grid G i,j The interface distance function φ at the n+1th time step i,j (x,n+1) and enter the forecast of the n+1th time step, where Δτ is the time step, sign() is the sign function, φ0=φ i,j (x,1).

2. The underwater explosion multiphase flow prediction method according to claim 1, characterized in that: Get the grid G i,j The interface distance function φ at the n+1th time step i,j (x,n+1) includes: according to The interface distance function φ is calculated i,j (x,n+1), where and φ i,j Any parameter φ in (x,n) i,j have: in, and, Where Δx is the grid G i,j The width of , ε = αΔx, α is a parameter.

3. The underwater explosion multiphase flow prediction method according to claim 2, characterized in that: in, in, For any parameters a and b, the function 4. The underwater explosion multiphase flow prediction method according to claim 1, characterized in that: Solve to obtain the motion interface Γ(n) and fluid coverage area Ω at the nth time step f (n) Includes: Determine any grid G i,j The inner i,j The point x where (x,n)=0 is located on the motion interface, and the point x that makes φ i,j The point x where (x,n)>0 is located in the fluid coverage area. The motion interface Γ(n) of the nth time step is obtained by combining all the points x in all grids that are located on the motion interface. The fluid coverage area Ω of the nth time step is obtained by combining all the points x in all grids that are located in the fluid coverage area. f (n).

5. The underwater explosion multiphase flow prediction method according to claim 1, characterized in that: The multiphase flow prediction results at the nth time step include: The governing equations for Euler fluids On any grid Ω k Integrate and introduce the divergence formula, F = (F(U), G(U)) on the grid Ω i There are continuous first-order partial derivatives, grid Ω k The boundary of is composed of block-smooth directed closed curves. The unit external normal vector of any point on the boundary of the directed closed curve is Then there exists a semi-discrete finite volume format: Among them, U k is the grid Ω k The conserved variable U at any point x in the computational domain is ρ is the fluid density, u is the fluid velocity at point x along the x direction, v is the fluid velocity at point x along the y direction, p is the fluid pressure, E is the total energy, and t represents time; The grid Ω k boundary Points on Split into the integral sum of each line segment on the directed closed curve, and simplify the semi-discrete finite volume format to obtain the semi-discrete equation: Among them, K kr is the grid Ω k The border The projection of the numerical flux on the rth line segment in the directed closed curve along its unit external normal vector, grid Ω k The border The directed closed curve of includes L line segments in total; The semi-discrete equation is solved by the third-order TVD-Runge-Kutta method to obtain the fluid coverage area Ω at the nth time step. f The conservation variable U of (n) is used to obtain the multiphase flow prediction result at the n+1th time step.

Images