A computational method for propagation of underwater explosion shock waves in non-uniform acoustic flow fields

Abstract

The present invention provides a method for calculating the propagation of underwater explosion shock waves in a non-uniform acoustic flow field. The method uses an axisymmetric calculation model of underwater explosion loads in a non-homogeneous fluid based on the local discontinuous Galerkin method and continuity conditions to simulate the propagation of underwater explosion shock waves in a non-homogeneous fluid. By constructing a calculation model specifically for non-homogeneous fluids and combining continuity conditions, an effective simulation of the propagation process of underwater explosion shock waves in a non-homogeneous flow field containing a discontinuous sound velocity interface and a sound velocity gradient is achieved. The present invention can accurately analyze the influence of the sound velocity jump layer interface on the propagation path and load peak of the underwater explosion shock wave, and provide more reliable underwater explosion load data for submarines operating near the sound velocity jump layer. The present invention can effectively study the influence of the flow field sound velocity gradient on the underwater explosion load and cavitation characteristics, and optimize the design and construction plan of protective ships and marine structures based on more accurate underwater explosion load input data.

Description

Technical Field

[0001] The invention belongs to the technical field of explosion shock wave calculation based on computer data processing, and in particular relates to a calculation method for propagation of underwater explosion shock waves in a non-uniform acoustic flow field. Background Art

[0002] The flow field background of existing numerical calculation models for underwater explosions is mostly homogeneous fluid, and the influence of fluid heterogeneity on underwater explosion loads is not considered. However, sonic jump layers are common in real marine environments, and submarines often operate in this area. The change in sound speed within the sonic jump layer makes it difficult for sonar to detect the position of submarines, and the shock wave of underwater explosions may pass through the sonic jump layer area. In previous studies, the analysis of underwater explosions is often based on the idealized homogeneous flow field assumption, ignoring the complex sound speed distribution in the ocean. This simplified treatment leads to deviations between the calculation results and the actual situation in many application scenarios related to underwater explosions, such as underwater engineering protection. For example, it is impossible to accurately estimate the changes in the transmission and effect of shock loads caused by sonic jump layers when underwater explosions attack submarines, which affects the accuracy of tactical decisions; in marine engineering construction, because the sound speed heterogeneity is not considered, the reliability and safety of the protection design of facilities that may be threatened by underwater explosions cannot be fully guaranteed. Therefore, it is necessary to study the influence of the non-uniformity of flow field sound speed on the shock load and cavitation characteristics of underwater explosions to fill the gap in the existing technology in this regard. Summary of the invention

[0003] In response to the above problems, the invention aims to solve the technical problem that the existing numerical calculation models for underwater explosions are mostly based on a homogeneous fluid background and do not fully consider the impact of fluid heterogeneity on underwater explosion loads, resulting in deviations between the calculation results and the actual situation in many practical application scenarios involving marine engineering protection.

[0004] The present invention provides a method for calculating the propagation of underwater explosion shock waves in a non-uniform acoustic flow field, which is characterized by comprising the following processes:

[0005] Step 1, establish the flow field range and size, and discretize the flow field into a rectangular grid; based on the corresponding rectangular grid node and unit information data, split the rectangular unit and number the new node and unit to obtain the triangular unit flow field data containing the number information;

[0006] Step 2: Based on the triangular unit flow field data, the incident wave loading surface, explosion point and source point positions are set, and then the Geers-Hunter model and the far-field underwater explosion pressure calculation formula are used to calculate the time-load data at the explosion point and the dynamic pressure data on the discrete triangular unit flow field incident wave loading surface;

[0007] Step 3, using the LDG method to derive the second-order wave equation for the propagation of underwater explosion shock waves in discrete triangular unit flow fields, and obtain the control equation for calculating the propagation of underwater explosion shock waves in discrete triangular unit flow fields;

[0008] Step 4, based on the obtained control equation, the calculation of each integral term in the control equation is converted into a Gaussian integral formula calculation through coordinate transformation, and the sound velocity value at each Gaussian integral point is calculated;

[0009] Step 5, calculating the numerical flux of the discrete triangular unit flow field where the sound velocity is continuous, discontinuous and at the unit boundary of the flow field;

[0010] Step 6: Calculate the approximate solutions of auxiliary variables 1 and 2 in the control equation. Auxiliary variable 1 is the dynamic pressure of the fluid on the space. r Take the derivative of the direction and multiply by the speed of sound c , auxiliary variable 2 is the dynamic pressure of the fluid on the space z Take the derivative of the direction and multiply by the speed of sound c ;

[0011] Step 7, calculating the second derivative of the flow field dynamic pressure with respect to time;

[0012] Step 8, calculating the dynamic pressure results of the discrete triangular unit flow field, introducing a pressure truncation model to deal with the cavitation effect during the calculation, and adding artificial volume viscosity to the calculated dynamic pressure results;

[0013] Step 9, based on the calculation results of steps 3 to 8, obtain the physical quantities of all units of the discrete triangular unit flow field at the initial moment, and use the fourth-order Runge-Kutta method to calculate the physical quantities of all units of the discrete triangular unit flow field at the next moment, so as to simulate the propagation characteristics of the underwater explosion shock wave in the inhomogeneous flow field.

[0014] Preferably, in step 1, first, the flow field range and flow field size are established in ABAQUS finite element software, and the entire flow field is discretized into rectangular grids using ABAQUS finite element software, each rectangular grid is a unit; each rectangular unit is split into 4 triangular units, and the newly generated nodes and units are numbered and Ω Represents the discretized triangular unit flow field.

[0015] Preferably, in step 2, the positions of the discrete triangular unit flow field explosion point, source point and incident wave loading surface are set, and for the shock wave generated by the underwater explosion of explosives, the Geers-Hunter model is used to calculate the time-load data of the explosive explosion at the discrete triangular unit flow field explosion point. p t After the calculation is completed, the time-load data at the explosion point p tImport the program, located in the discrete triangular unit flow field Ω Any point x on the incident wave loading surface j The calculation formula for the gradient of the dynamic pressure at is:

[0016] (1)

[0017] (2)

[0018] (3)

