A method for predicting underwater explosion bubble motion based on the discontinuous Galukin method

Abstract

The present application discloses a method for predicting underwater explosion bubble motion based on the discontinuous Galerkin method, which relates to the field of explosion technology. The method comprises: setting the viscous term and the source term in the unsteady viscous uniform two-phase flow Euler equation to 0, multiplying them by a smooth function, and integrating them over the fluid unit, introducing the divergence theory to construct a system control equation, and then using basis functions to perform polynomial expansion on the flux and the smooth function respectively, and then inserting the latter into the system control equation to obtain the discrete control equation corresponding to each basis function, and then obtaining the polynomial expansion form of the flux, and then inserting the latter into the semi-discrete equation of the unsteady viscous uniform two-phase flow Euler equation to solve and obtain the two-phase flow state at each time step, thereby predicting the propagation of the shock wave in the early stage of underwater explosion bubble motion and the bubble motion process, realizing the interface capture of the underwater near-field explosion process, and providing solid theoretical support for the in-depth study of the 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 discontinuous Galiuge method. Background Art

[0002] As ocean 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 explosions, the prediction of the propagation of the shock wave and the bubble movement process 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). 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 based on the BEM method. 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 structural interfaces, the BEM method requires cutting and stitching the grid, which reduces the solution accuracy and makes the algorithm complex.

