A method for predicting underwater explosion bubble motion based on grid-adaptive discontinuous galaxy fusion

Abstract

The present application discloses a method for predicting underwater explosion bubble motion based on grid adaptive discontinuous galvanometer, which relates to the field of underwater explosion technology. In the framework of underwater explosion bubble motion prediction based on discontinuous galvanometer, the method determines the solution resolution of each grid unit according to the underwater explosion bubble motion prediction result at each grid adaptive adjustment moment, and adjusts the grid division result based on the adjustment target of improving the solution resolution to enter the solution at the next forecast moment. The grid adaptive adjustment method can effectively improve the prediction accuracy, thereby achieving higher-precision result capture in the process of underwater explosion shock wave propagation and bubble motion, providing important support for interface capture of underwater near-field explosion process.

Description

Technical Field

[0001] The present application relates to the field of underwater explosion technology, and in particular to a method for predicting underwater explosion bubble motion based on grid-adaptive discontinuous galaxy. 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) 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.

[0006] (3) The discontinuous Galerkin method can effectively deal with discontinuous phenomena such as shock waves while maintaining numerical stability. It can accurately and effectively predict the propagation of shock waves in the early stages of underwater explosion bubble movement and the bubble movement process. However, its solution accuracy is still not ideal. Summary of the Invention

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

[0008] A method for predicting underwater explosion bubble motion based on grid-adaptive discontinuous Galiujin, the method comprising:

[0009] Based on the Euler equation of fluid and the EOS state equation, the governing equations of compressible fluid are constructed. The computational domain is discretized into multiple grid cells to obtain the grid division results at the first forecast time and determine the initial and boundary conditions of the computational domain.

[0010] Initialize the parameter t=1. For any t-th forecast time, use the discontinuous Galileo method to solve the compressible fluid governing equation based on the grid division result at the t-th forecast time to obtain the underwater explosion bubble motion forecast result at the t-th forecast time. The underwater explosion bubble motion forecast result at any t-th forecast time includes the fluid state vector of each grid cell in the grid division result at the t-th forecast time.

[0011] When the t-th forecast time does not belong to the grid adaptive adjustment time, the grid division result of the t-th forecast time is kept unchanged and used as the grid division result of the t+1-th forecast time. According to the underwater explosion bubble motion forecast result of the t-th forecast time and the grid division result of the t+1-th forecast time, the solution of the t+1-th forecast time is entered;

[0012] When the t-th forecast moment belongs to the grid adaptive adjustment moment, the solution resolution of each grid unit in the grid division result of the t-th forecast moment is determined according to the underwater explosion bubble motion forecast result at the t-th forecast moment, and the grid division result of the t-th forecast moment is adjusted based on the adjustment target of improving the solution resolution to obtain the grid division result of the t+1-th forecast moment, and the underwater explosion bubble motion forecast result at the t-th forecast moment is adjusted based on the grid division result of the t+1-th forecast moment; and the solution of the t+1-th forecast moment is entered according to the adjusted underwater explosion bubble motion forecast result at the t-th forecast moment and the obtained grid division result at the t+1-th forecast moment.

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

[0014] The present application discloses a method for predicting underwater explosion bubble motion based on grid adaptive discontinuous galvanometer. Within the framework of underwater explosion bubble motion prediction based on discontinuous galvanometer, the method dynamically adjusts the grid unit division during the prediction process with the adjustment target of improving the solution resolution. The grid adaptive method can effectively improve the prediction accuracy, thereby achieving higher-precision result capture during the underwater explosion shock wave propagation and bubble motion process, providing important support for interface capture of underwater near-field explosion processes. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0016] Figure 2 This is a flow chart of obtaining the solution resolution of all grid cells according to an embodiment of the present application.

[0017] Figure 3 Schematic diagram of the three fission directions of a grid cell in the form of a cube.

[0018] Figure 4 It is a schematic diagram of the fission and merging of grid cells in an example.

[0019] Figure 5 The figure is a comparison chart of the fluid density predicted by the grid-adaptive discontinuous GF method of the present application and the traditional fixed-grid discontinuous GF method in an example.

[0020] Figure 6 The figure is a comparison diagram of the fluid velocity predicted by using the grid-adaptive discontinuous GA-flow method of the present application and the traditional fixed-grid discontinuous GA-flow method in an example.