[0019] where represents the gradient of the dynamic pressure on the incident wave loading surface, and represents the vector position coordinates of the explosive source point and the explosion point, respectively, represents the vector position coordinates of the point on the incident wave loading surface, R 0 and R j They represent the distance between the explosive source and the explosion point and the distance between any point on the incident wave loading surface and the explosive source, respectively. for p t The derivative with respect to time.

[0020] Preferably, the step 3 is specifically:

[0021] The propagation of underwater explosion shock waves in fluid satisfies the second-order wave equation, and the calculation formula of the second-order wave equation is:

[0022] (4)

[0023] In the formula p is the dynamic pressure of the fluid, Δ represents the Laplace operator, represents the second derivative of the fluid's dynamic pressure with respect to time, c The fluid in the flow field Ω The speed of sound inside;

[0024] The three-dimensional flow field model is simplified into an axisymmetric model. In the axisymmetric model in the cylindrical coordinate system, formula (4) is rewritten as:

[0025] (5)

[0026] In the formula r and z are the radial coordinate and the axial coordinate, c ( r , z ) represents the flow field in the axisymmetric model ΩThe speed of sound at any point in The dynamic pressure of the fluid r Direction partial derivative, The dynamic pressure of the fluid r Direction to find the second-order partial derivative, The dynamic pressure of the fluid z Direction to find the second-order partial derivative;

[0027] Auxiliary variables Introduced into the axisymmetric model in the cylindrical coordinate system, the axisymmetric model of the second-order wave equation in the cylindrical coordinate system is expressed as:

[0028] (6)

[0029] (7)

[0030] Represents the derivative operation symbol, The dynamic pressure of the fluid on the space r Take the partial derivative of the direction and multiply it by the speed of sound c ; The dynamic pressure of the fluid on the space z Take the partial derivative of the direction and multiply it by the speed of sound c ;

[0031] The second-order wave equation of the above axisymmetric model is multiplied by the test function ψ , ϕ and φ , and integrating the discrete units, we get:

[0032] (8)

[0033] (9)

[0034] (10)

[0035] In the formula K Representing flow field Ω Discrete triangular elements in ; Represents the product of auxiliary variable 1 and the speed of sound r Direction partial derivative; The product of auxiliary variable 2 and the speed of sound is z Direction partial derivative; and Respectively represent the fluid dynamic pressure r and z Find partial derivatives;

[0036] The control equation for the propagation of underwater explosion shock waves in the discrete triangular unit flow field is obtained by the partial integration operation rule and Gaussian formula. The control equation is expressed as:

[0037] (11)

[0038] (12)

[0039] (13)

[0040] In the formula Represents the fluid dynamic pressure in discrete triangular elements K The approximate solution of and Respectively represent the auxiliary variables in the discrete triangle unit K The approximate solution of ∂K represents the boundary of the discrete triangle unit, is the unit external normal vector to the cell boundary, and Respectively represent the test function y right r and z Find the partial derivative, Represents the test function f right r Find the partial derivative, Represents the test function j right z Find the partial derivative, and They represent the flow field sound velocity r and z Directional derivative, , and is the numerical flux at the cell boundary.

[0041] Preferably, the step 4 is specifically:

[0042] The P2 polynomial is used to describe the spatial distribution of discrete triangle elements. and In the discrete triangle unit, it is expressed as:

[0043] (14)

[0044] (15)

[0045] In the formula , and Indicates k The approximate solution of degrees of freedom ( k= 1,...,6), represents the basis function. For the P2 polynomial, is a quadratic function, when the point k The value of the basis function at the above six points is 1, while the value of the basis function at the remaining five points is 0;

[0046] The triangular element in the cylindrical coordinate system is projected into the isoparametric element coordinate system. The three nodes of the triangular element in the cylindrical coordinate system are projected into the isoparametric element coordinate system to correspond to (-1,-1), (1,-1) and (-1,1) respectively. K The position of any point in the can be expressed by the following coordinate transformation expression:

[0047] (16)

[0048] Where ( ξ , η ) represents the position coordinates of the projection point in the isoparametric coordinate system, , and It is k The position coordinates of the three nodes of each unit in the cylindrical coordinate system;

[0049] Project the line element in the cylindrical coordinate system to the isoparametric element coordinate system. The two nodes of a line in the cylindrical coordinate system are projected to the isoparametric element coordinate system and correspond to -1 and 1 respectively. The coordinate transformation expression of any point on the line in the cylindrical coordinate system to the isoparametric element coordinate system is:

[0050] (17)

[0051] In the formula and They are respectively i The two nodes of an edge;

[0052] Calculate Gaussian integration points The expression of the sound speed at is:

[0053] (18)

[0054] In the formula It is k The sound velocity value at the centroid of each unit is Indicates k The centroid coordinates of the cells, represents the coordinates of the Gaussian integration points.

[0055] Preferably, the step 5 is specifically:

[0056] Calculate the numerical flux of the discrete triangular unit flow field where the sound velocity is continuous, discontinuous, and at the unit boundary of the flow field boundary and When the speed of sound c ( r , z ) is continuous, then the sound velocity value at the boundary of the triangular unit is continuous. and The calculation expression is as follows:

[0057] (19)

[0058] (20)

[0059] In the formula, {} represents the average value of the physical quantities on both sides of the unit boundary, [] represents the jump calculation of the physical quantities on the unit boundary in the normal direction, represents the average value of the product of the approximate solution of fluid dynamic pressure and the speed of sound on both sides of the unit boundary, represents the average value of the product of the auxiliary variable approximate solution and the sound speed on both sides of the unit boundary, It represents the jump calculation of the product of the approximate solution of the fluid dynamic pressure on both sides of the unit boundary and the speed of sound in the normal direction. Represents the jump calculation of the product of the auxiliary variable approximate solution and the speed of sound on both sides of the unit boundary in the normal direction.

[0060] The expression is as follows:

[0061] (twenty one)

[0062] (twenty two)

[0063] (twenty three)

[0064] (twenty four)