[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 Euler grid method 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 underwater explosion bubble pulsation by correcting the non-reflecting boundary. However, due to its Lagrangian characteristics, the meshless SPH method poses great challenges in the simulation of strongly compressible problems, because the distance between the initially distributed explosive particles will increase sharply after 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 underwater explosion bubble motion, and are unable to accurately and effectively predict the propagation of shock waves in the early stage of underwater explosion bubble motion and the bubble motion process. Summary of the Invention

[0008] In response to the above-mentioned problems and technical needs, this application proposes a method for predicting underwater explosion bubble motion based on the discontinuous Galiuge method. The technical solution of this application is as follows:

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

[0010] The viscous and source terms in the unsteady viscous uniform two-phase flow Euler equation are set to zero, multiplied by the smoothing function ω(r), and then integrated over the fluid unit Ω. The divergence theory is introduced into the flux in the unsteady viscous uniform two-phase flow Euler equation to construct the system governing equation; where r represents the position vector of the spatial point;

[0011] The flux U and smooth function ω(r) in the Euler equation for unsteady viscous uniform two-phase flow 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 form of the smooth function ω(r) and the k-th basis function φ k The corresponding coefficient, 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 the flux U and substituted into the semi-discrete equation of the unsteady viscous uniform two-phase flow Euler equation to obtain the numerical solution of the flux U in the unsteady viscous uniform two-phase flow Euler equation; the motion of underwater explosion bubbles is predicted based on 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 Gallicin method. The method is based on the Euler equation for unsteady viscous uniform two-phase flow and introduces the discontinuous Gallicin 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 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 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 of the fluid pressure distribution curves 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 this application are compared. DETAILED DESCRIPTION

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

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

[0024] Step 1: Based on the Euler equation of unsteady viscous uniform two-phase flow and the divergence theory, the system control equation is constructed.

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

[0026]

[0027] in, F u is the flux in the x direction of the viscous term, 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 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 units in this dissipative region is considered to be a mixture of gas and liquid phases, and the volume fractions of the gas and liquid phases in the fluid satisfy α1+α2=1. ρ is the total density of the fluid and ρ=ρ1α1+ρ2α2. u is the x-component of the fluid velocity vector u, and v is the y-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, and γ2 = 4.4 and B2 = 6×10 8 Pa.

[0034] Substituting equations (2), (3), and (4) above 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 as follows by the fluid pressure p, the fluid velocity vector u, and the volume fractions of the gas and liquid phases:

[0035]

[0036] From the above introduction, it can be seen that the flux U in the Euler equation for unsteady viscous uniform two-phase flow contains the parameters of the two-phase flow state. Therefore, solving the numerical solution of flux U at each time step can determine the change process of the two-phase flow state in the computational domain, 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 is large, the viscosity effect is weak, and the source term is also weak. Therefore, the viscosity term and 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 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 spatial point.

[0040] Further introducing the divergence theory for 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 unit Ω. 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-direction component of the smooth function ω(r), ω(r) y is the y-direction component of the smooth function ω(r).

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

[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] 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 form of the smooth function ω(r) and the k-th 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 smoothing 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 Partial derivative in the x direction, φ k,y is the basis function φ k Partial derivative in the y direction.

[0050] Introducing 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, and the coefficients c kThe solution expression is applied to the polynomial expansion form of the 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 Euler equation for unsteady viscous uniform two-phase flow 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 F(U) represents the center point of the right boundary of the fluid unit (i, j) along the x direction, 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 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 sides of the grid (i, j) along the y direction.

[0062] Substituting this into the semi-discrete equation and solving it yields a numerical solution to the flux U in the Euler equation for unsteady viscous uniform two-phase flow, thereby determining the two-phase flow state indicated by the numerical solution to flux U. By performing the above process at each time step, the numerical solution to flux U in the Euler equation for unsteady viscous uniform two-phase flow at each time step can be obtained. Based on the two-phase flow state indicated by the numerical solution to flux U at each time step, the motion 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 were used to calculate the propagation process of the shock wave generated by an underwater explosion to predict the bubble movement. The initial discontinuous wave was located at x = 0.02m, and after 10us, the shock wave front reached 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 galvanometer method. Figure 3 The figure shows the comparison of the distribution curves of fluid pressure calculated by the finite volume method and the discontinuous Gauss 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 discontinuous Gallantry method of this application and the commercial software AUTODYN were used to calculate the detonation wave generation and propagation process of a two-dimensional underwater charge. The explosion process of TNT explosives in two-dimensional space under rigid wall boundary conditions was considered and compared with the results of AUTODYN-2D. The calculation domain is [-0.1, 0.1] × [-0.1, 0.1] m 2 , using a 200×200 grid for discretization. The calculation results of the 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 in the figure, the calculation results of the fluid density obtained by the method of this application at 5μs, 10μs, 15μs, 20μs, and 25μs are as follows Figure 5 As shown. Figure 4 and Figure 5 It can be seen that the results of AUTODYN-2D are obviously asymmetric in the distribution of fluid density, while the results obtained by this application have better central symmetry. The curves of the calculated results of fluid density at two typical measurement points in the calculation domain over time are as follows: Figure 6 As shown, from Figure 6 As can be seen, when the first blast wave front passes through these two typical measurement points, the calculation results of the present application's method show a shorter rise time and higher peak density and pressure than the results obtained by AUTODYN-2D at 5.9μs and 8.8μs. After the blast wave reflects from the front rigid wall, the fluid density curve calculated by AUTODYN-2D shows severe unphysical oscillations at approximately 23.8μs and 28.4μs, while the calculation results of the present application's discontinuous Galerkin method are significantly better. These comparative results show that the present application's discontinuous Galerkin method outperforms AUTODYN-2D in terms of calculation accuracy and stability.

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

Claims

1. A method for predicting underwater explosion bubble movement based on the discontinuous Galiuge method, characterized in that: The underwater explosion bubble motion prediction method comprises: The viscous and source terms in the unsteady viscous uniform two-phase flow Euler equation are set to zero, multiplied by the smoothing function ω(r), and then integrated over the fluid unit Ω. The divergence theory is introduced into the flux in the unsteady viscous uniform two-phase flow Euler equation to construct the system governing equation; where r represents the position vector of the spatial point; The flux U and smooth function ω(r) in the Euler equation for unsteady viscous uniform two-phase flow 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 form of the smooth function ω(r) and the k-th basis function φ k The corresponding coefficient, 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 the flux U and substituted into the semi-discrete equation of the unsteady viscous uniform two-phase flow Euler equation to obtain the numerical solution of the flux U in the unsteady viscous uniform two-phase flow Euler equation; the motion of underwater explosion bubbles is predicted based on 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, wherein The Euler equation for unsteady viscous uniform two-phase flow is: in, F u is the flux of the viscous term in the x direction, G v is the flux of the viscous 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 Euler equation for unsteady viscous uniform two-phase flow is converted to u and G v And the source term S(U) is set to 0 and 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-direction component of the smooth function ω(r), ω(r) y is the y-direction component of the smooth function ω(r).

3. The underwater explosion bubble motion prediction method according to claim 2, wherein: 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 further 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, wherein: 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 Among them, 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 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, wherein: Rewrite the discrete control equation corresponding to the kth basis function 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, wherein: The semi-discrete equation of the Euler equation for unsteady viscous uniform two-phase flow 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 F(U) represents the center point of the right boundary of the fluid unit (i, j) along the x direction, 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 sides of the fluid unit (i, j) along the x-direction; G i,j+1 / 2 G(U) represents 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 sides of the grid (i, j) along the y direction.

Images