[0021] Figure 7 The figure is a comparison diagram of fluid pressure predicted by using the grid-adaptive discontinuous GA-flow method of the present application and the traditional fixed-grid discontinuous GA-flow method in an example.

[0022] Figure 8 The figure is a comparison chart of the specific internal energy predicted by using the grid-adaptive discontinuous Kaliujin method of the present application and the traditional fixed-grid discontinuous Kaliujin method in an example.

[0023] Figure 9 yes Figure 5 Schematic diagram of local refinement of grid cells near x=0.2.

[0024] Figure 10 yes Figure 5 Schematic diagram of local refinement of grid cells near x=0.94. DETAILED DESCRIPTION

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

[0026] This application discloses a method for predicting underwater explosion bubble motion based on grid adaptive intermittent galvanic gold. In the process of using the intermittent galvanic gold method to solve the underwater explosion bubble motion prediction method, a grid adaptive method is added to optimize the traditional intermittent galvanic gold method prediction method to improve the prediction accuracy. Please refer to Figure 1 The flow chart shown in FIG. 1 shows, the underwater explosion bubble motion prediction method includes the following steps:

[0027] Step 1: Based on the fluid Euler equation and the EOS state equation, the compressible fluid control equation is constructed.

[0028] For compressible, inviscid fluid dynamics, the flow of the fluid is governed by the Euler equations, which can be viewed as a simplified version of the Navier-Stokes equations and expressed in terms of material derivatives as follows:

[0029]

[0030] The three equations in the above formula (1) are respectively the mass conservation equation, the momentum conservation equation and the energy conservation equation. Where t represents time, ρ represents fluid density, u represents fluid velocity, p represents fluid pressure, and e represents specific internal energy (referring to the internal energy per unit mass). is the gradient operator.

[0031] For computational fluid dynamics, the above fluid Euler equation (1) is rewritten in a conservation form and expressed as (three-dimensional):

[0032]

[0033] Where u, v, and w are the components of the fluid velocity u in the x, y, and z directions, respectively. E is the total energy per unit volume, which is equivalent to the sum of the potential energy generated by the internal pressure and the kinetic energy of the gas, and thus:

[0034]

[0035] In order to construct a relatively generalized numerical discretization format for a class of nonlinear conservation laws, the above Euler equation can be written in the following vector form:

[0036]

[0037] where Q is the conserved state vector and F(Q), G(Q) and H(Q) are the nonlinear flux of the fluid and

[0038] The system is closed by the EOS state equation. This application adopts the following enhanced state equation:

[0039] p=(γ-1)ρe+γB (5)

[0040] Among them, γ and B are material constants. For monatomic gas, γ = 1.4 and B = 0. The speed of sound can be calculated by

[0041] The fluid Euler equation shown in formula (4) is combined with the EOS state equation shown in formula (5) to construct the compressible fluid control equation.

[0042] Step 2: Discretely divide the computational domain into multiple grid cells to obtain the grid division results for the first forecast moment and determine the initial conditions and boundary conditions of the computational domain. When performing grid division, the grid cells obtained can be of any shape, which provides flexibility for handling complex geometries. For example, it is common to divide the computational domain of a two-dimensional plane into triangles or quadrilaterals, and the grid cells obtained by dividing the computational domain of a three-dimensional space into hexahedrons or tetrahedrons. In practical applications, the computational domain of underwater explosion bubble motion forecasting scenarios is mostly three-dimensional space, and hexahedral grid cells are generally used for division.

[0043] Step 3, initialize the parameter t=1. For any t-th forecast time, use the discontinuous Galileo method to solve the compressible fluid control equation based on the grid division result at the t-th forecast time to obtain the underwater explosion bubble motion forecast result at the t-th forecast time.

