Underwater explosion bubble movement forecasting method based on discontinuous Golgu method

Abstract

The invention discloses an underwater explosion bubble motion forecasting method based on a discontinuous Cujin method, and relates to the technical field of explosion, and the method comprises the steps: carrying out the zero setting of a viscosity item and a source item in an unsteady viscosity uniform two-phase flow Euler equation, multiplying a smooth function, carrying out the integration on a fluid unit, introducing a divergence theory, and constructing a system control equation, polynomial extension is carried out on the flux and the smooth function by using the primary functions, and then the polynomial extension is substituted into a system control equation to obtain discrete control equations corresponding to the primary functions; the polynomial extension form of the flux is substituted into a semi-discrete equation of an unsteady viscosity uniform two-phase flow Euler equation for solving, so that the two-phase flow state of each time step can be obtained, and the propagation of shock waves in the early stage of underwater explosion bubble movement and the bubble movement process can be forecasted; interface capture in the underwater near-field explosion process can be achieved, and a solid theoretical support is provided for deep research of underwater explosion impact dynamics and bubble dynamics characteristics.

Description

Technical Field

[0001] The present application relates to the field of explosion technology, and in particular to a method for predicting underwater explosion bubble movement based on the intermittent Galukin method. Background Art

[0002] As marine security issues become increasingly prominent, countries around the world are competing to develop advanced underwater weapons with high speed, large charge and precision guidance. These advanced underwater weapons pose a great threat to underwater ships. Therefore, studying the underwater near-field explosion characteristics of advanced underwater weapons is of great significance to the damage and safety protection design of underwater ship structures.

[0003] In the research of underwater near-field explosion, the prediction of the propagation of the shock wave and the movement of bubbles generated by the explosion has always been a difficult problem. Currently, the following methods are commonly used:

[0004] (1) The most widely used method is the boundary element method (BEM). Based on the BEM method, Professor Zhang Aman's team has conducted extensive research on the three-dimensional bubble motion characteristics under complex boundaries such as near-free surfaces, rigid walls, and elastic walls. The relationship between bubble motion characteristics and characteristic parameters is given, laying the foundation for the damage of bubble loads to ship structures. However, when dealing with the penetration problem after the bubble jet and the coupling effect between bubbles and complex structure interfaces, the boundary element method needs to cut and sew the grid, which will reduce the solution accuracy and make the algorithm complicated.

[0005] (2) Many scholars have also used the grid method based on Euler expression to study the evolution of bubble pulsation. For example, Liu et al. used the Euler finite element method (EFEM) to study the bubble pulsation problem of free-field underwater explosion. However, the grid method based on Euler expression requires an additional interface capture algorithm to realize the tracking of multiphase interfaces.

[0006] (3) The meshless smoothed particle hydrodynamics (SPH) method can automatically track complex material interfaces due to its Lagrangian characteristics, and has been applied by many scholars to solve underwater explosion problems. For example, Liu et al. calculated the two-dimensional underwater explosion process based on the SPH method. Chen Juan et al. proposed a new numerical research method for large deformation, high inhomogeneity, deformable boundaries and free surface problems based on the SPH method. Wang Pingping and Zhang Aman et al. eliminated the influence of shock wave reflection on the pulsation of underwater explosion bubbles by correcting the non-reflecting boundary. However, due to its Lagrangian characteristics, the meshless SPH method has great challenges in the simulation of strongly compressible problems, because the distance between the initially distributed explosive particles will increase sharply after the detonation and expansion, resulting in chaotic particle distribution and a decrease in calculation accuracy, and finally the calculation cannot be sustained.

[0007] In summary, the various models currently widely used have their own defects in predicting the movement of underwater explosion bubbles, and are unable to accurately and effectively predict the propagation of the early shock wave of underwater explosion bubble movement and the bubble movement process. Summary of the invention

[0008] In view of the above problems and technical requirements, this application proposes a method for predicting underwater explosion bubble movement based on the intermittent Galekin method. The technical solution of this application is as follows:

[0009] A method for predicting underwater explosion bubble movement based on the intermittent Galogen method, the method comprising:

