Method for controlling physical quantity of underwater explosion near naval vessel based on multiphase flow model

Abstract

The invention relates to a method for controlling physical quantity of underwater explosion near a naval vessel based on a multiphase flow model, and belongs to an anti-impact technology of underwater explosion near the naval vessel. According to the method, an additional correction equation related to total mixed energy is created, and a model is expanded into seven equation models. A mixed energy correction equation is introduced, and a more accurate gas state equation is combined, so that the difference between the predicted thermodynamic states under the shock wave condition is solved; an algorithm program of an improved six-equation model is provided on an unstructured grid system, a homogeneous hyperbolic equation is solved by using a second-order MUSCL-Hancock format based on least square reconstruction and a Barth-Jespersen limiter and a two-phase flow HLLC Riemann solver, and then an instantaneous pressure relaxation equation is solved by using a Newton-Raphson iteration method. By improving the six-equation model in the diffusion interface method, the interface capture precision and the numerical algorithm stability are improved.

Description

Technical Field

[0001] The present invention relates to underwater explosion impact resistance technology near a ship, and in particular to a full flow field physical quantity control method of an underwater explosion near a ship based on an unstructured grid compressible multiphase flow model, belonging to the field of underwater explosion numerical control. Technical Background

[0002] When a high explosive charge is detonated in water, strong underwater shock waves propagate in all directions and generate large explosion bubbles at high pressure and temperature. This can cause severe damage to any nearby structures, such as naval surface ships or submarines. Underwater explosions near surface ships are a highly complex phenomenon involving the detonation of the explosive, the interaction of the underwater shock with the free surface, the bulk cavitation generated by the rarefaction waves reflected from the free surface, the expansion and contraction of the explosive bubble and its interaction with the free surface, the water jet generated by the collapse of the bubble near the structure, the nonlinear deformation of the structure, and the interaction between the fluid and the structure. In addition, the speed of underwater shock propagation is much faster than the speed of bubble motion. This means that the underwater explosion phenomena near ships have different time scales and should be divided into two stages: (i) the propagation and interaction of the underwater shock wave initially occurs in a short time interval, and (ii) the expansion and contraction of the explosive bubble occurs in a longer time interval. Therefore, long-term calculations with very small time intervals are required to observe both phenomena simultaneously. This is a challenging task for numerical algorithms, as it requires the ability to accurately capture the wave propagation while being efficient and robust for such long-term calculations.

[0003] Over the past two decades, many researchers have used commercial software such as DYNA and USA to study the structural response and impact loads of underwater explosions near ships. Although commercial software easily provides solutions to three-dimensional complex underwater explosion problems, it is well known that it is somewhat inaccurate compared with experimental or theoretical results; for example, the inaccuracy of commercial specifications in dealing with wave interactions at moving interfaces (such as free surfaces and bubble surfaces) may lead to the emergence of pseudo-cavitation and small energy losses inside bubbles. In recent years, in order to overcome the problem of numerical inaccuracy, related research at home and abroad has developed two types of numerical methods for multiphase flow simulation. The first method is called the sharp-interface method (SIM), which treats the interface as a sharp discontinuity. The second method is called the diffuse interface method (DIM), which treats the interface as a numerical diffusion region, similar to the contact discontinuity in aerodynamics. Sharp interface methods such as volume of fluid and level set methods have attracted attention due to their simple transition condition specifications and the ability to effectively track interfaces. However, such methods are traditionally limited by low-speed flows, low density-viscosity ratios, and mass conservation errors. The diffusion interface law is mainly based on hyperbolic multiphase flow models, such as the Baer-Nunziato (BN) seven-equation model and its variants. The same algorithm can be applied to the pure single fluid and mixed multiphase flow regions. Therefore, unlike the sharp interface method, the diffusion interface method does not require special treatment of the interface. In addition, it can dynamically generate interfaces that do not exist initially, and can handle interfaces that separate pure single fluids from mixed multiphase flows.

[0004] Based on the BN seven-equation model, Saurel and Abgrall developed a single-component inviscid seven-equation model by assuming that all fluids are continuous. They used multiple sets of equilibrium equations to treat the fluid phases separately, which describe the motion state of each fluid, including two pressures, two velocities, and two temperatures. In addition to the equilibrium equations, this method also uses a transport equation to monitor the volume fraction changes of different fluids. However, since this method needs to solve a large set of equations, there are 12 equations in the case of two-dimensional two-phase flow and 14 equations in three-dimensional space, the computational cost is extremely high. To solve this problem, Kapila et al. developed a simplified two-phase flow five-equation model that uses a single pressure, a single velocity, and two temperatures, and phase equilibrium is enforced. The model demonstrates its ability to simulate wave propagation at different interfaces of compressible fluids and in compressible mixtures. However, the limitation of the five-equation model is that due to the use of a single pressure, a mixture state equation (EOS) with a non-monotonic sound speed must be established, which leads to challenges in maintaining the positivity of the volume fraction and the accuracy of the propagating waves across the interface.