[0044] When using the discontinuous Galerkin method for solving problems, a set of basis functions is selected within each grid cell to represent the local approximation of the solution. The basis functions are usually polynomials, and their order determines the accuracy of the method. Different grid cells can use basis functions of different orders to achieve local mesh refinement. The constructed compressible fluid control equation is then multiplied by a test function on both sides and integrated over the entire computational domain, thereby converting the partial differential equation into a weak form. This can reduce the smoothness requirements of the solution and allow discontinuous solutions to be properly handled. For the interfaces between adjacent grid cells, numerical fluxes need to be defined to handle the solutions at the interfaces. Numerical fluxes are constructed based on physical conservation laws and are used to approximate the true flux values ​​to ensure the continuity and conservation of the solution at the cell boundaries. For transient problems, an appropriate time integration scheme needs to be selected to advance the time step. Commonly used methods include explicit and implicit time integration schemes, and the specific choice depends on the characteristics of the problem and the required computational efficiency. Ultimately, the above steps will lead to a system of linear or nonlinear algebraic equations with unknown solution coefficients. For linear systems, an iterative or direct solver can be used directly; for nonlinear systems, an iterative method (such as the Newton-Raphson method) is usually used to gradually solve the problem, and then the underwater explosion bubble motion forecast result at the t-th forecast time is obtained. This application does not elaborate on the specific solution method of the discontinuous Jialiujin method.

[0045] The underwater explosion bubble motion prediction result obtained at any t-th prediction moment includes the fluid state vector of each grid cell in the grid division result at the t-th prediction moment, that is, the conservation state vector Q in the above formula (4). The fluid state vector can be written as x represents the grid cell Coordinates within, φ n ( ) is the grid unit The nth basis function, χ n (t) is the basis function φ n ( ) Modal expansion coefficient at the tth forecast time, N p is a grid cell The total number of basis functions included. According to the solution characteristics of the discontinuous Galerkin method, the fluid state vector of each grid cell is continuous, but discontinuities are allowed between two grid cells, so that the conserved state vector Q in the entire computational domain is a sheet-like continuous function in the finite element space.

[0046] Coordinates represent grid cells in basis functions The natural coordinates of the grid cell in the two-dimensional plane are x = (ξ, η), and the natural coordinates of the grid cell in the three-dimensional space are x = (ξ, η, ζ). This application uses normalized Jacobi polynomials to construct basis functions, and takes the ξ component of the natural coordinate x as an example for explanation. The same applies to η and ζ:

[0047] The initial n-th order Jacobian polynomial can be calculated using the Rodrigues formula:

[0048]

[0049] Among them, parameter α>-1, parameter β>-1. Jacobi polynomials satisfy the following orthogonality conditions:

[0050]

[0051] Among them, ω (α，β) (ξ)=(1-ξ) α (1+ξ) β is the weight function of the orthogonal polynomial, and m is the order parameter. When m=n, δ mn =1, otherwise δ mn =0. is a constant:

[0052]

[0053] Γ() is the quadrature function. Γ(m+α+1) represents the quadrature from 1 to m+α+1, and the same applies to other quadrature methods.

[0054] In the one-dimensional case, the orthogonal basis functions are constructed by using the normalized Jacobian polynomials:

[0055]

[0056] Among them, the normalized Jacobi polynomial The derivative of the normalized Jacobi polynomial can be expressed as:

[0057]

[0058] Combined with the orthogonality conditions of the Jacobi polynomials, the orthogonality properties of the basis functions can be obtained as follows:

[0059]

[0060] Considering that for one-dimensional isoparametric units, ξ = [-1, 1] is generally taken, so the above equation shows that the basis function is orthogonal on the isoparametric unit. For simplicity, α = β = 0, so ω (α，β) (ξ)=1 and The Jacobi polynomials simplify to Langende polynomials:

[0061]

[0062] Then we get the derivative of the basis function:

[0063]

[0064] In quadrilateral and hexahedral mesh cells, the basis functions can be constructed directly through tensor products. For quadrilateral mesh cells in a two-dimensional plane, the basis functions are:

[0065]

[0066] For a hexahedral mesh element in three-dimensional space, the basis function is:

[0067]

[0068] Among them, 0≤i, j, k≤N represents the highest order of the normalized Jacobi polynomial, and N determines the convergence rate of the numerical algorithm. The above basis functions are orthogonal on the isoparametric unit. The isoparametric quadrilateral unit is a square in the two-dimensional plane and has a natural coordinate space. The isoparametric hexahedral element is a cube in three-dimensional space and The derivatives of the basis functions can be found using the chain rule.

[0069] The construction of orthogonal basis functions on triangles and tetrahedrons is slightly more complicated. The Prori-Koornwinder-Dubiner (PKD) polynomials constructed from Jacobi polynomials are used as basis functions, which are expressed as follows:

[0070] For a triangular mesh element in a two-dimensional plane, the basis function is:

[0071]

[0072] For a tetrahedral mesh element in three-dimensional space, the basis function is:

[0073]

[0074] The above basis functions are orthogonal on the isoparametric unit. The isoparametric triangular unit is a triangle in two dimensions and the natural coordinate space is The isoparametric tetrahedral element is a tetrahedron in three-dimensional space and its natural coordinate space is Similarly, the derivatives of the basis functions can be found using the chain rule.

[0075] Step 4: When the t-th forecast moment does not belong to the grid adaptive adjustment moment, the grid division result of the t-th forecast moment is kept unchanged and used as the grid division result of the t+1-th forecast moment. According to the underwater explosion bubble motion forecast result at the t-th forecast moment and the grid division result at the t+1-th forecast moment, the solution of the t+1-th forecast moment is entered.

[0076] In the traditional underwater explosion bubble motion prediction process using intermittent galaxy, after completing the grid division in the initial state, the grid unit division result is always kept unchanged to perform iterative solution prediction. The grid division result will affect the solution prediction process and the prediction accuracy. However, the present application adds a grid adaptive adjustment moment in the underwater explosion bubble motion prediction process. When the t-th prediction moment does not belong to the grid adaptive adjustment moment, the next prediction moment is directly entered for iterative solution as in the conventional practice. When the t-th prediction moment belongs to the grid adaptive adjustment moment, the grid division result will be adaptively adjusted. There is at least one grid adaptive adjustment moment in the entire underwater explosion bubble motion prediction process, that is, at least one grid adaptive adjustment is performed according to the following step 5.

[0077] Step 5: When the t-th forecast moment belongs to the grid adaptive adjustment moment, the solution resolution of each grid unit in the grid division result of the t-th forecast moment is determined according to the underwater explosion bubble motion forecast result at the t-th forecast moment, and the grid division result of the t-th forecast moment is adjusted based on the adjustment goal of improving the solution resolution to obtain the grid division result of the t+1-th forecast moment, and the underwater explosion bubble motion forecast result at the t-th forecast moment is adjusted based on the grid division result at the t+1-th forecast moment.

[0078] Determining the solution resolution of each grid cell in the grid division result at the t-th forecast time includes the following steps, please refer to Figure 2 The flowchart shown:

[0079] (1) Traverse each grid cell in the grid division result of the t-th forecast time in turn.

[0080] (2) For any k-th grid cell traversed, the solution resolution λ of the k-th grid cell is updated according to the fluid density in the fluid state vector of the k-th grid cell. k ,include:

[0081] (2.a) Calculate the density gradient of the fluid density ρ at different coordinates in the fluid state vector of the k-th grid cell in each fission direction of the k-th grid cell, and calculate the solution resolution in the corresponding direction based on the density gradient at different coordinates in each fission direction.

[0082] The fission directions of the grid cells can be determined according to the shape of the initially divided grid cells. For example, the common cubic grid cells contain three mutually perpendicular fission directions η1, η2, and η3. Calculate the density gradients ρ1, ρ2, and ρ3 of the fluid density ρ at each coordinate in the fluid state vector of the kth grid cell in the fission directions η1, η2, and η3. Similarly, the density gradients of the fluid density ρ at other coordinates in the fission directions η1, η2, and η3 can be obtained, thereby determining the density gradient ranges of different coordinates in the kth grid cell in each fission direction.

[0083] The density gradient ranges of all mesh cells in the entire computational domain along each fission direction are then calculated. The relative range of the density gradient range of the kth mesh cell in each fission direction relative to the density gradient range of the entire computational domain along the current fission direction is then compared to determine the solution resolution of the kth mesh cell in the current fission direction. The larger the relative range of the density gradient range of the kth mesh cell in each fission direction relative to the density gradient range of the entire computational domain along the current fission direction, the higher the corresponding solution resolution. The specific correspondence can be customized. For example, if the density gradient range of the entire computational domain along fission direction η1 is [1,100], and the density gradient range of the kth mesh cell along fission direction η1 is [40,50], the solution resolution along η1 is determined to be 1. If the density gradient range of the kth mesh cell along fission direction η1 is [80,90], the solution resolution along η1 is determined to be 2. The specific correspondence can be customized. In this way, the solution resolutions Θ1, Θ2, and Θ3 of the kth grid unit in the fission directions η1, η2, and η3 can be obtained respectively.