[0010] The viscosity term and source term in the unsteady viscous uniform two-phase flow Euler equation are set to 0 and then multiplied by the smooth function ω(r) and integrated over the fluid unit Ω. The divergence theory is introduced into the flux in the unsteady viscous uniform two-phase flow Euler equation to obtain the system control equation; where r represents the position vector of the space point;

[0011] The flux U and smooth function ω(r) in the unsteady viscous uniform two-phase flow Euler equation are expanded into polynomials using K+1 basis functions: φ k is the kth basis function, c k is the polynomial expansion of the flux U with the kth basis function φ k The corresponding coefficient, ω k is the polynomial expansion of the smooth function ω(r) with the kth basis function φ k The corresponding coefficients, K is an integer parameter;

[0012] Substitute the flux U and the polynomial expansion form of the smooth function ω(r) in the unsteady viscous uniform two-phase flow Euler equation into the system control equation to obtain the discrete control equation corresponding to each basis function. k The corresponding discrete control equation is rewritten as a partial differential equation to obtain the coefficient c k The solution expression of ;

[0013] Each coefficient c k The solution expression is applied to the polynomial expansion form of flux U and substituted into the semi-discrete equation of unsteady viscous uniform two-phase flow Euler equation, and the numerical solution of flux U in the unsteady viscous uniform two-phase flow Euler equation is obtained; the motion of underwater explosion bubbles is predicted according to the two-phase flow state indicated by the numerical solution of the flux in the unsteady viscous uniform two-phase flow Euler equation at each time step.

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

[0015] The present application discloses a method for predicting the motion of underwater explosion bubbles based on the discontinuous Galkin method. The method is based on the Euler equations for unsteady viscous uniform two-phase flow and introduces the discontinuous Galkin method for motion prediction. The method can predict the propagation of shock waves in the early stage of underwater explosion bubble motion and the bubble motion process, and can realize the interface capture of the underwater near-field explosion process, providing a solid theoretical support for the in-depth study of underwater explosion impact dynamics and bubble dynamics characteristics. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] Figure 1 It is a flow chart of a method for predicting underwater explosion bubble movement according to an embodiment of the present application.

[0017] Figure 2 The figure is a comparison diagram of the distribution curves of the fluid density predicted by the method of the present application and the finite volume method in an underwater explosion example.

[0018] Figure 3 yes Figure 2 A comparison diagram of the distribution curves of fluid pressure predicted by the method of the present application and the finite volume method in an underwater explosion example.

[0019] Figure 4 This is a result diagram of fluid density predicted by AUTODYN-2D in another underwater explosion example.

[0020] Figure 5 yes Figure 4 Result diagram of fluid density predicted by the method of the present application in an underwater explosion example.

[0021] Figure 6 yes Figure 4 and Figure 5In the underwater explosion example, the distribution curves of fluid density at two typical measuring points in the calculation domain predicted by AUTODYN-2D and the method of the present application are compared. DETAILED DESCRIPTION

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

[0023] This application discloses a method for predicting underwater explosion bubble motion based on the intermittent galaxy method. Figure 1 The schematic flow chart shown in FIG. 1 shows, the underwater explosion bubble motion prediction method includes:

[0024] Step 1: The system control equation is constructed based on the Euler equation for unsteady viscous uniform two-phase flow combined with the divergence theory.

[0025] The system control equation of the present application is constructed based on the unsteady viscous uniform two-phase flow Euler equation, which is:

[0026]

[0027] in, F u is the flux of the viscosity term in the x direction, G v is the viscous flux in the y direction, The x direction and the y direction are two coordinate axes of a pre-established plane coordinate system.