[0005] Therefore, in order to solve this problem, Saurel et al. then proposed a simplified six-equation model for simulating multiphase flow problems involving interfaces, cavitation, and mixture shock waves. Saurel's six-equation model uses two pressures and only one velocity to accurately maintain the pressure non-equilibrium effect in the limit of rigid velocity relaxation. As described by Kapila et al., analytical, numerical, and permeation experiments on the BN model have demonstrated that the velocity relaxation length is very small, so it is natural to consider a simplified model with a single velocity field. This model has better characteristics than Kapila's five-equation model because it can easily maintain the positivity of the volume fraction and the monotonic mixed sound speed. However, since the original six-equation model uses a non-conservative internal energy equation, there are differences in the prediction of the thermodynamic state under shock wave conditions, which affects the calculation accuracy of cross-interface wave propagation. In addition, the six-equation model was originally used only for the simulation of simple multiphase flow problems, such as shock wave-bubble, shock wave-droplet, and shock wave-water-cylinder interactions. Although the six-equation model shows promise for simulating underwater explosions and shock interactions, it has not been fully studied according to surveys. Summary of the invention

[0006] The technical problem to be solved by the present invention is to provide a method for controlling the physical quantities of the entire flow field of an underwater explosion near a ship based on an unstructured grid compressible multiphase flow model in response to the problems and shortcomings of the prior art. The method has good algorithm accuracy under shock wave conditions and low spatial calculation cost of a structured grid system.

[0007] The object of the present invention is achieved by the following technical solutions. The present invention is a method for controlling the physical quantity of underwater explosion near a ship based on a multiphase flow model, which is characterized in that the steps are as follows:

[0008] Method Model: An additional correction equation related to the total mixing energy is created to expand the six-equation two-phase flow simplified model into the following seven-equation model, namely the volume fraction equation (1), mass conservation equations (2) and (3), momentum conservation equation (4), energy conservation equations (5) and (6), and the total mixing energy correction equation (7):

[0009]

[0010] Among them, u, α m , ρ m , p m and e m (m = 1, 2) represents velocity, phase volume fraction, phase density, phase pressure and phase internal energy; the relaxation parameter μ represents the dynamic compaction viscosity, which is used to represent the pressure work involved in the pressure relaxation process; the global variables and mixing parameters are defined by the average volume as follows:

[0011] α2=1-α1, ρ=α1ρ1+α2ρ2, p=α1p1+α2p2, Y m =(α m ρ m ) / ρ, e=Y1e1+Y2e2,

[0012] Where E is the specific total mixing energy, Y is the mass fraction, and the interface pressure p I It is calculated from the asymptotic limit of the interface pressure in the seven-equation symmetric nonequilibrium model:

[0013]

[0014] Among them, Z m =ρ m c m represents the acoustic impedance of phase m, where ρ m is the density of phase m, c m is the speed of sound of phase m.

[0015] The above-mentioned method for controlling the physical quantity of underwater explosion near a ship based on a multiphase flow model has a further preferred technical solution: the method adopts an improved rigid gas state equation SG-EOS to describe the motion characteristics of multiple liquid or gas phase fluids:

[0016]

[0017] where γm is the specific heat ratio; C v，m is the constant volume heat capacity; p ∞，m is a pressure parameter used to determine the degree of hardness relative to an ideal gas, with larger values ​​indicating high incompressibility; while e ∞，m represents the internal energy shift associated with the phase transition; subsequently, the phase velocity c m It is expressed as:

[0018]

[0019] The mixed sound speed of this model is expressed as:

[0020]

[0021] where Y m is the mass fraction of the mth phase (Y1+Y2=1);

[0022] According to the Rankine-Hugoniot condition, the flow characteristics (ρ, u, p, e) after the shock wave are expressed as:

[0023] ρ ref u s =ρ(u s -u),#(13)

[0024] pp ref =ρ ref u s u,#(14)

[0025]

[0026] After algebraic calculations of equations (13) to (15) and equations (9) and (10), the impact velocity u is obtained: s The quadratic equation is as follows:

[0027]

[0028] Where, the reference sound speed c ref Defined as:

[0029]

[0030] Finally, by considering the positive solutions in equation (16), we obtain u s and the following relationship between u:

[0031]

[0032] By fitting the shock Hugoniot data to Eq. (18) and obtaining p from Eq. (17) ∞ , we can estimate γ and cref The value of liquid water SG-EOS parameter c is obtained by this method. ref =1485m / s,γ=5.45 and p ∞ =4.045×10 8 Pa;

[0033] The shock Hugoniot curves calculated from equation (18) are compared with the experimental data; the relative error between the numerical solution and the experimental data is defined as:

[0034]

[0035] in, and They represent the shock wave velocity of the experimental data and numerical solution at a specific location respectively.

[0036] The above-mentioned method for controlling physical quantities of underwater explosions near ships based on a multiphase flow model has a further preferred technical solution: equation models (1) to (7) are solved using the following FSM algorithm:

[0037]

[0038] Here, L h represents the hyperbolic operator, L s represents the source operator;

[0039] The hyperbolic operators of the system are expressed in conservative or non-conservative form and discretized using the Godunov scheme. At the mesh interface, the numerical quantities are calculated using the two-fluid HLLC approximate Riemann solver. When dealing with the six-equation model of multiphase mixed flow, it is necessary to ensure that the interfaces between different fluids have equal pressures.

[0040] The Pelanti-Shyue transient pressure relaxation method is used as part of the source operator solver in the system;

[0041] (3) Hyperbolic Operator

[0042] The Godunov format is used to discretize the conservative equations as follows:

[0043]