[0084] (2.b) Calculate the maximum solution resolution of the kth grid cell in each fission direction as the overall resolution of the kth grid cell. For example, in the above example, the maximum solution resolution Θ1, Θ2, and Θ3 is taken as the overall resolution of the kth grid cell.

[0085] It should be noted that the fluid state vector includes other parameters in addition to the fluid density. This application mainly considers the density gradient to obtain the overall resolution of the kth grid unit because the grid adaptive adjustment effect achieved in this way is better. In fact, the gradients of other parameters can also be calculated to obtain the overall resolution of the kth grid unit, such as calculating the gradient of the fluid pressure to obtain the solution resolution, etc., which are all calculated in the same way. This application does not limit this.

[0086] (2.c) If the kth grid cell has not yet been solved, the calculated overall resolution is used as the solution resolution for the kth grid cell. If the kth grid cell has already been solved, the larger of the calculated overall resolution and the existing solution resolution is used as the solution resolution for the kth grid cell.

[0087] (3) According to the relative position relationship between other grid cells and the k-th grid cell, the solution resolution λ of the k-th grid cell is used k Update the solution resolution of each other grid cell. From this step, we can see that when traversing the k-th grid cell, not only the solution resolution of the k-th grid cell itself is calculated and updated, but the solution resolution of other grid cells is also updated. Similarly, when traversing other grid cells, the solution resolution of the k-th grid cell may also be updated. This is also the reason why the k-th grid cell may already have a solution resolution when traversing the k-th grid cell in (2.c) above.

[0088] Determine that the grid cells that share a common node with the k-th grid cell and are directly adjacent to the k-th grid cell are the first-level adjacent cells of the k-th grid cell. For example, when the grid cell is a cubic structure, there are 6 grid cells in front, back, left, right, top and bottom that share a common node with the k-th grid cell in the three fission directions η1, η2, and η3. These 6 grid cells are the first-level adjacent cells of the k-th grid cell. Continuing to extrapolate based on this principle, for any integer parameter 1≤g≤λ k , determine that the other directly adjacent grid cells that share the same nodes with the g-th level neighboring cells of the k-th grid cell are the g+1-th level neighboring cells of the k-th grid cell, and then k -g updates the solution resolution of the g-th level neighboring unit of the k-th grid unit. This includes: for any g-th level neighboring unit of the k-th grid unit, when the solution resolution of the g-th level neighboring unit has not been calculated, k -g is the solution resolution of the g-th level neighboring unit. When the g-th level neighboring unit has a solution resolution, k The larger value of the solution resolution of -g and the g-th level neighboring unit is used as the updated solution resolution of the g-th level neighboring unit.

[0089] (4) If all grid cells in the grid division result at the t-th forecast time have not been traversed, continue traversing the next grid cell that has not been traversed and repeat the above process. After all grid cells in the grid division result at the t-th forecast time have been traversed, the solution resolution of all grid cells is finally obtained.

[0090] After obtaining the solution resolution of each grid cell in the grid division result at the t-th forecast time, it is compared with the solution resolution of the grid cell at the last grid adaptive adjustment time and grid adaptive adjustment is performed, wherein the first forecast time is used as the last grid adaptive adjustment time of the first grid adaptive adjustment time. For any grid cell in the grid division result at the t-th forecast time, the comparison result of its solution resolution with the grid cell at the last grid adaptive adjustment time has three cases:

[0091] (1) Case 1: If the solution resolution of the grid cell at the current grid adaptive adjustment moment is higher than the solution resolution of the grid cell in the corresponding area at the last grid adaptive adjustment moment, it means that the solution resolution of the area where the grid cell is located has been improved. In this case, the grid cell in the grid division result at the t-th forecast moment is further split into multiple grid cells to improve the accuracy.

[0092] In one embodiment, when a grid unit is further fissioned into multiple grid units, the grid unit is equally divided into several levels along each fission direction to obtain multiple grid units; the more the solution resolution of the grid unit is improved compared to the solution resolution of the grid unit in the corresponding area at the last grid adaptive adjustment moment, the more levels the grid unit is divided into along each fission direction.

[0093] In addition, when splitting a grid cell into multiple grid cells, it is also necessary to correspondingly adjust the underwater explosion bubble motion forecast result at the grid cell at the t-th forecast time to obtain the fluid state vector of each grid cell after the fission, so that the adjusted underwater explosion bubble motion forecast result at the t-th forecast time matches the grid division result at the t+1-th forecast time. This includes: assigning the fluid state vector of the grid cell before the fission contained in the underwater explosion bubble motion forecast at the t-th forecast time to each of the fissioned grid cells.