[0028] In the above formula, g is the acceleration due to gravity, and t is the time. α1 is the volume fraction of the gas phase, and ρ1 is the density of the gas phase. α2 is the volume fraction of the liquid phase, and ρ2 is the density of the liquid phase. In the results calculated by the dissipative interface tracking algorithm, the material interface is distributed in multiple fluid units. The fluid in the fluid unit at this dissipative region is regarded as a mixed product of the gas phase and the liquid phase, and the volume fractions of the gas phase and the liquid phase in the fluid satisfy α1+α2=1. ρ is the total density of the fluid and ρ=ρ1α1+ρ2α2. u is the x-direction component of the fluid velocity vector u, and v is the y-direction component of the fluid velocity vector u. p is the fluid pressure. d is the characteristic length, for example, the pipe diameter can be used as the characteristic length. μ is the viscosity of the fluid and μ=α1μ1+α2μ2, where μ1 is the viscosity coefficient of the gas phase and μ2 is the viscosity coefficient of the liquid phase.

[0029] E is the internal energy per unit mass of fluid and e is the specific internal energy, u is the fluid velocity vector and u·u=u 2 +v 2 Based on the isobaric assumption and the Stiffened state equation, it can be determined that the fluid pressure p satisfies:

[0030] p=(γ-1)ρe-γB (2)

[0031] In the above formula (2), γ and B are the parameters of the fluid in the Stiffened state equation, and can be expressed using the gas phase and liquid phase as follows:

[0032]

[0033] In the above equations (3) and (4), γ1 and B1 are the parameters of the gas phase in the Stiffened state equation, and γ2 and B2 are the parameters of the liquid phase in the Stiffened state equation, and they are usually known parameters. For example, γ1 = 1.4 and B1 = 0 for air, γ2 = 4.4 and B2 = 6×10 8 Pa.

[0034] Substituting the above equations (2), (3) and (4) into the expression of the internal energy E per unit mass of fluid, we can obtain that the internal energy per unit mass of fluid is expressed by the fluid pressure p, the fluid velocity vector u and the volume fraction of the gas and liquid phases as follows:

[0035]

[0036] From the above introduction, it can be seen that the flux U in the unsteady viscous uniform two-phase flow Euler equation contains the parameters of the two-phase flow state. Therefore, by solving the numerical solution of the flux U at each time step, the change process of the two-phase flow state in the calculation domain can be determined, thereby predicting the movement of underwater explosion bubbles.

[0037] For underwater explosion compressible fluid, since the fluid velocity is very large during the explosion, the Re property takes a large value, the viscosity effect is weak, and the effect of the source term is also weak. Therefore, the viscosity term and the source term are ignored within the error range. Therefore, the viscosity term F in the unsteady viscous uniform two-phase flow Euler equation of formula (1) is u and G v And the source term S(U) is set to 0 and then multiplied by the smoothing function ω(r) and volume integrated over the fluid unit Ω to obtain:

[0038]

[0039] Among them, U t is the partial derivative of U with respect to time, F x (U) is the partial derivative of F(U) in the x direction, G y (U) is the partial derivative of G(U) in the y direction. ω(r) can be any smooth function, and r represents the position vector of the space point.

[0040] Further introducing the divergence theory into the flux U, the system control equation is further constructed from equation (6):

[0041] ∫ Ω ω(r)Ut dΩ+∫ Γ ω(r)[F(U)·n x +G(U)·n y ]dΓ-∫ Ω [ω(r) x F(U)+ω(r) y G(U)]dΩ=0(7)

[0042] Where Γ is the boundary of the fluid cell Ω. x is the x-component of the local unit external normal vector n of the boundary Γ, n y is the y-component of the local unit external normal vector n of the boundary Γ. ω(r) x is the x-component of the smooth function ω(r), ω(r) y is the y-component of the smooth function ω(r).

[0043] Step 2: Use K+1 basis functions to perform polynomial expansion on the flux U and smooth function ω(r) in the unsteady viscous uniform two-phase flow Euler equation.

[0044] K is an integer parameter, P 1 The K+1 basis functions in the space can be expressed as {φ0,φ1,φ2…φ K}, in one example, take 27 basis functions {φ0, φ1, φ2…φ 26}, the specific basis function can be selected according to actual needs. Then the polynomial expansion form of flux U and smooth function ω(r) can be written as:

[0045]

