Multiphase flow prediction method for underwater explosion based on level-set function

Abstract

This application discloses a method for underwater explosion multiphase flow prediction based on level-set function, which relates to the field of explosion numerical simulation. The method uses an interface distance function to capture the two-phase flow interface and determine the fluid coverage domain. Then, the third-order TVD-Runge-Kutta method is used to solve the Euler governing equation within the fluid coverage area to predict the multiphase flow. The updated formula of the interface distance function is given. The time discretization of the interface distance function and the time discretization of the Euler governing equation adopt different strategies. It can improve the overall solution calculation efficiency on the basis of achieving good prediction accuracy, thereby having good performance in prediction accuracy, precision, efficiency and stability. It can provide solid theoretical support for in-depth research on underwater explosion dynamics. The method can be easily extended to the precise capture and simulation of complex interface topological motion processes, thus providing an excellent solution for engineering applications.

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 military and 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] Numerical simulations of underwater explosions employ a variety of methods. The most widely used include the Lagrange method, the ALE method, and the front tracking method. While each method has its own unique characteristics, they all suffer from varying degrees of shortcomings and limitations. The Lagrange method is particularly well-suited for computations involving small deformations, providing clear material interfaces. However, for large deformation fluid problems, mesh re-meshing is required, which can lead to computational interruptions caused by mesh self-intersection. The ALE method can be considered for large deformations, but mesh movement is difficult to implement for multiple large deformation interfaces or when interfaces undergo topological changes, potentially introducing errors and violating the conservation properties of the numerical scheme. The front tracking method typically uses the classic Euler method for calculations outside the vicinity of interfaces, employing a specially designed method near interfaces. However, its drawback is that it struggles to handle variations in interface topology and requires more computational resources as the number of interfaces increases. Existing numerical prediction models struggle to achieve both high accuracy and efficiency, resulting in suboptimal prediction results. Summary of the Invention

[0004] In response to the above-mentioned problems and technical needs, this application proposes a method for predicting underwater explosion multiphase flow 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 along the i-th column in the x-direction and the j-th column in the y-direction in the computational domain Interface distance function at the nth time step Solve to get the motion interface of the nth time step and fluid coverage area ; Among them, n, i, j are integer parameters, It's a grid Any point within

[0007] Fluid coverage area at the nth time step 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;

[0008] According to the grid Interface distance function at the nth time step , using the second-order Runge-Kutta method to solve , get the grid Interface distance function at the n+1th time step And enter the forecast of the n+1th time step, where is the time step, is a sign function, .

[0009] Its further technical solution is to obtain the grid Interface distance function at the n+1th time step include:

[0010] according to Calculate the interface distance function , among which, for and Any parameter in have:

[0011]

[0012] in,

[0013]

[0014] and,

[0015]

[0016] in, It's a grid The width, , , As a parameter.

[0017] Its further technical solution is: , ;

[0018] in,

[0019]

[0020]

[0021] For any parameter and ,function ,and .

[0022] Its further technical solution is to solve the motion interface of the nth time step and fluid coverage area include:

[0023] Determine any grid Internal point Located on the motion interface, determine point Located in the fluid coverage area, all points on the moving interface in all grids are integrated Get the motion interface of the nth time step , all points in all grids that are within the fluid coverage area Get the fluid coverage area at the nth time step .

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

[0025] Euler fluid governing equations In any grid Integrate above and introduce the divergence formula, In the grid There are continuous first-order partial derivatives, grid The boundary of the directed closed curve 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 scheme:

[0026]

[0027] in, It's a grid Any point inside Conserved variables at , any point in the computational domain Conserved variables at , , , is the fluid density, Yes The fluid velocity along the x direction at , Yes The fluid velocity along the y direction at , is the fluid pressure, is the total energy, Indicates time;

[0028] The grid 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:

[0029]

[0030] in, It's a grid The border The first of the directed closed curves The numerical flux on the line segment is projected onto its unit external normal vector, the grid The border The directed closed curves of line segments;

[0031] The semi-discrete equation is solved using the third-order TVD-Runge-Kutta method to obtain the fluid coverage area at the nth time step. Conserved variables To obtain the multiphase flow prediction results at the n+1th time step.

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

[0033] This application discloses a method for underwater explosion multiphase flow prediction based on a level-set function. This method uses an interface distance function based on the level-set method to capture the two-phase flow interface. The third-order TVD-Runge-Kutta method is used to solve the Euler control equation in the fluid coverage area to obtain the multiphase flow prediction results. The time discretization and update formula of the level-set function are given. The time discretization of the level-set function and the time discretization of the Euler control equation adopt different strategies. It can improve the overall solution calculation efficiency on the basis of achieving good prediction accuracy, thereby having good performance in prediction accuracy, precision, efficiency and stability. It can provide solid theoretical support for in-depth research on underwater explosion dynamics. This method can be easily extended to the precise capture and simulation of complex interface topological motion processes, thus providing an excellent solution for engineering applications. BRIEF DESCRIPTION OF THE DRAWINGS