[0094] For example, in an example, taking the grid unit as a cube, one case is to use Figure 4 The grid unit shown in (a) is divided into two equal levels in each fission direction, thereby obtaining 8 grid units as shown in FIG. Figure 4 Or in another case, Figure 4 The grid cells shown in (a) are divided into 4 levels in each fission direction, thereby obtaining 64 grid cells as shown in FIG. Figure 4 As shown in (c) in .

[0095] (2) Case 2: If the solution resolution of the grid cell at the current grid adaptive adjustment moment is lower than the solution resolution of the grid cell in the corresponding area at the previous grid adaptive adjustment moment, it means that the solution resolution of the area where the grid cell is located has decreased. In this case, the grid cell in the grid division result at the t-th forecast moment and several adjacent grid cells are merged into one grid cell.

[0096] In one embodiment, when merging multiple grid cells in the grid division result at the t-th forecast time into one grid cell, the grid cell and multiple adjacent grid cells along each fission direction are merged into one grid cell. The greater the degree of decrease in the solution resolution of the grid cell compared to the solution resolution of the grid cell in the corresponding area at the last grid adaptive adjustment time, the greater the number of grid cells merged for the grid cell along each fission direction.

[0097] Similarly, when merging multiple grid cells into one grid cell, it is also necessary to correspondingly adjust the underwater explosion bubble motion forecast results of these multiple grid cells at the t-th forecast time to obtain the fluid state vector of the merged grid cell, so that the adjusted underwater explosion bubble motion forecast result at the t-th forecast time matches the grid division result at the t+1-th forecast time. This includes assigning the average value of the fluid state vectors of the multiple grid cells before merging, which are included in the underwater explosion bubble motion forecast at the t-th forecast time, to the merged grid cell.

[0098] For example, in one embodiment, a grid cell and its adjacent grid cell along each fission direction are merged into one grid cell, thereby Figure 4 The eight grid cells shown in (b) are merged into one grid cell. Figure 4 As shown in (d) in another example, the grid unit and its three adjacent grid units along each fission direction are merged into one grid unit, so that Figure 4 The 64 grid cells shown in (c) are merged into one grid cell as shown in Figure 4 As shown in (f) in FIG. In another example, a grid unit and its adjacent grid unit along each fission direction are merged into one grid unit, thereby Figure 4 The 64 grid cells shown in (c) are merged into 8 grid cells as shown in Figure 4 As shown in (g)

[0099] (3) When the solution resolution of the grid cell at the current grid adaptive adjustment moment is consistent with the solution resolution of the grid cell in the corresponding area at the last grid adaptive adjustment moment, it means that the solution resolution of the area where the grid cell is located is basically unchanged, and the grid cell is kept unchanged.

[0100] Then, according to the adjusted underwater explosion bubble movement forecast result at the t-th forecast moment and the obtained grid division result at the t+1-th forecast moment, the solution for the t+1-th forecast moment is entered. At this time, the grid division result at the t+1-th forecast moment is often different from the grid division result at the t-th forecast moment, and is more conducive to improving the forecast accuracy. In the actual forecast process, a grid adaptive adjustment moment is set every several forecast moments in the underwater explosion bubble movement forecast process, and the grid division result is dynamically adjusted at each grid adaptive adjustment moment, that is, after the grid units are merged at a grid adaptive adjustment moment, they may continue to be merged or may fission again. After the grid units are fissioned at a grid adaptive adjustment moment, they may continue to fission or may be merged again, and are dynamically adjusted according to the actual solution resolution of the forecast process. For example, Figure 4 In the example shown, Figure 4 After the grid unit shown in (a) in the figure is split into 8 grid units shown in (b) at one grid adaptive adjustment moment, it may continue to split into 64 grid units shown in (e) at the next grid adaptive adjustment moment, or it may be re-merged into one grid unit shown in (d) at the next grid adaptive adjustment moment.