[0046] where, for any integer parameter k∈[0,K], φ k is the kth basis function, c k is the polynomial expansion of the flux U with the kth basis function φ k The corresponding coefficient, ω k is the polynomial expansion of the smooth function ω(r) with the kth basis function φ k The corresponding coefficients. Since the smooth function ω(r) is known, there can be multiple {ω 0 ,ω 1 ,ω 2 …ω K The coefficient combination of} can perform polynomial expansion on ω(r) based on the basis function, and generally a set of linearly independent combinations can be selected to achieve this.

[0047] Step 3: Substitute the polynomial expansion form (8) of the flux U and the smooth function ω(r) in the unsteady viscous uniform two-phase flow Euler equation into the system control equation (6), and the discrete control equations corresponding to the K+1 basis functions can be obtained. The discrete control equation corresponding to the kth basis function is:

[0048] ∫ Ω φ k c k dΩ+∫ Γ φ k h n (U - ,U + )dΓ-∫ Ω [φ k,x F(U)+φ k,y G(U)]dΩ=0 (9)

[0049] In the above formula (9), φ k,x is the basis function φ k The partial derivative in the x direction, φ k,y is the basis function φ k Partial derivative in the y direction.

[0050] Introducing the numerical flux h n (U - ,U + ) is used to characterize the flux in the normal direction of the fluid unit boundary:

[0051]

[0052] Where q = un x +vn y Represents the projection of velocity on the external normal.

[0053] Step 4: Set any k-th basis function φ k The corresponding discrete control equation is rewritten as a partial differential equation to obtain the coefficient c k The solution expression of . The discrete control equation corresponding to the k-th basis function of equation (10) is rewritten into the form of a partial differential equation:

[0054]

[0055] Where M represents the mass matrix, and any element in the mass matrix M with m1 rows and m2 columns is Integer parameter m1∈[0,K], integer parameter m2∈[0,K].

[0056] Step 5: The coefficient c shown in (11) is obtained through the above steps k The solution expression is related to the two-phase flow state. kThe solution expression of is applied to the polynomial expansion form of flux U (8), and then substituted into the semi-discrete equation of the Euler equation for unsteady viscous uniform two-phase flow.

[0057] The semi-discrete equation of the unsteady viscous uniform two-phase flow Euler equation is:

[0058]

[0059] Among them, U ij It represents the flux U at the i-th fluid unit (i, j) in the x-direction and the j-th fluid unit (i, j) in the y-direction.

[0060] F i+1 / 2,j represents the F(U) of the center point of the right boundary of the fluid unit (i, j) along the x direction, F i-1 / 2,j F(U) represents the center point of the left boundary of the fluid unit (i, j) along the x direction, Δx i Represents the distance between the center points of the two side boundaries of the fluid unit (i, j) along the x-direction.

[0061] G i,j+1 / 2 G(U) represents the position of the center point of the upper boundary of the fluid unit (i, j) along the y direction, G i,j-1 / 2 G(U), Δy, represents the position of the center point of the lower boundary of the fluid unit (i, j) along the y direction j Represents the distance between the center points of the two side boundaries of the grid (i, j) along the y direction.

[0062] Substituting into the semi-discrete equation and solving it, we can get the numerical solution of flux U in the unsteady viscous uniform two-phase flow Euler equation, and thus determine the two-phase flow state indicated by the numerical solution of flux U. The above process is used to solve at each time step, so that the numerical solution of flux U in the unsteady viscous uniform two-phase flow Euler equation at each time step can be obtained. According to the two-phase flow state indicated by the numerical solution of flux U at each time step, the movement of underwater explosion bubbles can be predicted.

[0063] In addition, after obtaining the numerical solution of the flux U at each time step, the sound velocity of the fluid at the current time step can be further calculated. And the viscosity of the fluid at the current time step μ=α1μ1+α2μ2, where μ1 is the viscosity coefficient of the gas phase and μ2 is the viscosity coefficient of the liquid phase.