[0065] Where + represents the unit K Physical quantity on, – represents the unit K Edge e The physical quantity of another adjacent unit, , , sign is a mathematical symbol function, and Respectively represent units K Edge e The approximate solution of the dynamic pressure at the unit K Edge e Adjacent units on the edge e The approximate solution for the dynamic pressure at is, and Respectively represent units K Edge e The auxiliary variable approximate solution at and the unit K Edge e Adjacent units on the edge e The auxiliary variable approximate solution at and Respectively represent units K Edge e The normal vector and unit K Edge e Adjacent cell edges e The normal vector of It is the edge e Length;

[0066] When the speed of sound c ( r , z ) is discontinuous, the sound velocity value at the triangular unit boundary is discontinuous; the numerical flux at the unit boundary needs to satisfy the continuity condition and ; and Respectively represent units K Edge e The sound velocity value at the unit K Edge e Adjacent units on the edge e The sound velocity value at the point where the sound velocity is discontinuous is calculated as follows:

[0067] (25)

[0068] (26)

[0069] In the formula Represents the sound velocity value of the unit on either side of the sound velocity discontinuity interface;

[0070] and The calculation of is as follows:

[0071] (27)

[0072] (28)

[0073] In the formula and denote the approximate solutions of the numerical flux of dynamic pressure and auxiliary variables at the unit boundary, respectively. represents the average value of the approximate solution of the fluid dynamic pressure on both sides of the unit boundary, represents the average value of the auxiliary variable approximate solution on both sides of the unit boundary, It represents the jump calculation of the approximate solution of the fluid dynamic pressure on both sides of the unit boundary in the normal direction. Represents the jump calculation of the approximate solution of the auxiliary variables on both sides of the unit boundary in the normal direction;

[0074] Calculate the numerical flux at the boundary of the discrete triangular unit flow field. The calculation of the numerical flux depends on the boundary conditions. The solutions for different boundary conditions are as follows:

[0075] (29)

[0076] (30)

[0077] In the formula represents the numerical flux of the product of the approximate solution of the fluid dynamic pressure and the speed of sound at the unit boundary conditions, represents the numerical flux of the product of the approximate solution of the auxiliary variable and the speed of sound at the cell boundary, represents the density of the fluid, represents the bulk modulus of the fluid, represents the time derivative of the approximate solution of the fluid dynamic pressure at the non-reflecting boundary element, and represent the values ​​of the approximate solutions of auxiliary variables 1 and 2 at the free boundary elements, respectively.

[0078] Preferably, in step 6, the auxiliary variable 1 in the control equation is calculated, that is, :

[0079] (31)

[0080] Calculate the auxiliary variable 2 in the control equation, that is :

[0081] (32)

[0082] (33)

[0083] Where M is the mass matrix, M -1 represents the inverse matrix of the mass matrix, l m ( r , z )and l n ( r , z ) indicates the m and n Basis functions ( m , n= 1,...,6).

[0084] Preferably, the second-order derivative of the flow field dynamic pressure with respect to time is calculated in step 7, specifically:

[0085] calculate :

[0086] (34)

[0087] In the formula It represents the approximate solution of the second-order derivative of the fluid dynamic pressure with respect to time.

[0088] Preferably, the pressure cutoff model in step 8 processes the cavitation effect specifically as follows:

[0089] When the absolute pressure in the fluid is lower than the cavitation limit, cavitation will occur. The program will identify the cavitation area and mark the cavitation area. The pressure truncation model in the flow field is expressed as:

[0090] (35)

[0091] (36)

[0092] In the formula is the pseudo-dynamic pressure, is the displacement of the fluid particles, is the hydrostatic pressure of the fluid, is the cavitation limit; the cavitation limit is set to 0, indicating that the fluid cannot transmit negative pressure;

[0093] Adding artificial bulk viscosity to the dynamic pressure of the fluid , The calculation expression is as follows:

[0094] (37)

[0095] In the formula b 1 is the damping coefficient, represents the characteristic length of the triangular element, represents the area of ​​the triangular unit, is the volume strain rate.

[0096] Preferably, in step 9, the fourth-order Runge-Kutta method is used to calculate the physical quantity at the next moment, specifically:

[0097] In step 9, the fourth-order Runge-Kutta method is used to calculate the physical quantity at the next moment, specifically:

[0098] Introducing auxiliary variables , , and , then formulas (11)-(13) are expressed as:

[0099] (38)

[0100] (39)

[0101] (40)

[0102] get:

[0103] ;

[0104] Combined with initial conditions , ;

[0105] The physical quantity of the next time step is solved by the fourth-order Runge-Kutta method, which is expressed as:

[0106] (41)

[0107] (42)

[0108] (43)

[0109] (44)

[0110] (45)

[0111] In the formula for t m The fluid dynamic pressure at the moment and the first derivative of the fluid dynamic pressure with respect to time; for t m +Δ t The fluid dynamic pressure at the moment and the first derivative of the fluid dynamic pressure with respect to time; , , , They respectively represent the fluid dynamic pressure calculated at different times and under different auxiliary variable conditions and the first-order derivative of the fluid dynamic pressure with respect to time.

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

[0113] The present invention adopts the high-order local discontinuous Galerkin (LDG) method to solve the wave equation satisfied by the underwater explosion pressure load. Compared with the traditional acoustic finite element calculation method, the high-order LDG method is more suitable for dealing with strong discontinuity problems such as underwater explosions, and can more accurately capture the propagation process of underwater explosion pressure loads, effectively improving the calculation accuracy of underwater explosion loads;

[0114] The present invention can successfully simulate the propagation of underwater explosion shock waves in heterogeneous fluids. By constructing a computational model specifically for heterogeneous fluids and combining continuity conditions, an effective simulation of the propagation process of underwater explosion shock waves in a heterogeneous flow field containing a discontinuous sound velocity interface and a sound velocity gradient is achieved. This overcomes the defect that the existing technology is mostly based on the background of homogeneous fluids for calculation and cannot accurately reflect the impact of the non-uniformity of fluid sound velocity on the underwater explosion load in the real marine environment. It provides key technical support for in-depth research on underwater explosion problems in complex marine environments;