[0044] Among them, A i is the area of ​​the ith computational grid, L s is the length s of the grid border, is the numerical flux at the mesh boundaries computed by the Riemann solver; this flux is defined as: F n,s =(F,G)·n s where n s =(nx ,n y ) s is the unit normal vector on the mesh boundary s;

[0045]

[0046] where u n =(u, v)·n; the superscript “^” indicates the solution of the Riemann problem on the grid boundary;

[0047] (1.1) MUSCL-Hancock method

[0048] In order to achieve second-order accuracy in time and space, the original variables W = (α1, ρ1, ρ2, u, v, p1, p2) T The MUSCL-Hancock method was used; the governing equation (1) was expressed in a non-conservative form using the original variables as follows:

[0049]

[0050] The coefficient matrix is ​​as follows:

[0051]

[0052] In order to suppress spurious oscillations near the solution gradient, a slope limiter is introduced in the reconstruction process; then the slope limiter φ i ∈[0, 1] performs piecewise linear reconstruction on the data:

[0053]

[0054] where r is the position vector; the original variable The gradient of is estimated using the least squares method:

[0055]

[0056] Among them, Ni represents the number of neighbors of grid i; the least squares weight vector is defined as follows:

[0057]

[0058] Here the coefficients are:

[0059]

[0060] The Barth-Jespersen slope limiter was adopted, which requires the piecewise linear reconstruction of the data to be bounded by the maximum and minimum values ​​on the grid and its neighbors;

[0061] First, the maximum and minimum values ​​of the grid and its neighbors are calculated as follows:

[0062]

[0063] Secondly, the slope limiter φ of each reconstruction point s is obtained i，s :

[0064]

[0065] Finally, the final value of the limiter φ i Calculated as φ i，s Minimum value of:

[0066]

[0067] Finally, the Riemann problem on each mesh interface is computed from the left and right evolution data on the corresponding interface:

[0068]

[0069] (4) Source operator

[0070] The source operator solver system uses the Pelanti-Shyue transient pressure relaxation method implemented based on saturation constraints as shown below:

[0071]

[0072] Where (αρ) m is constant during relaxation; the system can be replaced by a single equation for the unknown relaxation pressure Pr; with the help of SG-EOS, the density equation becomes:

[0073]

[0074] Next, in order to achieve the conservation of mixing energy, a mechanical balance p is applied. r =p I =p1=p2; Combining equations (34) and (35), the following equation is obtained:

[0075]

[0076] Wherein, the superscript "0" indicates the state before the relaxation step; the instantaneous pressure relaxation equation is solved by the Newton-Raphson iteration method; once the relaxation pressure Pr is obtained, the phase density and volume fraction can be calculated;

[0077] The corrected mixture pressure is defined as follows:

[0078]

[0079] After correcting for pressure, the internal energy of each phase is initialized again to its respective EOS in the next time step, i.e., e m =e m (p c , p m ).

[0080] Compared with the prior art, the present invention has the following beneficial effects:

[0081] The present invention improves the six-equation model in the diffusion interface method, introduces the mixed energy correction equation, and combines it with a more accurate gas state equation to solve the difference between the predicted thermodynamic states under shock wave conditions; an algorithm program of the improved six-equation model is proposed on an unstructured grid system, and the second-order MUSCL-Hancock format based on least squares reconstruction and Barth-Jespersen limiter and the two-phase flow HLLC Riemann solver are used to solve the homogeneous hyperbolic equation, and then the Newton-Raphson iteration method is used to solve the instantaneous pressure relaxation equation. Then, the above model is numerically tested on this basis, and finally it is applied to the research field of underwater explosion near ships.

[0082] 1. Efficient and stable numerical algorithm development:

[0083] Explore efficient and stable numerical algorithms for underwater explosion multiphase flow problems near ships, solve underwater explosion multiphase flow problems in complex environments, and complete the solver code writing. Optimize the calculation process through parallel computing technology (Open Multiprocessing+Message Passing Interface, OpenMP+MPI), establish a new compressible multiphase flow model to improve calculation efficiency and stability.

[0084] 2. Additional correction equation associated with mixing energy:

[0085] A new additional correction equation associated with the mixing energy is introduced for Saurel's six-equation model, combined with a more accurate rigid gas equation (SG-EOS) to resolve the difference between the predicted thermodynamic states under shock wave conditions. By comparing the shock wave Hugoniot curve, more accurate SG-EOS parameters are determined to achieve a higher precision fit to the experimental results.

[0086] 3. Numerical algorithms on unstructured grid systems:

[0087] The model is extended to simulate multidimensional problems related to underwater explosions near ships in complex environments, and a numerical algorithm program for the improved six-equation model is developed on an unstructured grid system. The second-order MUSCL-Hancock format based on least squares reconstruction and Barth-Jespersen limiter is used for discretization, and the two-phase flow HLLC Riemann solver is used to solve the homogeneous hyperbolic equations to address the shortcomings in multiple dimensions.

[0088] 4. Interface capture and calculation in complex geometric environments:

[0089] The Immersed Boundary Method (IBM) is combined to expand the computational domain to any complex terrain to ensure accurate capture of the fluid interface in a complex geometric environment. The accuracy of interface capture is improved by improving the six-equation model in the diffusion interface method, and the Newton-Raphson iteration method is used to solve the instantaneous pressure relaxation equation to ensure the stability of the calculation process.

[0090] (1) Accurate simulation of multiphase flow phenomena of underwater explosions near ships:

[0091] Traditional multiphase flow models are not accurate and reliable enough when dealing with complex flow phenomena (such as explosion shock waves, bubble expansion and collapse, etc.). The present invention proposes an improved six-equation model and combines it with a more accurate gas state equation to improve the simulation accuracy of these complex multiphase flow phenomena.

[0092] (2) Application of unstructured grids in complex geometric environments:

[0093] The fluid dynamics problem of underwater explosion near a ship has complex geometric boundary conditions, which are difficult to adapt to with traditional structured grids. This invention develops a numerical algorithm suitable for unstructured grids, enabling the model to perform accurate calculations in a complex geometric environment, thereby more realistically simulating the effects of explosions around ships.

[0094] (3) Capturing and numerical stability of multiphase flow interfaces:

[0095] In underwater explosion simulation, accurately capturing fluid interfaces (such as air-water interfaces) is one of the key issues. Traditional methods have difficulties in dealing with fluid interfaces, which easily leads to numerical instability. The present invention improves the accuracy of interface capture and the stability of the numerical algorithm by improving the six-equation model in the diffusion interface method, adopting the second-order MUSCL-Hancock format based on least squares reconstruction and Barth-Jespersen limiter and the two-phase flow HLLC Riemann solver.

[0096] (1) Improved accuracy of underwater explosion simulation: The improved six-equation model proposed in the present invention combined with a more accurate gas state equation can more accurately simulate multiphase flow phenomena of underwater explosions near ships, such as explosion shock waves, bubble expansion and collapse, etc.

[0097] (2) Applicable to complex geometric environments: The application of this technology on unstructured grid systems enables it to perform accurate calculations in complex geometric environments, thereby more realistically simulating the effects of explosions around ships.

[0098] (3) Improving the accuracy and numerical stability of interface capture: The accuracy of interface capture and the stability of the numerical algorithm are improved by improving the six-equation model in the diffusion interface method and adopting the second-order MUSCL-Hancock format and two-phase flow HLLC Riemann solver based on least squares reconstruction and Barth-Jespersen limiter.

[0099] (4) Breaking the foreign technology monopoly and promoting independent technology: This invention not only promotes the research on underwater explosion problems near ships, the development of related technologies such as hull structure design, explosion impact damage prediction and personnel safety protection, but also breaks the foreign technology monopoly in computational fluid dynamics (CFD) commercial software and promotes the widespread application of independent technology in this field. BRIEF DESCRIPTION OF THE DRAWINGS

[0100] Figure 1 A flowchart for accurately solving the physical quantities of the entire flow field of underwater explosions near ships based on the FSM algorithm of the unstructured grid system;

[0101] Figure 2 The figure is a schematic diagram of the computational unit and its boundary of a quadrilateral unstructured grid;

[0102] Figure 3 The distribution diagram of the flow field mixing density, velocity and liquid volume fraction at T = 0.24ms for the one-dimensional water-gas two-phase shock tube problem test;

[0103] Figure 4 The present invention simulates the shock wave pressure time history curve at two measurement targets at a distance of 0.0 cm and 5.0 cm from the rigid bottom surface below the explosion center;

[0104] Figure 5 The gas phase volume fraction distribution diagram of the flow field at different times for simulating the underwater explosion problem near the water surface in the present invention;

[0105] Figure 6 The flow field pressure distribution diagram at different times for simulating the underwater explosion problem near the water surface of the present invention;

[0106] Figure 7 A calculation grid distribution diagram and a mathematical model configuration diagram for simulating the surroundings of a ship in the present invention;

[0107] Figure 8 The pressure distribution diagram of the underwater explosion flow field near the ship in the first stage of the present invention simulation;

[0108] Fig. 9 The present invention simulates the shock wave pressure time history curve at two measurement targets at a straight line distance of 2 and 2.5 m along the explosion center;

[0109] Fig.10 This is the gas phase volume fraction distribution diagram of the underwater explosion flow field near the ship in the second stage of simulation of the present invention. DETAILED DESCRIPTION

[0110] The implementation process of the present invention is further described in detail below in conjunction with the accompanying drawings and specific implementation methods.

[0111] Example 1, a control method for physical quantities of underwater explosion near a ship based on a multiphase flow model, the flow chart is as follows Figure 1 :

[0112] 1 Computational model

[0113] 1.1 Governing equations

[0114] In order to compensate for the deviation of the original six-equation two-phase flow model in predicting the thermodynamic state under the shock wave environment, the present invention innovatively introduces an additional correction equation that is closely related to the overall mixing energy. This correction mechanism shows significant effectiveness on both sides of the interface approaching the single-phase flow limit. Based on this, the six-equation simplified model is expanded into the following seven-equation model, namely, the volume fraction equation (1), the mass conservation equations (2) and (3), the momentum conservation equation (4), the energy conservation equation

[0115] (5) and (6), and the mixed total energy correction equation (7):

[0116]

[0117] Among them, u, α m , ρ m , p m and e m (m=1, 2) represents velocity, phase volume fraction, phase density, phase pressure and phase internal energy; the relaxation parameter μ represents the dynamic compaction viscosity, which is used to represent the pressure work involved in the pressure relaxation process;

[0118] Global variables and blending parameters are defined by averaging the volumes as follows:

[0119] α2=1-α1, ρ=α1ρ1+α2ρ2, p=α1p1+α2p2, Y m =(α m ρm ) / ρ, e=Y1e1+Y2e2,

[0120] Where E is the specific total mixing energy, Y is the mass fraction, and the interface pressure p I It is calculated from the asymptotic limit of the interface pressure in the seven-equation symmetric nonequilibrium model:

[0121]

[0122] Among them, Z m =ρ m c m represents the acoustic impedance of phase m, where ρ m is the density of phase m, c m is the speed of sound of phase m.

[0123] 1.2 Equation of state

[0124] This method uses the improved rigid gas state equation SG-EOS to describe the motion characteristics of various liquid or gas phase fluids:

[0125]

[0126] where γ m is the specific heat ratio; C v，m is the constant volume heat capacity; p ∞，m is a pressure parameter used to determine the degree of hardness relative to an ideal gas, with larger values ​​indicating high incompressibility; while e ∞，m represents the internal energy shift associated with the phase transition; subsequently, the phase velocity c m It is expressed as:

[0127]

[0128] The mixed sound speed of this model is expressed as:

[0129]

[0130] where Y m is the mass fraction of the mth phase (Y1+Y2=1);

[0131] In equation (11), the monotonicity of the mass fraction is particularly important, which is in sharp contrast to the non-monotonic nature of the Wood speed of sound in the Kapila five-equation model. It is precisely because the six-equation model has this unique monotonicity that it has shown obvious advantages in the numerical approximation process. In the invention, the inventor proposed an innovative method to optimize the parameters of the SG-EOS state equation: by using the impact Hugoniot curve, a more accurate fit of the experimental results is achieved. According to the Rankine-Hugoniot condition, the flow characteristics (p, u, p, e) after the shock wave can be expressed as:

[0132] ρ ref u s =ρ(u s -u),#(13)

[0133] pp ref =ρ ref u s u,#(14)

[0134]

[0135] After algebraic calculations of equations (13) to (15) and equations (9) and (10), the impact velocity u is obtained: s The quadratic equation is as follows:

[0136]

[0137] Where, the reference sound speed c ref Defined as:

[0138]

[0139] Finally, by considering the positive solutions in equation (16), we obtain u s and the following relationship between u:

[0140]

[0141] By fitting the shock Hugoniot data to Eq. (18) and obtaining p from Eq. (17) ∞ , we can estimate γ and c ref The value of liquid water SG-EOS parameter c is obtained by this method. ref =1485m / s,γ=5.45 and p ∞ =4.045×10 8 Pa;

[0142] The shock Hugoniot curves calculated from equation (18) are compared with the experimental data; the relative error between the numerical solution and the experimental data is defined as:

[0143]

[0144] in, and They represent the shock wave velocity of the experimental data and numerical solution at a specific location respectively.

[0145] 2 Numerical methods

[0146] The control equations (1) to (7) are solved using the following FSM algorithm:

[0147]

[0148] Here, L h represents the hyperbolic operator, L s represents the source operator;

[0149] The hyperbolic operators of the system can be expressed in conservative or non-conservative form. The hyperbolic operators of the system can be expressed in conservative or non-conservative form. Here, the Godunov scheme is used for discretization. At the grid interface, the numerical quantities are calculated using the two-fluid HLLC approximate Riemann solver. When dealing with the six-equation model of multiphase mixed flow, it is necessary to ensure that the interfaces between different fluids have equal pressures.

[0150] The Pelanti-Shyue transient pressure relaxation method is used as part of the source operator solver in the system.

[0151] 2.1 Hyperbolic Operators

[0152] On a 2D quadrilateral unstructured grid, such as Figure 2 As shown, the Godunov format is used to discretize the conservation equations as follows:

[0153]

[0154] Among them, A i is the area of ​​the ith computational grid, L s is the length s of the grid border, is the numerical flux at the mesh boundaries computed by the Riemann solver; this flux is defined as: F n，s =(F, G)·n s where n s =(n x , n y ) s is the unit normal vector on the mesh boundary s;

[0155]

[0156] where u n =(u, v)·n; the superscript “^” indicates the solution of the Riemann problem on the grid boundary;

[0157] (1.1) MUSCL-Hancock method

[0158] In order to achieve second-order accuracy in time and space, the original variables W = (α1, ρ1, ρ2, u, v, p1, p2) T The MUSCL-Hancock method was used; the governing equation (1) was expressed in a non-conservative form using the original variables as follows:

[0159]

[0160] The coefficient matrix is ​​as follows:

[0161]

[0162] In order to suppress spurious oscillations near the solution gradient, a slope limiter is introduced in the reconstruction process; then the slope limiter φ i ∈[0, 1] performs piecewise linear reconstruction on the data:

[0163]

[0164] where r is the position vector; the original variable The gradient of is estimated using the least squares method:

[0165]

[0166] Among them, Ni represents the number of neighbors of grid i; the least squares weight vector is defined as follows:

[0167]

[0168] Here the coefficients are:

[0169]

[0170] The Barth-Jespersen slope limiter was adopted, which requires the piecewise linear reconstruction of the data to be bounded by the maximum and minimum values ​​on the grid and its neighbors;

[0171] First, the maximum and minimum values ​​of the grid and its neighbors are calculated as follows:

[0172]

[0173] Secondly, the slope limiter φ of each reconstruction point s is obtainedi，s :

[0174]