[0034] Figure 1 This is a flow chart of a method for predicting multiphase flow of underwater explosions according to an embodiment of the present application.

[0035] Figure 2 This is a schematic diagram of the computational domain and the initial state of the shock wave in a test case of a double-Mach reflection strong discontinuity problem.

[0036] 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 results of Titarev.

[0037] Figure 4 The present invention is a method for predicting underwater explosion multiphase flow, which is used to simulate an underwater near-free surface explosion evolution process and predict the evolution cloud diagram of fluid density, fluid pressure and moving interface.

[0038] Figure 5 The present invention uses the underwater explosion multiphase flow prediction method of the present invention to simulate the evolution process of an underwater near-wall explosion bubble water jet, and obtains an evolution cloud diagram of fluid density, fluid pressure and motion interface. DETAILED DESCRIPTION

[0039] The specific implementation of this application will be further described below with reference to the accompanying drawings.

[0040] This application discloses a method for predicting underwater explosion multiphase flow based on level-set function. Figure 1 As shown in the flowchart, the underwater explosion multiphase flow prediction method includes:

[0041] First, the computational domain is meshed. The computational domain for underwater explosion scenarios is a rectangular computational domain of predetermined size. A 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 and y directions to obtain several rectangular grids. The grids along the i-th column in the x direction and the j-th column in the y direction in the computational domain are defined as grids. . Where i and j are integer parameters.

[0042] This application uses the level-set function to predict multiphase flow and define any grid The interface distance function at the nth time step is ,in, It's a grid At any point within , the initial value of the integer parameter n is 1.

[0043] When n=1, the interface distance function of each grid at the nth time step is initialized as , then according to any grid Interface distance function at the nth time step Solve to get the motion interface of the nth time step and fluid coverage area ,include:

[0044] Determine any grid Internal point Located on the motion interface, determine point Located in the fluid coverage area, determine such that point Located in the gas coverage area. Comprehensive all points on the motion interface in all grids Get the motion interface of the nth time step , all points in all grids that are within the fluid coverage area Get the fluid coverage area at the nth time step , all points in all grids within the gas coverage area Get the gas coverage area at the nth time step , that is, the computational domain of the nth time step can be defined as:

[0045]

[0046] In addition, when determining the fluid coverage area at the nth time step After that, the fluid coverage area at the nth time step 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.

[0047] The governing equations for Euler fluids are:

[0048]

[0049] Among them, any point in the computational domain Conserved variables at , , , is the fluid density, Yes The fluid velocity along the x direction at , Yes The fluid velocity along the y direction at , is the fluid pressure, is the total energy, Indicates time.

[0050] The above Euler fluid control equation is applied to any grid Integrate above and introduce the divergence formula, In the grid There are continuous first-order partial derivatives, grid The boundary of the directed closed curve 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 scheme:

[0051]

[0052] in, It's a grid Any point inside Conserved variables at .

[0053] The grid boundary Points on The split is divided into the integral sum of each line segment on the directed closed curve, so that the above semi-discrete finite volume format is simplified to obtain the semi-discrete equation:

[0054]

[0055] in, It's a grid The border The first of the directed closed curves The numerical flux on the line segment is projected onto its unit external normal vector, the grid The border The directed closed curves of Line segments. Numerical flux Gaussian integral and numerical flux can be used The 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.

[0056] By solving the above semi-discrete equation using the third-order TVD-Runge-Kutta method, we can obtain the fluid coverage area at the nth time step: Conserved variables To obtain the multiphase flow prediction result at the n+1th time step. When solving, the above semi-discrete equation can be rewritten as:

[0057]

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

[0059]

[0060] in, is the time step. The key to the above discretization method is to solve the numerical flux in the normal direction of the grid boundary. 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.

[0061] 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, according to the motion interface of the nth time step The unit external normal vector and curvature at each point can be used to calculate the interface distance function of each grid at the next time step, but the motion interface calculated in this way The interface distance function of the nearby grid in the next time step does not strictly satisfy 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 any grid by solving the following function update formula Interface distance function at the n+1th time step :

[0062]

[0063] in, is the time step, which is the same as the time step for solving the Euler fluid governing equations It has nothing to do with the grid size, so you can choose the value based on the grid size. , It's a grid width. , is a sign function, When greater than 0, otherwise .

[0064] 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 computing resource efficiency and improve computing efficiency while ensuring solution accuracy. The time discretization of the second-order Runge-Kutta method is as follows:

[0065]

[0066] in, is an intermediate variable in the solution process. and Any parameter in have:

[0067]

[0068] Among them,

[0069]

[0070] So we can get:

[0071]

[0072] in,

[0073]

[0074] In the above formula, It's a grid width. , As a parameter, generally .

[0075]

[0076]

[0077] in,

[0078]

[0079]

[0080] In the above formula,

[0081]

[0082] 、 、 、 The same calculation can be done. For any parameter and ,function . Representation Grid The adjacent grid on the left along the x direction The interface distance function at the nth time step, Representation Grid The adjacent grid on the right along the x direction The interface distance function at the nth time step, Representation Grid The adjacent grid above along the y direction The interface distance function at the nth time step, Representation Grid The adjacent grid below along the y direction Interface distance function at the nth time step.

[0083] Through the above formula, each grid can be updated Interface distance function at the n+1th time step Then enter the calculation of the next time step and repeat the above process to make numerical forecasts.

[0084] 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 reflecting 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 walls adopt the reflecting 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 for the numerical simulation are obtained based on the normal shock wave relationship. The relationship between the shock wave before and after is as follows:

[0085]

[0086] 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 time T=0.2s, the shock wave propagates to the lower right corner of 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) in the figure, Titarev's reference simulation results show the pressure wave front distribution as follows Figure 3 As shown in (b), compared with Figure 3 As can be seen from (a) and (b) in the figure, 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.

[0087] 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. Considering the reflection of the pool bottom, it is considered as total reflection in the calculation. The underwater explosion multiphase flow prediction method of this 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 below. Figure 4 As shown by Figure 4As can be seen, part of the initial shock wave acts on the pool bottom and reflects back to the gas-liquid interface, generating a compression wave inside the bubble and a rarefaction wave. The other part propagates toward the interface and reflects back as a rarefaction wave near the free surface. This rarefaction wave in turn acts on the gas-liquid interface, transmitting a rarefaction wave into the bubble, while a compression wave is reflected back into the water. Simulating near-surface underwater explosions has been a hot topic in recent years in the field of interface capture research. This example demonstrates that the proposed multiphase flow prediction method for underwater explosions can be well applied to this type of problem.

[0088] In a simulation example of the evolution of a bubble jet caused by an underwater near-wall explosion, a solid wall boundary is set at the left boundary of the computational domain at x=0, a symmetric boundary is set at the lower boundary, and non-reflecting boundaries are set at the upper and right boundaries. A bubble of known radius is located 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 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 to predict such problems.

[0089] The above description 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 variations directly derived or imagined 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 scope of protection of the present application.

Claims

1. A method for underwater explosion multiphase flow prediction 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 along the i-th column in the x-direction and the j-th column in the y-direction in the computational domain Interface distance function at the nth time step Solve to get the motion interface of the nth time step and fluid coverage area ; Among them, n, i, j are integer parameters, It's a grid Any point within Fluid coverage area at the nth time step 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; According to the grid Interface distance function at the nth time step , using the second-order Runge-Kutta method to solve , get the grid Interface distance function at the n+1th time step And enter the forecast of the n+1th time step, where is the time step, is a sign function, .

2. The underwater explosion multiphase flow prediction method according to claim 1, characterized in that: Get the grid Interface distance function at the n+1th time step include: according to Calculate the interface distance function , among which, for and Any parameter in have: in, and, in, It's a grid The width, , As a parameter.

3. The underwater explosion multiphase flow prediction method according to claim 2, characterized in that: in, in, For any parameter and ,function .

4. The underwater explosion multiphase flow prediction method according to claim 1, characterized in that: Solve to get the motion interface of the nth time step and fluid coverage area include: Determine any grid Internal point Located on the motion interface, determine point Located in the fluid coverage area, all points on the moving interface in all grids are integrated Get the motion interface of the nth time step , all points in all grids that are within the fluid coverage area Get the fluid coverage area at the nth time step .

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: Euler fluid governing equations In any grid Integrate above and introduce the divergence formula, In the grid There are continuous first-order partial derivatives, grid The boundary of the directed closed curve 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 scheme: in, It's a grid Any point inside Conserved variables at , any point in the computational domain Conserved variables at , , , is the fluid density, Yes The fluid velocity along the x direction at , Yes The fluid velocity along the y direction at , is the fluid pressure, is the total energy, Indicates time; The grid 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: in, It's a grid The border The first of the directed closed curves The numerical flux on the line segment is projected onto its unit external normal vector, the grid The border The directed closed curves of line segments; The semi-discrete equation is solved using the third-order TVD-Runge-Kutta method to obtain the fluid coverage area at the nth time step. Conserved variables To obtain the multiphase flow prediction results at the n+1th time step.

Images