[0115] The present invention comprehensively studies the influence of non-uniformity of sound velocity in the flow field on the underwater explosion shock wave and cavitation characteristics. Through the organic combination of the above-mentioned calculation method, non-homogeneous fluid simulation and cavitation effect processing, a complete research system is formed, which can deeply explore the intrinsic connection and interaction mechanism between non-homogeneous flow field factors such as sound velocity jump layers and underwater explosion shock wave propagation characteristics and cavitation phenomena, and provide new technical means and theoretical basis for optimizing submarine combat strategies in complex marine environments and designing more effective marine engineering protection structures. BRIEF DESCRIPTION OF THE DRAWINGS

[0116] In order to more clearly illustrate the technical solutions of the present invention or the prior art, the following briefly introduces the drawings required for use in the embodiments or the prior art descriptions. Obviously, the following description is only one embodiment of the present invention, and a person skilled in the art can obtain other drawings based on these drawings without creative work.

[0117] Figure 1 It is a flow chart of the overall calculation process of the present invention;

[0118] Figure 2 It is a schematic diagram of the splitting of a rectangular unit of the present invention;

[0119] Figure 3 The present invention simplifies the three-dimensional flow field model into an axisymmetric model schematic diagram;

[0120] Figure 4 Schematic diagram of the transformation from cylindrical coordinate system to isoparametric coordinate system: (a) two-dimensional transformation (b) one-dimensional transformation;

[0121] Figure 5 It is a schematic diagram of the propagation of underwater explosion shock waves in a non-homogeneous flow field;

[0122] Figure 6 This is the cavitation diagram of underwater explosion shock wave in non-homogeneous flow field. DETAILED DESCRIPTION

[0123] The axisymmetric calculation model of underwater explosion load in heterogeneous fluid established by the local discontinuous Galerkin (LDG) method and continuity conditions in the present invention has many significant advantages. First, the model overcomes the limitations of the traditional homogeneous fluid assumption, provides a powerful tool for in-depth exploration of underwater explosion phenomena in complex marine environments, and promotes the further development and improvement of theoretical research and practical applications related to underwater explosions. Second, it can accurately simulate the sound velocity discontinuous interface of the flow field, accurately analyze the influence of the sound velocity jump layer interface on the propagation path and load peak of the underwater explosion shock wave, thereby providing more reliable underwater explosion load data for submarine combat operations near the sound velocity jump layer. Third, it can effectively study the influence of the flow field sound velocity gradient on the underwater explosion load and cavitation characteristics, which is helpful in the design of marine engineering protection, based on more accurate underwater explosion load input data to optimize the design and construction plan of protective ships and marine structures, enhance the ability of protective facilities to resist underwater explosion threats, and improve the safety and stability of marine engineering facilities.

[0124] The following will further introduce the specific implementation method as well as the technical difficulties and inventive points of this invention in combination with this design example.

[0125] This paper proposes a method for calculating the propagation of underwater explosion shock waves in non-uniform acoustic flow fields. This method uses the discontinuous Galerkin (LDG) method to solve the dynamic pressure of the flow field, and then studies the propagation characteristics of the far-field underwater explosion shock wave in the non-uniform flow field. By comparing the flow field dynamic pressure calculated by this method with the analytical solution of standing waves in non-homogeneous flow fields, the results show that the calculation results of this method can achieve third-order accuracy.

[0126] Combination Figure 1 to Figure 6 and Table 1, a calculation method for propagation of underwater explosion shock waves in a non-uniform acoustic flow field according to the present invention, the steps are as follows.

[0127] Step 1: Based on the flow field range, size and corresponding rectangular mesh node and unit information data (stored in the exported inp file) established and discretized by ABAQUS finite element software, the rectangular unit is split and the new nodes and units are numbered to obtain the triangular unit flow field data results containing numbering information. This provides basic data for further calculations based on the numbered triangular unit flow field.

[0128] First, the flow field range and flow field size are established in ABAQUS finite element software. The entire flow field is discretized into rectangular grids using ABAQUS finite element software, and each rectangular grid is a unit. The discrete flow field is exported as an inp file, where the exported inp file includes the node numbers and coordinate positions of each node in the flow field and the node composition of each unit. After the program reads the inp file exported by ABAQUS, it combines Figure 2 , splitting each rectangular unit into 4 triangular units. The program will number the newly generated nodes and units and use Ω Represents the discretized triangular unit flow field. The data of the flow field will be used as input data for subsequent calculation processes.

[0129] Step 2: This step sets the incident wave loading surface, explosion point and source point based on the discrete triangular unit flow field data obtained in step 1, and then uses the Geers-Hunter model and the far-field underwater explosion pressure calculation formula to calculate the time-load data at the explosion point and the dynamic pressure data on the incident wave loading surface of the discrete triangular unit flow field. These calculated data will provide the necessary data basis for subsequent operations such as underwater explosion shock wave propagation analysis in the flow field.

[0130] The present invention adopts a calculation model of far-field underwater explosion, where the explosive is located in a discrete triangular unit flow field. Ω After the explosive explodes, the underwater explosion shock wave enters the discrete triangular unit flow field through the incident wave loading surface Ω . Combined Figure 3 First, the discrete triangle unit flow field explosion point and source point and the position of the incident wave loading surface are set in the program. For the shock wave generated by the underwater explosion of explosives, the present invention adopts the Geers-Hunter model to calculate. Figure 3 , calculate the time-load data of explosive explosion at the explosion point in the discrete triangular unit flow field p t After the calculation is completed, the time-load data at the explosion point p t Import into the program. The program reads the time-load data of the explosion point of the explosive explosion in the discrete triangular unit flow field p t . The flow field in the discrete triangular unit Ω Any point x on the incident wave loading surface j The calculation formula for the gradient of the dynamic pressure at is:

[0131] (1)

[0132] (2)

[0133] (3)

[0134] Mode where x represents the gradient of the dynamic pressure on the incident wave loading surface, s and x 0 Represents the vector position coordinates of the explosive source point and the explosion point, x j Represents the vector position coordinates of the point on the incident wave loading surface, R 0 and R j They represent the distance between the explosive source and the explosion point and the distance between any point on the incident wave loading surface and the explosive source, respectively. for p t The derivative with respect to time.