[0175] Finally, the final value of the limiter φ i Calculated as φ i，s Minimum value of:

[0176]

[0177] Finally, the Riemann problem on each mesh interface is computed from the left and right evolution data on the corresponding interface:

[0178]

[0179] (4) Source operator

[0180] The source operator solver system uses the Pelanti-Shyue transient pressure relaxation method implemented based on saturation constraints as shown below:

[0181]

[0182] Where (aρ) m is constant during relaxation; the system can be replaced by a single equation for the unknown relaxation pressure Pr; with the help of SG-EOS, the density equation becomes:

[0183]

[0184] Next, in order to achieve the conservation of mixing energy, a mechanical balance p is applied. r =p I =p1=p2; Combining equations (34) and (35), the following equation is obtained:

[0185]

[0186] Wherein, the superscript "0" indicates the state before the relaxation step; the instantaneous pressure relaxation equation is solved by the Newton-Raphson iteration method; once the relaxation pressure Pr is obtained, the phase density and volume fraction can be calculated;

[0187] The corrected mixture pressure is defined as follows:

[0188]

[0189] After correcting for pressure, the internal energy of each phase is initialized again to its respective EOS in the next time step, i.e., e m =e m (p c , p m ).

[0190] 3 Numerical tests and examples

[0191] 3.1 Water-gas two-phase shock tube problem test

[0192] In order to verify the ability of the numerical scheme of the present invention to accurately capture the shock wave and interface position, the test considered a 1m long tube filled with high pressure liquid water on the left side and air on the right side. Using the above proposed air and liquid water SG-EOS parameters Y1 = 1.4, p ∞,1 =0Pa and Y2=5.45, p ∞,2 =4.045×10 8 Pa. Table 1 lists the initial conditions for the left and right sides. Figure 3 The distribution curves of the mixed density, velocity and liquid volume fraction of the flow field at T = 0.24ms are shown. This case uses 1000 grids for testing, and the CFL number is equal to 0.7. The test results show that after correction of the mixed total energy equation, the results obtained by the model of the present invention are very consistent with the exact solution of the Euler equation. The shock wave velocity and interface between the fluids are accurate, and there is no oscillation near the interface.

[0193] Table 1 Left and right initial conditions for one-dimensional water-gas two-phase shock tube problem

[0194] Location <![CDATA[α1]]> <![CDATA[ρ1(kg / m 3 )]]> <![CDATA[ρ2(kg / m 3 )]]> u(m / s) <![CDATA[p1(MPa)]]> <![CDATA[p2(MPa)]]> left <![CDATA[10 -6 ]]> 1 <![CDATA[10 3 ]]> 0 <![CDATA[10 3 ]]> <![CDATA[10 3 ]]> right <![CDATA[1-10 -6 ]]> 1 <![CDATA[10 3 ]]> 0 0.1 0.1

[0195] 3.2 Testing of cavitation phenomena in near-water explosions

[0196] In order to further verify the ability of the numerical scheme of the present invention to simulate the cavitation phenomenon of multiphase flow, the inventors considered a case of underwater explosion near the water surface, which has been studied before. In the calculation area of ​​[0m, 4m]×[0m, 2.5m], the explosion bubble has a radius of 0.12m, is located at 0.15m underwater, and the air-water interface is located at 1.35m. The bottom boundary condition is a solid-wall boundary condition, while the other boundary conditions are transmissive boundaries with zero gradient. The number of quadrilateral unstructured grids is equal to 398588, and the CFL number is equal to 0.4. Table 3 lists the initial conditions of water, air and explosion bubbles. Before starting the simulation, in order to verify the accuracy of the numerical scheme of the present invention in two-dimensional space, the inventors additionally considered the verification case of a two-dimensional underwater explosion near a rigid interface proposed in the prior art. In this verification case, except that the pressure of the bubble is 829MPa, the other initial condition parameters are the same as those in Table 2. The inventors set two measurement targets at 0.0cm and 5.0cm from the rigid bottom surface below the explosion center. Figure 4The shock wave pressure time history curve shown shows that the calculated numerical results are in good agreement with the reference data. The shock wave pressure in the water decays rapidly with time and distance from the explosion center. The shock wave then interacts with the bottom surface, with maximum pressure pulse peaks of approximately 0.5 GPa and 0.45 GPa, respectively. The second pressure peak in the figure corresponds to the cavitation bubble collapse region. Figure 5 The contour diagram of the air volume fraction of the explosion problem near the free water surface at different times of T = 0.2, 0.4, and 0.8 ms is shown. It can be clearly observed that the explosion gas gradually expands and drives the air-water interface to form a low-pressure area between the water surface and the cavitation phenomenon evolution. Figure 6 The simulation results of the pressure contours at the corresponding time step are shown. The shock wave reflected from the bottom, the rarefaction wave reflected from the air-water interface, and the weak shock wave transmitted to the air can be clearly seen. The solid line in the figure represents the air-water interface, and the dashed line represents the bubble shape. In addition, the refraction of the shock wave when it contacts the bottom surface can be seen, forming a Mach effect.

[0197] Table 2 Initial conditions of water, air and explosion bubbles for underwater explosion problem near water surface

[0198]

[0199] 3.3 Underwater explosion simulation near ships