[0064] In one example, the discontinuous Galilean method of the present application and the traditional finite volume method are used to calculate the propagation process of the shock wave generated by the underwater explosion to predict the bubble movement. The initial discontinuous wave is located at x=0.02m, and after 10us, the shock wave front reaches x=0.0764m. Figure 2The figure shows the comparison results of the distribution curves of fluid density calculated by the finite volume method and the discontinuous Galileo method. Figure 3 The figure shows the comparison of the distribution curves of fluid pressure calculated by the finite volume method and the discontinuous Galvan method. Figure 2 and Figure 3 It can be seen that under the same calculation grid, the wavefront calculated using the discontinuous Galerkin method of the present application is sharper, the shock wave peak is higher, and the calculation effect is better.

[0065] In another example, the intermittent Galilean method of this application and the commercial software AUTODYN are used to calculate the detonation wave generation and propagation process of two-dimensional underwater charges. The explosion process of TNT explosives in two-dimensional space under rigid wall boundary conditions is considered and compared with the results of AUTODYN-2D. The calculation domain is [-0.1, 0.1] × [-0.1, 0.1] m 2 , a 200×200 grid was used for discretization. The calculation results of fluid density obtained by AUTODYN-2D at 5μs, 10μs, 15μs, 20μs, and 25μs are shown in the figure below. Figure 4 As shown, the calculation results of the fluid density obtained by the method of the present application at 5μs, 10μs, 15μs, 20μs, and 25μs are shown in FIG. Figure 5 As shown. Figure 4 and Figure 5 It can be seen that in the distribution result diagram of fluid density, the results of AUTODYN-2D are obviously asymmetric, while the calculation results obtained in this application have better central symmetry. The curves of the calculated results of fluid density at two typical measuring points in the calculation domain over time are shown in Figure 2. Figure 6 As shown, from Figure 6 It can be seen that when the first explosion wave passes through these two typical measurement points, the calculation results of the method of the present application show a shorter rising edge and higher peak density and pressure than the results obtained by AUTODYN-2D at 5.9μs and 8.8μs. After the explosion wave is reflected from the front rigid wall, at about 23.8μs and 28.4μs, the curve of the fluid density calculated by AUTODYN-2D shows serious non-physical oscillations, while the calculation results of the intermittent Galvanometer method of the present application are significantly better. These comparison results show that the calculation accuracy and stability of the intermittent Galvanometer method of the present application are superior to AUTODYN-2D.

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

Claims

1. A method for predicting underwater explosion bubble movement based on the intermittent Galekin method, characterized in that: The underwater explosion bubble motion prediction method comprises: The viscosity term and source term in the unsteady viscous uniform two-phase flow Euler equation are set to 0 and then multiplied by the smooth function ω(r) and integrated over the fluid unit Ω. The divergence theory is introduced into the flux in the unsteady viscous uniform two-phase flow Euler equation to obtain the system control equation; where r represents the position vector of the space point; The flux U and smooth function ω(r) in the unsteady viscous uniform two-phase flow Euler equation are expanded into polynomials using K+1 basis functions: φ k is the kth basis function, c k is the polynomial expansion of the flux U with the kth basis function φ k The corresponding coefficient, ω k is the polynomial expansion of the smooth function ω(r) with the kth basis function φ k The corresponding coefficients, K is an integer parameter; Substitute the flux U and the polynomial expansion form of the smooth function ω(r) in the unsteady viscous uniform two-phase flow Euler equation into the system control equation to obtain the discrete control equation corresponding to each basis function. k The corresponding discrete control equation is rewritten as a partial differential equation to obtain the coefficient c k The solution expression of ; Each coefficient c k The solution expression is applied to the polynomial expansion form of flux U and substituted into the semi-discrete equation of unsteady viscous uniform two-phase flow Euler equation, and the numerical solution of flux U in the unsteady viscous uniform two-phase flow Euler equation is obtained; the motion of underwater explosion bubbles is predicted according to the two-phase flow state indicated by the numerical solution of the flux in the unsteady viscous uniform two-phase flow Euler equation at each time step.