[0135] Step 3, based on the dynamic pressure data on the incident wave loading surface of the discrete triangular unit flow field calculated in step 2, the LDG method is used to derive the basic equations for the propagation of underwater explosion shock waves in the discrete triangular unit flow field. Thus, the control equations that can be used to calculate the propagation of underwater explosion shock waves in the discrete triangular unit flow field are obtained. The calculated control equations provide a specific calculation method basis for the subsequent calculation of unit dynamic pressure and shock wave propagation characteristics in the flow field.

[0136] The fluid disturbed by the underwater explosion shock wave can be regarded as an inviscid, compressible, and adiabatic fluid. Step 2 has obtained the dynamic pressure on the incident wave loading surface. The propagation of the underwater explosion shock wave in the fluid satisfies the second-order wave equation, and the dynamic pressure at other positions in the flow field can be obtained.

[0137] S3-1. The calculation formula of the second-order wave equation is:

[0138] (4)

[0139] In the formula p is the dynamic pressure of the fluid, Δ represents the Laplace operator, represents the second derivative of the fluid's dynamic pressure with respect to time, c The fluid in the flow field Ω The speed of sound inside;

[0140] S3-2, Combination Figure 3 , the present invention simplifies the three-dimensional flow field model into an axisymmetric model. In the axisymmetric model under the cylindrical coordinate system, calculation formula 4 can be written as:

[0141] (5)

[0142] In the formula rand z are the radial coordinate and the axial coordinate, c ( r , z ) represents the flow field in the axisymmetric model Ω The speed of sound at any point in The dynamic pressure of the fluid r Direction partial derivative, The dynamic pressure of the fluid r Direction to find the second-order partial derivative, The dynamic pressure of the fluid z Direction to find the second-order partial derivative;

[0143] S3-3. Auxiliary variables Introduced into the axisymmetric model in the cylindrical coordinate system, the axisymmetric model of the second-order wave equation in the cylindrical coordinate system is expressed as:

[0144] (6)

[0145] (7)

[0146] Represents the derivative operation symbol, q 1 The dynamic pressure of the fluid on the space r Take the partial derivative of the direction and multiply it by the speed of sound c ; q 2 The dynamic pressure of the fluid on the space z Take the partial derivative of the direction and multiply it by the speed of sound c ;

[0147] S3-4, multiply the second-order wave equation of the above axisymmetric model by the basis function ψ , ϕ and φ , and integrating the discrete units, we can get:

[0148] (8)

[0149] (9)

[0150] (10)

[0151] In the formula K Representing flow field Ω Discrete triangular elements in ; Represents the product of auxiliary variable 1 and the speed of sound r Direction partial derivative; The product of auxiliary variable 2 and the speed of sound is z Direction partial derivative; and Respectively represent the fluid dynamic pressure r and z Find partial derivatives;

[0152] S3-5. The control equation for the propagation of underwater explosion shock waves in discrete triangular unit flow fields can be obtained through the rule of partial integration and Gaussian formula. The control equation can be expressed as:

[0153] (11)

[0154] (12)

[0155] (13)

[0156] In the formula Represents the fluid dynamic pressure in discrete triangular elements K The approximate solution of and Respectively represent the auxiliary variables in the discrete triangle unit K The approximate solution of ∂K represents the boundary of the discrete triangle unit, is the unit external normal vector to the cell boundary, and Respectively represent the test function y right r and z Find the partial derivative, Represents the test function f right r Find the partial derivative, Represents the test function j right z Find the partial derivative, and They represent the flow field sound velocity r and z Directional derivative, , and is the numerical flux at the cell boundary.

[0157] Step 4, based on the control equation of the underwater explosion shock wave propagation in the discrete triangular unit flow field in step 3. The calculation accuracy is improved by multiplying each physical quantity by the basis function, and the calculation of each integral term in equations 11-13 in step 3 is converted into a Gaussian integral formula calculation through coordinate transformation, and the sound speed value at each Gaussian integral point is calculated. This provides an important basis for the subsequent calculation of the shock wave propagation characteristics in the flow field and other related calculations.

[0158] Further solve equations 11-13 in step 3.

[0159] S4-1. For the LDG method, a more accurate solution can be obtained by increasing the order of the polynomial. In the present invention, the P2 polynomial is used to describe the spatial distribution of discrete triangle units. In formulas 11-13, and In the discrete triangle unit, it is expressed as:

[0160] (14)

[0161] (15)

[0162] In the formula , and Indicates k The approximate solution of degrees of freedom ( k = 1,...,6), represents the basis function. For the P2 polynomial, is a quadratic function, when the point k The value of the basis function at the above six points is 1, while the value of the basis function at the remaining five points is 0;

[0163] S4-2. Combination Figure 4 , project the triangular element in the cylindrical coordinate system into the isoparametric element coordinate system. The three nodes of the triangular element in the cylindrical coordinate system projected into the isoparametric element coordinate system correspond to (-1,-1), (1,-1) and (-1,1) respectively. K The position of any point in the can be expressed by the following coordinate transformation expression:

[0164] (16)

[0165] Where ( ξ , η ) represents the position coordinates of the projection point in the isoparametric coordinate system, , and It is k The position coordinates of the three nodes of each unit in the cylindrical coordinate system;

[0166] S4-3. Combination Figure 4 , project the line element in the cylindrical coordinate system into the isoparametric element coordinate system, and the two nodes of a line in the cylindrical coordinate system are projected into the isoparametric element coordinate system to correspond to -1 and 1 respectively. The coordinate transformation expression of any point on the line in the cylindrical coordinate system to the isoparametric element coordinate system is:

[0167] (17)

[0168] In the formula and They are respectively i The two nodes of an edge;

[0169] S4-4. After completing the isoparametric coordinate transformation, the integral terms in formulas 11-13 in step 3 can be calculated using the Gaussian integral formula. Calculate the Gaussian integral point The expression of the sound speed at is:

[0170] (18)

[0171] In the formula It is k The sound velocity value at the centroid of each unit is Indicates k The centroid coordinates of the cells, represents the coordinates of the Gaussian integration points.