[0200] Finally, the inventors focused on the simulation research of underwater explosions around ships. This research is one of a series of important topics to explore the damage caused by underwater explosions to naval surface ships. At the current stage of research, the inventors temporarily regard ships as rigid and fixed structures. Although it is known that underwater impact loads are sufficient to cause deformation of the hull and even more serious damage, in order to accurately simulate these complex phenomena, it is necessary to construct a coupling model that can reflect the interaction between fluids and solids, which is also the focus of the inventors' subsequent research. The core of the present invention is to comprehensively simulate the complete time history of underwater explosions near ships, covering the entire process from the initial propagation of shock waves and their interaction with multiple interfaces (such as bubbles, water surface, and hull structure), to the subsequent expansion and contraction of bubbles, until the final violent impact of the water jet on the hull.

[0201] Figure 7The distribution diagram of the computational grid around the ship and the configuration diagram of the mathematical model are shown. The ship is 10m wide, 10m deep, and has a draft of 3m; a high-pressure and high-density explosion bubble with a radius of 1m is placed 3m to the left of the center of the ship 6m below the water surface. In the computational domain of [-30m, 30m]×[-40m, 20m], all boundaries are transmissive boundary conditions. The number of quadrilateral unstructured grids is equal to 877760, and the CFL number is equal to 0.4. Except that the bubble pressure is 2MPa, the initial conditions and constant parameters are the same as the test problem of near-surface explosion cavitation phenomenon described in Section 3.2. Figure 8 The pressure distribution diagram of the flow field near the ship at different times in the first stage of the underwater explosion (T = 0 ~ 20ms) is shown, and the underwater shock wave propagation and its interaction with the water surface, hull surface and bubbles at the beginning of the explosion: At the moment the underwater explosion begins, a high-intensity circular shock wave quickly spreads in the water. At about T = 2ms, this underwater shock wave comes into contact with the hull structure, triggering the shock wave reflection phenomenon. At about 4ms, the reflected shock wave meets the bubble and interweaves into a complex waveform area. When the main underwater shock wave hits the interface between air and water (marked by the white solid line), it will generate strong reflection rarefaction waves, which then return to the water, and at the same time, the free water surface begins to rise. Fig. 9 The following is a time history curve of the shock wave pressure at two measurement targets at a straight line distance of 2 and 2.5 m from the explosion center. The pressure rapidly increases to a peak value of about 2 MPa, and then drops to about 1.1 MPa. As the reflected shock wave contacts the bubble interface again, the flow field pressure value at the measurement target returns to about 1.4 MPa, and finally gradually decreases to the ambient pressure level of the surrounding waters. Fig.10 The distribution diagram of the gas phase volume fraction in the flow field near the ship at different times in the second stage of the underwater explosion (T = 20-1000ms) is shown. After the explosion, the bubble expands and contracts, and then interacts with the air-water interface and the hull structure: before T = 20ms, the size of the bubble is almost the same as the initial size. At T = 100ms and 200ms, it is clearly seen that the bubble expands and the air-water interface begins to rise. At T = 600ms, the expanded bubble hits the surface of the hull, and the air-water interface continues to rise. At T = 900ms, the over-expanded bubble begins to shrink, and the rising height of the air-water interface exceeds the height of the ship. Over time, the bubble shrinks further, which will cause the water jet to hit the ship structure.

[0202] 4 Conclusion

[0203] The present invention uses an improved six-equation compressible multiphase flow model to numerically simulate the underwater explosion problem near a ship. The model uses a second-order MUSCL-Hancock format based on least squares reconstruction and Barth-Jespersen limiter and a two-phase flow HLLC Riemann solver to solve homogeneous hyperbolic equations on a quadrilateral unstructured grid, and then uses the Newton-Raphson iteration method to solve the instantaneous pressure relaxation equation. The calculation model is tested on the water-gas two-phase shock tube problem and the cavitation phenomenon of near-surface explosion, which shows that it has a good capture ability for shock waves, interface positions and cavitation phenomena. Finally, the present invention considers the underwater explosion problem near a ship, and places a high-pressure and high-density bubble with an initial explosion pressure of 2Mpa on the left side of the hull below the water surface. Through calculation, the two main physical characteristics of the underwater explosion are accurately simulated: the propagation and interaction of the explosion shock wave in a short time interval (0-20ms) in the first stage; the water jet effect when the bubble expands and collapses in the second stage of oscillation (20-1000ms) in a longer time interval. The inventors found that in the early stage of underwater explosion, the impact pressure of the shock wave on the ship structure is very large, but its duration is very short. The duration of the water jet pressure in the later stage is much longer than the initial explosion impact. The calculation method and research results have reference and guiding value for improving the understanding of the anti-shock mechanism of underwater explosion near ships and its in-depth research.

Claims

1. A method for controlling physical quantities of underwater explosions near a ship based on a multiphase flow model, characterized in that: The steps are as follows: Method Model: An additional correction equation related to the total mixing energy is created to expand the six-equation two-phase flow simplified model into the following seven-equation model, namely the volume fraction equation (1), mass conservation equations (2) and (3), momentum conservation equation (4), energy conservation equations (5) and (6), and the total mixing energy correction equation (7): Among them, u, α m , ρ m , p m and e m (m=1, 2) represents velocity, phase volume fraction, phase density, phase pressure and phase internal energy; the relaxation parameter μ represents the dynamic compaction viscosity, which is used to represent the pressure work involved in the pressure relaxation process; Global variables and blending parameters are defined by averaging the volumes as follows: α2=1-α1,ρ=α1ρ1+α2ρ2,p=α1p1+α2p2,Y m =(a m r m ) / ρ, e=Y1e1+Y2e2, Where E is the specific total mixing energy, Y is the mass fraction, and the interface pressure p I It is calculated from the asymptotic limit of the interface pressure in the seven-equation symmetric nonequilibrium model: Among them, Z m =ρ m c m represents the acoustic impedance of phase m, where ρ m is the density of phase m, c m is the speed of sound of phase m.