[0101] In a numerical test of a one-dimensional SOD discontinuity problem, the initial conditions are set as follows: the flow field density of the left pipe is 1, the flow field velocity is 0, and the flow field pressure is 1; the flow field density of the right pipe is 0.125, the flow field velocity is 0, and the flow field pressure is 0.1. The fluid is an ideal gas, γ = 1.25, the computational domain is [0,1], the initial grid number is 200, and the calculation is performed using the fixed grid discontinuity Galin method and the grid adaptive discontinuity Galin method of this application. The comparison results of the flow field density, flow field velocity, flow field pressure and specific internal energy distribution curves at time T = 0.25 are shown in the figure below. Figure 5-Figure 8 As shown, the red curve is the forecast result obtained by using the grid adaptive discontinuous Jialiujin method of this application, and the black curve is the forecast result obtained by using the fixed grid discontinuous Jialiujin method. Figure 5-Figure 8 It can be seen that the prediction results obtained by using the grid-adaptive discontinuous Kaliujin method have significantly improved the capture accuracy of both density discontinuities and velocity and pressure discontinuities, reducing numerical dispersion without significantly increasing the amount of calculation. In the process of using the grid-adaptive discontinuous Kaliujin method of this application for forecasting, the total number of grids at T=0.25 reached 241, an increase of 20.5% compared to the initial 200 grids, and the increased grids are mainly distributed in Figure 5-Figure 8 The density gradient on the left side is large, mainly concentrated near x = 0.2 and x = 0.94. The local density effect at these two locations is as follows: Figure 9 Shown and Figure 10 shown.

[0102] 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 motion based on grid-adaptive discontinuous galaxy, characterized in that: The underwater explosion bubble motion prediction method comprises: Based on the Euler equation of fluid and the EOS state equation, the governing equations of compressible fluid are constructed. The computational domain is discretized into multiple grid cells to obtain the grid division results at the first forecast time and determine the initial and boundary conditions of the computational domain. Initializing parameter t=1, for any t-th forecast time, using the discontinuous Galileo method to solve the compressible fluid governing equation based on the grid division result at the t-th forecast time, obtain the underwater explosion bubble motion forecast result at the t-th forecast time; wherein the underwater explosion bubble motion forecast result at any t-th forecast time includes the fluid state vector of each grid cell in the grid division result at the t-th forecast time; When the t-th forecast time does not belong to the grid adaptive adjustment time, the grid division result of the t-th forecast time is kept unchanged and used as the grid division result of the t+1-th forecast time. According to the underwater explosion bubble motion forecast result of the t-th forecast time and the grid division result of the t+1-th forecast time, the solution of the t+1-th forecast time is entered; When the t-th forecast moment belongs to the grid adaptive adjustment moment, the solution resolution of each grid unit in the grid division result of the t-th forecast moment is determined according to the underwater explosion bubble motion forecast result at the t-th forecast moment, and the grid division result of the t-th forecast moment is adjusted based on the adjustment target of improving the solution resolution to obtain the grid division result of the t+1-th forecast moment, and the underwater explosion bubble motion forecast result at the t-th forecast moment is adjusted based on the grid division result of the t+1-th forecast moment; and the solution of the t+1-th forecast moment is entered according to the adjusted underwater explosion bubble motion forecast result at the t-th forecast moment and the obtained grid division result at the t+1-th forecast moment.

2. The underwater explosion bubble motion prediction method according to claim 1, wherein The grid division result at the t-th forecast moment is adjusted based on the adjustment target of improving the solution resolution to obtain the grid division result at the t+1-th forecast moment, including any grid unit in the grid division result at the t-th forecast moment: When the solution resolution of the grid unit at the current grid adaptive adjustment moment is higher than the solution resolution of the grid unit in the corresponding area at the previous grid adaptive adjustment moment, the grid unit in the grid division result at the t-th forecast moment is split into multiple grid units; wherein the first forecast moment is used as the previous grid adaptive adjustment moment of the first grid adaptive adjustment moment; When the solution resolution of the grid cell at the current grid adaptive adjustment moment is lower than the solution resolution of the grid cell in the corresponding area at the previous grid adaptive adjustment moment, the grid cell and several adjacent grid cells in the grid division result at the t-th forecast moment are merged into one grid cell; When the solution resolution of the grid unit at the current grid adaptive adjustment moment is consistent with the solution resolution of the grid unit in the corresponding area at the last grid adaptive adjustment moment, the grid unit is kept unchanged.