[0172] Step 5: This step is based on the control equation of underwater explosion shock wave propagation in discrete triangular unit flow field in step 3 and the calculation method of integral term and the sound velocity value at Gaussian integral point in step 4. This step calculates the numerical flux of the discrete triangular unit flow field at the unit boundary where the sound velocity is continuous, discontinuous and at the flow field boundary. and This provides important data input for subsequent calculations of shock wave propagation characteristics in the flow field and other related calculations.

[0173] S5-1, when the speed of sound c ( r , z ) is continuous, then the sound velocity value at the boundary of the triangular unit is continuous. For the equations 11-13 in step 3 and The calculation expression is as follows:

[0174] (19)

[0175] (20)

[0176] In the formula, {} represents the average value of the physical quantities on both sides of the unit boundary, and [] represents the jump calculation expression of the physical quantity of the unit boundary in the normal direction is as follows:

[0177] (twenty one)

[0178] (twenty two)

[0179] (twenty three)

[0180] (twenty four)

[0181] Where + represents the unit K Physical quantity on, – represents the unit K Edge e The physical quantity of another adjacent unit, , , sign is a mathematical symbol function, and Respectively represent units K Edge e The approximate solution of the dynamic pressure at the unit K Edge e Adjacent units on the edge e The approximate solution for the dynamic pressure at is, and Respectively represent units K Edge e The auxiliary variable approximate solution at and the unit K Edge e Adjacent units on the edge e The auxiliary variable approximate solution at and Respectively represent units K Edge e The normal vector and unit K Edge e Adjacent cell edges e The normal vector of It is the edge e Length;

[0182] S5-2, when the speed of sound c ( r , z ) is discontinuous, the sound velocity value at the triangular unit boundary is discontinuous; the numerical flux at the unit boundary needs to satisfy the continuity condition and ; and Respectively represent units K Edge e The sound velocity value at the unit K Edge e Adjacent units on the edge e The sound velocity value at the point where the sound velocity is discontinuous is calculated as follows:

[0183] (25)

[0184] (26)

[0185] In the formula c iRepresents the sound velocity value of the unit on either side of the sound velocity discontinuity interface;

[0186] and The calculation of is as follows:

[0187] (27)

[0188] (28)

[0189] S5-3. Calculate the numerical flux at the boundary of the discrete triangular unit flow field. The calculation of the numerical flux depends on the boundary conditions. The solutions for different boundary conditions are as follows:

[0190] (29)

[0191] (30)

[0192] In the formula represents the bulk modulus of the fluid, represents the time derivative of the approximate solution of the fluid dynamic pressure at the non-reflecting boundary element, and represent the values ​​of the approximate solutions of auxiliary variables 1 and 2 at the free boundary elements, respectively.

[0193] Step 6, based on the control equation of underwater explosion shock wave propagation in discrete triangular unit flow field in step 3, the calculation method of integral term and the sound velocity value at Gaussian integral point in step 4, and the numerical flux results at all unit boundaries calculated in step 5. This step calculates the formulas 12 and 13 in step 3. and The result of the calculation is the subsequent calculation Provided important data input.

[0194] S6-1. Calculation :

[0195] (31)

[0196] S6-2. Calculation :

[0197] (32)

[0198] (33)

[0199] Where M is the mass matrix, M -1 represents the inverse matrix of the mass matrix, l m ( r , z)and l n ( r , z ) indicates the m and n Basis functions ( m , n = 1,...,6).

[0200] Step 7, based on the control equation of underwater explosion shock wave propagation in the flow field derived in step 3, the calculation method of the integral term and the sound velocity value at the Gaussian integration point in step 4, the numerical flux results calculated at all unit boundaries in step 5, and the results in step 6 and This step calculates the second-order derivative of the flow field dynamic pressure with respect to time in formula 11 in step 3. The calculation results provide important input data for the subsequent Runge-Kutta cycle;

[0201] calculate :

[0202] (34).

[0203] Step 8: Based on the relevant data obtained from the calculations completed in steps 3-7, the dynamic pressure results of the discrete triangular unit flow field are calculated. On this basis, the pressure truncation model is introduced to deal with the cavitation effect, and artificial volume viscosity is added to the calculated dynamic pressure results. , thus providing a reliable basis for flow field characteristics analysis.

[0204] S8-1. The pressure cutoff model handles the cavitation effect. When the absolute pressure in the fluid is lower than the cavitation limit, cavitation will occur. The program will identify the cavitation area and mark the cavitation area. The pressure cutoff model in the flow field can be expressed as:

[0205] (35)

[0206] (36)

[0207] In the formula is the pseudo-dynamic pressure, is the displacement of the fluid particles, The hydrostatic pressure of the fluid, is the cavitation limit; the cavitation limit is set to 0, indicating that the fluid cannot transmit negative pressure;

[0208] S8-2. Add artificial volume viscosity. In order to eliminate the non-physical oscillation caused by shock waves and cavitation effects, and obtain more reasonable, accurate and stable dynamic pressure data results, artificial volume viscosity is added to the dynamic pressure of the fluid. , The calculation expression is as follows:

[0209] (37)

[0210] In the formula b 1 is the damping coefficient, represents the characteristic length of the triangular element, represents the area of ​​the triangular unit, is the volume strain rate.

[0211] Step 9, based on the relevant calculations and processed data completed in steps 3-8, the physical quantities of all units of the discrete triangular unit flow field at the initial moment can be obtained. The fourth-order Runge-Kutta method is used to calculate the physical quantities of all units of the discrete triangular unit flow field at the next moment. This step is repeated continuously to finally achieve the simulation of the propagation characteristics of underwater explosion shock waves in non-homogeneous flow fields, including the influence of the sound velocity discontinuous interface on the propagation path and peak value of underwater explosion shock waves, and the influence of the sound velocity gradient on the underwater explosion load peak value and flow field cavitation characteristics.

[0212] Time step update, using the fourth-order Runge-Kutta method to calculate the physical quantity at the next moment:

[0213] Introducing auxiliary variables , , and , then formulas (11)-(13) are expressed as:

[0214] (38)

[0215] (39)

[0216] (40)

[0217] (40)

[0218] (41)