2. The method for controlling physical quantities of underwater explosions near ships based on a multiphase flow model according to claim 1 is characterized in that: This method uses the improved rigid gas state equation SG-EOS to describe the motion characteristics of various liquid or gas phase fluids: where γ m is the specific heat ratio; C v，m is the constant volume heat capacity; p ∞，m is a pressure parameter used to determine the degree of hardness relative to an ideal gas, with larger values ​​indicating high incompressibility; while e ∞，m represents the internal energy shift associated with the phase transition; subsequently, the phase velocity c m It is expressed as: The mixed sound speed of this model is expressed as: where Y m is the mass fraction of the mth phase (Y1+Y2=1); According to the Rankine-Hugoniot condition, the flow characteristics (ρ, u, p, e) after the shock wave are expressed as: ρ ref in s =ρ(in s -u),#(13) pp ref =ρ ref in s in, #(14) After algebraic calculations of equations (13) to (15) and equations (9) and (10), the impact velocity u is obtained: s The quadratic equation is as follows: Where, the reference sound speed c ref Defined as: Finally, by considering the positive solutions in equation (16), we obtain u s and the following relationship between u: By fitting the shock Hugoniot data to Eq. (18) and obtaining p from Eq. (17) ∞ , we can estimate γ and c ref The value of liquid water SG-EOS parameter c is obtained by this method. ref =1485m / s,γ=5.45 and p ∞ =4.045×10 8 Pa; The shock Hugoniot curves calculated from equation (18) are compared with the experimental data; the relative error between the numerical solution and the experimental data is defined as: in, and They represent the shock wave velocity of the experimental data and numerical solution at a specific location respectively.

3. The method for controlling physical quantities of underwater explosions near ships based on a multiphase flow model according to claim 1, characterized in that: The control equations (1) to (7) are solved using the following FSM algorithm: Here, L h represents the hyperbolic operator, L s represents the source operator; The hyperbolic operators of the system are expressed in conservative or non-conservative form and discretized using the Godunov scheme. At the mesh interface, the numerical quantities are calculated using the two-fluid HLLC approximate Riemann solver. When dealing with the six-equation model of multiphase mixed flow, it is necessary to ensure that the interfaces between different fluids have equal pressures. The Pelanti-Shyue transient pressure relaxation method is used as part of the source operator solver in the system; (1) Hyperbolic Operator The Godunov format is used to discretize the conservative equations as follows: Among them, A i is the area of ​​the ith computational grid, L s is the length s of the grid border, is the numerical flux at the mesh boundaries computed by the Riemann solver; this flux is defined as: F n，s =(F, G)·n s Where n s =(n x , n y ) s is the unit normal vector on the mesh boundary s; where u n =(u, v)·n; the superscript "∧" indicates the solution of the Riemann problem on the grid boundary; (1.1) MUSCL-Hancock method In order to achieve second-order accuracy in time and space, the original variables W = (α1, ρ1, ρ2, u, v, p1, p2) T The MUSCL-Hancock method was used; the governing equation (1) was expressed in a non-conservative form using the original variables as follows: The coefficient matrix is ​​as follows: In order to suppress spurious oscillations near the solution gradient, a slope limiter is introduced in the reconstruction process; then the slope limiter φ i ∈[0, 1] performs piecewise linear reconstruction on the data: where r is the position vector; the original variable The gradient of is estimated using the least squares method: Among them, Ni represents the number of neighbors of grid i; the least squares weight vector is defined as follows: Here the coefficients are: The Barth-Jespersen slope limiter was adopted, which requires the piecewise linear reconstruction of the data to be bounded by the maximum and minimum values ​​on the grid and its neighbors; First, the maximum and minimum values ​​of the grid and its neighbors are calculated as follows: Secondly, the slope limiter φ of each reconstruction point s is obtained i，s : Finally, the final value of the limiter φi is calculated as the minimum value of φis: Finally, the Riemann problem on each mesh interface is computed from the left and right evolution data on the corresponding interface: (2) Source operator The source operator solver system uses the Pelanti-Shyue transient pressure relaxation method implemented based on saturation constraints as shown below: Where (αρ) m is constant during relaxation; the system can be replaced by a single equation for the unknown relaxation pressure Pr; with the help of SG-EOS, the density equation becomes: Next, in order to achieve the conservation of mixing energy, a mechanical balance p is applied. r =p l =p1=p2; Combining equations (34) and (35), the following equation is obtained: Wherein, the superscript "0" indicates the state before the relaxation step; the instantaneous pressure relaxation equation is solved by the Newton-Raphson iteration method; once the relaxation pressure Pr is obtained, the phase density and volume fraction can be calculated; The corrected mixture pressure is defined as follows: After correcting for pressure, the internal energy of each phase is initialized again to its respective EOS in the next time step, i.e., e m =e m (p c , p m ).

Images