3. The underwater explosion bubble motion prediction method according to claim 2, wherein: Adjusting the underwater explosion bubble motion prediction result at the t-th prediction moment based on the grid division result at the t+1-th prediction moment includes: When a grid cell in the grid division result at the t-th forecast time fissions into multiple grid cells in the grid division result at the t+1-th forecast time, the fluid state vector of the grid cell before fission contained in the underwater explosion bubble motion forecast at the t-th forecast time is assigned to each of the fissioned grid cells; When multiple grid cells in the grid division result at the t-th forecast moment are merged into one grid cell in the grid division result at the t+1-th forecast moment, the average value of the fluid state vectors of the multiple grid cells before the merger contained in the underwater explosion bubble motion forecast at the t-th forecast moment is assigned to the merged grid cell.

4. The underwater explosion bubble motion prediction method according to claim 2, wherein: When a grid unit in the grid division result at the t-th forecast moment is fissioned into multiple grid units, the grid unit is equally divided into several levels along each fission direction to obtain multiple grid units; the more the solution resolution of the grid unit is improved compared with the solution resolution of the grid unit in the corresponding area at the last grid adaptive adjustment moment, the more levels the grid unit is divided into along each fission direction.

5. The underwater explosion bubble motion prediction method according to claim 2, wherein: When merging multiple grid cells in the grid division result of the t-th forecast moment into one grid cell, the grid cell and multiple adjacent grid cells along each fission direction are merged into one grid cell; the more the solution resolution of the grid cell decreases compared to the solution resolution of the grid cell in the corresponding area at the last grid adaptive adjustment moment, the more grid cells are merged for the grid cell along each fission direction.

6. The underwater explosion bubble motion prediction method according to claim 1, wherein: The solution resolution of each grid cell in the grid division result at the t-th forecast time is determined based on the underwater explosion bubble motion forecast result at the t-th forecast time, including: Traverse each grid cell in the grid division result of the t-th forecast time in turn. For any k-th grid cell traversed, update the solution resolution of the k-th grid cell according to the fluid density in the fluid state vector of the k-th grid cell. , and based on the relative position relationship between other grid cells and the k-th grid cell, the solution resolution of the k-th grid cell is used Update the solution resolution of other grid cells until all grid cells in the grid division result of the t-th forecast time are traversed to obtain the solution resolution of all grid cells.

7. The underwater explosion bubble motion prediction method according to claim 6, wherein: The solution resolution of the kth grid cell is obtained by updating the fluid density in the fluid state vector of the kth grid cell. include: Calculate the density gradients of the fluid density at different coordinates in the fluid state vector of the kth grid unit in each fission direction of the kth grid unit, and calculate the solution resolution in the corresponding direction according to the density gradients at different coordinates in each fission direction; Calculate the maximum value of the solution resolution of the k-th grid cell in each fission direction as the overall resolution of the k-th grid cell; When the solution resolution of the k-th grid unit has not been calculated, the calculated overall resolution is used as the solution resolution of the k-th grid unit; when the solution resolution of the k-th grid unit has been calculated, the larger value of the calculated overall resolution and the existing solution resolution is used as the solution resolution of the k-th grid unit.

8. The underwater explosion bubble motion prediction method according to claim 6, wherein: The solution resolution of the kth grid cell is used according to the relative position relationship between other grid cells and the kth grid cell. Updating the solution resolution of other mesh cells includes: Determine the grid cells that share the same node with the k-th grid cell and are directly adjacent to it as the first-level adjacent cells of the k-th grid cell. For any integer parameter , determine that the other grid cells directly adjacent to the g-th level neighboring cells of the k-th grid cell are the g+1-th level neighboring cells of the k-th grid cell, according to Updates the solution resolution of the g-th neighboring cells of the k-th grid cell.

9. The underwater explosion bubble motion prediction method according to claim 8, wherein: The basis Updating the solution resolution of the g-th level neighboring cells of the k-th grid cell includes: For any g-th level neighboring unit of the k-th grid unit, when the resolution of the g-th level neighboring unit has not been calculated, As the solution resolution of the g-th level adjacent unit; when the g-th level adjacent unit has a solution resolution, and the solution resolution of the g-th level adjacent unit as the updated solution resolution of the g-th level adjacent unit.

10. The underwater explosion bubble motion prediction method according to claim 2, characterized in that: During the underwater explosion bubble movement prediction process, a grid adaptive adjustment moment is set every several prediction moments, and the grid division result is dynamically adjusted at each grid adaptive adjustment moment.

Images