[0219] Formula (38) can be expressed as:

[0220] ;

[0221] Combined with initial conditions , ;

[0222] The physical quantity of the next time step is solved by the fourth-order Runge-Kutta method, which is expressed as:

[0223] 1. From the equations , Calculated , and calculate ;

[0224] 2. From the system of equations , Calculated , , and calculate ;

[0225] 3. From the system of equations , calculate , , and calculate ;

[0226] 4. From the system of equations , calculate , , and calculate ;

[0227] 5. Final calculation .

[0228] In the formula for t m The fluid dynamic pressure at the moment and the first derivative of the fluid dynamic pressure with respect to time; for t m +Δ t The fluid dynamic pressure at the moment and the first derivative of the fluid dynamic pressure with respect to time; , , , They respectively represent the fluid dynamic pressure calculated at different times and under different auxiliary variable conditions and the first-order derivative of the fluid dynamic pressure with respect to time.

[0229] As shown in Table 1, the calculation method of the present invention is compared with the analytical solution of the propagation of standing waves in a non-homogeneous flow field. The calculation results of the method of the present invention can achieve third-order accuracy, which verifies the accuracy of the method. Figure 5 and 6 As shown, the propagation diagram of underwater explosion shock wave in the inhomogeneous flow field and the cavitation diagram of the inhomogeneous flow field are given at some typical moments.

[0230] Table 1 Error and accuracy of calculated results and standing wave analytical solution

[0231] Grid Resolution error Accuracy 40×20 8.87e-3 80×40 9.10e-4 3.28 160×80 1.35e-4 2.75 320×160 1.87e-5 2.85

[0232] The above description is only a preferred embodiment of the present application and is not intended to limit the present application. For those skilled in the art, the present application may have various modifications and variations. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present application shall be included in the protection scope of the present application.

[0233] Although the above describes the specific implementation methods of the present invention, it is not intended to limit the scope of protection of the present invention. Those skilled in the art should understand that various modifications or variations that can be made by those skilled in the art on the basis of the technical solution of the present invention without creative work are still within the scope of protection of the present invention.

Claims

1. A method for calculating the propagation of underwater explosion shock waves in a non-uniform acoustic flow field, characterized in that: The process includes: Step 1, establish the flow field range and size, and discretize the flow field into a rectangular grid; based on the corresponding rectangular grid node and unit information data, split the rectangular unit and number the new node and unit to obtain the triangular unit flow field data containing the number information; Step 2: Based on the triangular unit flow field data, the incident wave loading surface, explosion point and source point positions are set, and then the Geers-Hunter model and the far-field underwater explosion pressure calculation formula are used to calculate the time-load data at the explosion point and the dynamic pressure data on the discrete triangular unit flow field incident wave loading surface; Step 3, using the LDG method to derive the second-order wave equation for the propagation of underwater explosion shock waves in discrete triangular unit flow fields, and obtain the control equation for calculating the propagation of underwater explosion shock waves in discrete triangular unit flow fields; the propagation of underwater explosion shock waves in fluids satisfies the second-order wave equation, and the calculation formula of the second-order wave equation is: (1) In the formula p is the dynamic pressure of the fluid, Δ represents the Laplace operator, represents the second derivative of the fluid's dynamic pressure with respect to time, c The fluid in the flow field Ω The speed of sound inside; Auxiliary variables Introduced into the axisymmetric model in the cylindrical coordinate system, the axisymmetric model of the second-order wave equation in the cylindrical coordinate system is expressed as: (2) (3) Represents the derivative operation symbol, q 1 The dynamic pressure of the fluid on the space r Take the partial derivative of the direction and multiply it by the speed of sound c ; q 2 The dynamic pressure of the fluid on the space z Take the partial derivative of the direction and multiply it by the speed of sound c ; Multiply both sides of the equation by the test function ψ , ϕ and φ , the control equation of underwater explosion shock wave propagation in discrete triangular unit flow field is obtained by the partial integration operation rule and Gaussian formula. The control equation is expressed as: (4) (5) (6) In the formula Represents the fluid dynamic pressure in discrete triangular elements K The approximate solution of and Respectively represent the auxiliary variables in the discrete triangle unit K The approximate solution of ∂K represents the boundary of the discrete triangle unit, is the unit external normal vector to the cell boundary, and Represent the test function right r and z Find the partial derivative, Represents the test function right r Find the partial derivative, Represents the test function right z Find the partial derivative, and They represent the flow field sound velocity r and z Directional derivative, , and is the numerical flux at the cell boundary; Step 4: Based on the obtained control equation, the calculation of each integral term in the control equation is converted into a Gaussian integral formula calculation through coordinate transformation, and the sound velocity value at each Gaussian integral point is calculated; wherein the coordinate transformation and the calculation of the sound velocity value at the Gaussian integral point are specifically as follows: The triangular element in the cylindrical coordinate system is projected into the isoparametric element coordinate system. The three nodes of the triangular element in the cylindrical coordinate system are projected into the isoparametric element coordinate system to correspond to (-1,-1), (1,-1) and (-1,1) respectively. K The position of any point in the can be expressed by the following coordinate transformation expression: (7) Where ( ξ , η ) represents the position coordinates of the projection point in the isoparametric coordinate system, , and It is k The position coordinates of the three nodes of each unit in the cylindrical coordinate system; Project the line element in the cylindrical coordinate system to the isoparametric element coordinate system. The two nodes of a line in the cylindrical coordinate system are projected to the isoparametric element coordinate system and correspond to -1 and 1 respectively. The coordinate transformation expression of any point on the line in the cylindrical coordinate system to the isoparametric element coordinate system is: (8) In the formula and They are respectively i The two nodes of an edge; Calculate Gaussian integration points The expression of the sound speed at is: (9) In the formula It is k The sound velocity value at the centroid of each unit is Indicates k The centroid coordinates of the cells, represents the coordinates of the Gaussian integration points; Step 5, calculate the numerical flux of the discrete triangular unit flow field where the sound velocity is continuous, discontinuous and at the unit boundary of the flow field; when the sound velocity of the flow field is discontinuous, the specific calculation is: When the speed of sound c ( r , z ) is discontinuous, the sound velocity value at the triangular unit boundary is discontinuous; the numerical flux at the unit boundary needs to satisfy the continuity condition and ; and Respectively represent units K Edge e The sound velocity value at the unit K Edge e Adjacent units on the edge e The sound velocity value at the point where the sound velocity is discontinuous is calculated as follows: (10) (11) In the formula c i Represents the sound velocity value of the unit on either side of the sound velocity discontinuity interface; and The calculation of is as follows: (12) (13) In the formula and denote the approximate solutions of the numerical flux of dynamic pressure and auxiliary variables at the unit boundary, respectively. represents the average value of the approximate solution of the fluid dynamic pressure on both sides of the unit boundary, represents the average value of the auxiliary variable approximate solution on both sides of the unit boundary, It represents the jump calculation of the approximate solution of the fluid dynamic pressure on both sides of the unit boundary in the normal direction. Represents the jump calculation of the approximate solution of the auxiliary variables on both sides of the unit boundary in the normal direction; Step 6: Calculate the approximate solutions of auxiliary variables 1 and 2 in the control equation. Auxiliary variable 1 is the dynamic pressure of the fluid on the space. r Directional derivative and multiply by the speed of sound c , auxiliary variable 2 is the dynamic pressure of the fluid on the space z Directional derivative and multiply by the speed of sound c ; Step 7, calculating the second derivative of the flow field dynamic pressure with respect to time; Step 8, calculating the dynamic pressure results of the discrete triangular unit flow field, introducing a pressure truncation model to deal with the cavitation effect during the calculation, and adding artificial volume viscosity to the calculated dynamic pressure results; Step 9, based on the calculation results of steps 3 to 8, obtain the physical quantities of all units of the discrete triangular unit flow field at the initial moment, and use the fourth-order Runge-Kutta method to calculate the physical quantities of all units of the discrete triangular unit flow field at the next moment, so as to simulate the propagation characteristics of the underwater explosion shock wave in the inhomogeneous flow field.