2. The underwater explosion bubble motion prediction method according to claim 1, characterized in that: The Euler equation for unsteady viscous uniform two-phase flow is: in, F u is the flux of the viscosity term in the x direction, G v is the flux of the viscosity term in the y direction, α1 is the volume fraction of the gas phase, ρ1 is the density of the gas phase; α2 is the volume fraction of the liquid phase, ρ2 is the density of the liquid phase; ρ is the total density of the fluid and ρ=ρ1α1+ρ2α2, u is the x-component of the fluid velocity vector u, v is the y-component of the fluid velocity vector u, E is the internal energy per unit mass of the fluid, p is the fluid pressure, μ is the viscosity of the fluid, g is the acceleration due to gravity, and t is time; The viscosity term F in the unsteady viscous uniform two-phase flow Euler equation is u and G v And the source term S(U) is set to 0 and then multiplied by the smoothing function ω(r) and integrated over the fluid unit Ω to obtain: Among them, U t is the partial derivative of U with respect to time, F x (U) is the partial derivative of F(U) in the x direction, G y (U) is the partial derivative of G(U) in the y direction; Further introducing the divergence theory into the flux U, the system control equation is obtained as follows: ∫ Ω ω(r)U t dΩ+∫ Γ ω(r)[F(U)·n x +G(U)·n y ]dΓ-∫ Ω [ω(r) x F(U)+ω(r) y G(U)]dΩ=0 Where Γ is the boundary of the fluid unit Ω, n x is the x-component of the local unit external normal vector n of the boundary Γ, n y is the y-component of the local unit external normal vector n of the boundary Γ; ω(r) x is the x-component of the smooth function ω(r), ω(r) y is the y-component of the smooth function ω(r).

3. The underwater explosion bubble motion prediction method according to claim 2, characterized in that: Internal energy per unit mass of fluid e is the specific internal energy, u is the fluid velocity vector and u·u=u 2 +v 2 , the underwater explosion bubble motion prediction method also includes: Based on the isobaric assumption and the Stiffened state equation, the fluid pressure p = (γ-1)ρe-γB is determined; The internal energy per unit mass of fluid is obtained by sorting Among them, γ and B are the parameters of the fluid in the Stiffened state equation, and there is γ1 and B1 are the parameters of the gas phase in the Stiffened state equation, and γ2 and B2 are the parameters of the liquid phase in the Stiffened state equation.

4. The underwater explosion bubble motion prediction method according to claim 2, characterized in that: The discrete control equation corresponding to the kth basis function is obtained as follows: ∫ Ω f k c k dΩ+∫ Γ f k h n (U - ,U + )dΓ-∫ Ω [f k,x F(U)+φ k,y G(U)]dΩ=0 The flux in the normal direction of the fluid unit boundary is q=un x +vn y represents the projection of velocity on the external normal; φ k,x is the basis function φ k The partial derivative in the x direction, φ k,y is the basis function φ k Partial derivative in the y direction.

5. The underwater explosion bubble motion prediction method according to claim 4, characterized in that: The discrete control equation corresponding to the kth basis function is rewritten into the form of a partial differential equation: Where M represents the mass matrix, and any element in the mass matrix M with m1 rows and m2 columns is Integer parameter m1∈[0,K], integer parameter m2∈[0,K].

6. The underwater explosion bubble motion prediction method according to claim 2, characterized in that: The semi-discrete equation of the unsteady viscous uniform two-phase flow Euler equation is: Among them, U ij represents the flux U at the i-th fluid unit (i, j) in the x-direction and the j-th fluid unit (i, j) in the y-direction; F i+1 / 2,j represents the F(U) of the center point of the right boundary of the fluid unit (i, j) along the x direction, F i-1 / 2,j F(U) represents the center point of the left boundary of the fluid unit (i, j) along the x direction, Δx i Represents the distance between the center points of the two side boundaries of the fluid unit (i, j) along the x direction; G i,j+1 / 2 G(U) represents the position of the center point of the upper boundary of the fluid unit (i, j) along the y direction, G i,j-1 / 2 G(U), Δy, represents the position of the center point of the lower boundary of the fluid unit (i, j) along the y direction j Represents the distance between the center points of the two side boundaries of the grid (i, j) along the y direction.

Images