2. The method for calculating the propagation of underwater explosion shock waves in a non-uniform acoustic flow field according to claim 1, characterized in that: In the step 1, first, the flow field range and flow field size are established in the ABAQUS finite element software, and the entire flow field is discretized into rectangular grids using the ABAQUS finite element software, each rectangular grid is a unit; each rectangular unit is split into 4 triangular units, and the newly generated nodes and units are numbered and Represents the discretized triangular unit flow field.

3. The method for calculating the propagation of underwater explosion shock waves in a non-uniform acoustic flow field as claimed in claim 2, characterized in that: In step 2, the positions of the discrete triangular unit flow field explosion point, source point and incident wave loading surface are set, and the Geers-Hunter model is used to calculate the time-load data of the explosion of the explosive at the discrete triangular unit flow field explosion point for the shock wave generated by the underwater explosion of the explosive. p t After the calculation is completed, the time-load data at the explosion point p t Import the program, located in the discrete triangular unit flow field Ω Any point x on the incident wave loading surface j The calculation formula for the gradient of the dynamic pressure at is: (14) (15) (16) Mode represents the gradient of the dynamic pressure on the incident wave loading surface, and Represent the vector position coordinates of the explosive source point and the explosion point, Represents the vector position coordinates of the point on the incident wave loading surface, R 0 and R j They represent the distance between the explosive source and the explosion point and the distance between any point on the incident wave loading surface and the explosive source, respectively. ṗ t for p t The derivative with respect to time.

4. The method for calculating the propagation of underwater explosion shock waves in a non-uniform acoustic flow field according to claim 1, characterized in that: In step 6, the auxiliary variable 1 in the control equation is calculated, that is, : (17) Calculate the auxiliary variable 2 in the control equation, that is : (18) (19) Where M is the mass matrix, M -1 represents the inverse matrix of the mass matrix, l m ( r , z )and l n ( r , z ) indicates the m and n Basis Functions .

5. The method for calculating the propagation of underwater explosion shock waves in a non-uniform acoustic flow field as claimed in claim 4, characterized in that: The second-order derivative of the flow field dynamic pressure with respect to time is calculated in step 7, specifically: calculate : (20) In the formula It represents the approximate solution of the second-order derivative of the fluid dynamic pressure with respect to time.

6. The method for calculating the propagation of underwater explosion shock waves in a non-uniform acoustic flow field as claimed in claim 5, characterized in that: The pressure cutoff model in step 8 processes the cavitation effect specifically as follows: When the absolute pressure in the fluid is lower than the cavitation limit, cavitation will occur. The program will identify the cavitation area and mark the cavitation area. The pressure truncation model in the flow field is expressed as: (21) (22) In the formula is the pseudo-dynamic pressure, is the displacement of the fluid particles, is the hydrostatic pressure of the fluid, is the cavitation limit; the cavitation limit is set to 0, indicating that the fluid cannot transmit negative pressure; Added artificial bulk viscosity to the dynamic pressure of the fluid , The calculation expression is as follows: (23) In the formula b 1 is the damping coefficient, represents the characteristic length of the triangular element, represents the area of ​​the triangular unit, is the volume strain rate.

7. The method for calculating the propagation of underwater explosion shock waves in a non-uniform acoustic flow field as claimed in claim 6, characterized in that: In step 9, the fourth-order Runge-Kutta method is used to calculate the physical quantity at the next moment, specifically: In step 9, the fourth-order Runge-Kutta method is used to calculate the physical quantity at the next moment, specifically: Introducing auxiliary variables , , and , then formulas (4)-(6) are expressed as: (24) (25) (26) get: (27) Combined with initial conditions , ; The physical quantity of the next time step is solved by the fourth-order Runge-Kutta method, which is expressed as: （28） (29) (30) (31) (32) In the formula for t m The fluid dynamic pressure at the moment and the first derivative of the fluid dynamic pressure with respect to time; for t m +Δ t The fluid dynamic pressure at the moment and the first derivative of the fluid dynamic pressure with respect to time; , , , They respectively represent the fluid dynamic pressure calculated at different times and under different auxiliary variable conditions and the first-order derivative of the fluid dynamic pressure with respect to time.

Images