Underwater explosion bubble movement forecasting method based on grid self-adaptive discontinuous custard flow metal

Abstract

The invention discloses an underwater explosion bubble movement forecasting method based on grid self-adaption intermittent custard flow, and relates to the technical field of underwater explosion. The solution resolution of each grid unit is determined at each grid adaptive adjustment moment according to the underwater explosion bubble motion prediction result, the grid division result is adjusted based on the adjustment target for improving the solution resolution, and the solution of the next prediction moment is started, so that the prediction precision can be effectively improved through the grid adaptive adjustment method. Therefore, higher-precision result capture can be achieved in the underwater explosion shock wave propagation and bubble movement process, and important support is provided for interface capture in the 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 movement based on grid adaptive discontinuous galaxy. 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) 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.

[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 motion and the bubble motion process. However, its solution accuracy is still not ideal. Summary of the invention

[0007] In view of the above problems and technical requirements, 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 galaxy, the method comprising:

[0009] Based on the fluid Euler equation combined with the EOS state equation, the compressible fluid control equation is constructed, the computational domain is discretized into multiple grid cells to obtain the grid division result of the first forecast time and determine the initial conditions and boundary conditions of the computational domain;

[0010] Initialize the parameter t=1, for any t-th forecast time, use the discontinuous Galician 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; wherein, the underwater explosion bubble motion forecast result at any t-th forecast time includes the fluid state vector of each grid unit 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, and the solution of the t+1-th forecast time is entered according to the underwater explosion bubble movement forecast result of the t-th forecast time and the grid division result of the t+1-th forecast time;

[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 intermittent galvanometer. Under the framework of underwater explosion bubble motion prediction based on intermittent 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 process. BRIEF DESCRIPTION OF THE DRAWINGS

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

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

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

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

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

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

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

[0022] Figure 8 The figure is a comparison diagram of the specific internal energy predicted by using the grid adaptive intermittent Ka Liu Jin method of the present application and the traditional fixed grid intermittent Ka Liu Jin method in an example.

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

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

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

[0026] The present application discloses a method for predicting underwater explosion bubble motion based on grid adaptive intermittent galaxy. In the process of using the intermittent galaxy method to solve the underwater explosion bubble motion prediction, the method adds a grid adaptive method to optimize the traditional intermittent galaxy 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 comprises 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, 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. In which, 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 conservation form and expressed as (three-dimensional):

[0032]

[0033] Among them, 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 then:

[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 equation (4) is combined with the EOS state equation shown in equation (5) to obtain the compressible fluid control equation.

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

[0043] Step 3, initializing parameter t=1, for any t-th forecast time, using 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 solving the problem using the discontinuous Galerkin method, a set of basis functions are selected in 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 grid refinement. Then, the constructed compressible fluid control equation is multiplied by a test function on both sides and integrated in the entire computational domain, thereby converting the partial differential equation into a weak form. This can reduce the smoothness requirements of the solution and enable discontinuous solutions to be properly handled. For the interface between adjacent grid cells, it is necessary to define a numerical flux to handle the solution at the interface. The numerical flux is constructed based on the physical conservation law and is used to approximate the true flux value to ensure the continuity and conservation of the solution on the cell boundary. For transient problems, it is necessary to select a suitable time integration scheme 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. Finally, the above steps will lead to a linear or nonlinear algebraic equation system 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 at the t-th forecast amount. This application does not elaborate on the specific solution method of the intermittent Galiokin method.

[0045] The underwater explosion bubble motion prediction result obtained at any t-th prediction time includes the fluid state vector of each grid unit in the grid division result at the t-th prediction time, 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 The coordinates inside, φ 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 unit is continuous, but discontinuity is allowed between two grid units, so that the conservation state vector Q in the entire computational domain is a sheet-like continuous function in the finite element space.

[0046] The coordinates in the basis function represent the grid cells The natural coordinates of the grid unit in the two-dimensional plane are x = (ξ, η), and the natural coordinates of the grid unit 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 is true for η 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 orthogonal 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 is true for other quadrature methods.

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

[0055]

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

[0057]

[0058] Combining the orthogonal conditions of the Jacobi polynomials, the orthogonal 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 functions are orthogonal on the isoparametric units. For simplicity, α = β = 0, so ω (α，β) (ξ)=1 and 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 the natural coordinate space is 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 the natural coordinate space 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, and the solution of the t+1-th forecast moment is entered according to the underwater explosion bubble motion forecast result of the t-th forecast moment and the grid division result of the t+1-th forecast moment.

[0076] In the traditional underwater explosion bubble motion forecasting process using intermittent galujin, after the grid division is completed in the initial state, the grid unit division result is always kept unchanged to iteratively solve the forecast, and the grid division result will affect the solution forecasting process and the forecast accuracy. However, the present application adds a grid adaptive adjustment moment in the underwater explosion bubble motion forecasting process. When the t-th forecast moment does not belong to the grid adaptive adjustment moment, the next forecast moment is directly entered for iterative solution as in the conventional practice. When the t-th forecast 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 forecasting 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 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 at the t+1-th forecast moment.

[0078] Determining the solution resolution of each grid cell in the grid division result at the tth 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 at 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 kth grid unit in each fission direction of the kth grid unit, 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. The density gradients ρ1, ρ2 and ρ3 of the fluid density ρ at each coordinate in the fluid state vector of the kth grid unit in the fission directions η1, η2, η3 are calculated respectively. Similarly, the density gradients of the fluid density ρ at other coordinates in the fission directions η1, η2, η3 can be obtained, thereby determining the density gradient ranges of different coordinates in the kth grid unit in each fission direction.

[0083] Then, the density gradient range of all grid cells in the entire computational domain in each fission direction is counted, and the relative value interval of the density gradient range of the kth grid cell in each fission direction in the density gradient range of the entire computational domain in the current fission direction is compared to obtain the solution resolution of the kth grid cell in the current fission direction. The larger the relative value interval of the density gradient range of the kth grid cell in each fission direction in the density gradient range of the entire computational domain in the current fission direction, the larger the corresponding solution resolution. The specific corresponding relationship can be customized. For example, when the density gradient range of the entire computational domain in the fission direction η1 is [1,100], when the density gradient range of the kth grid cell in the fission direction η1 is [40,50], the solution resolution in the fission direction η1 is determined to be 1, and when the density gradient range of the kth grid cell in the fission direction η1 is [80,90], the solution resolution in the fission direction η1 is determined to be 2. The specific corresponding relationship 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 value of the solution resolution of the kth grid unit in each fission direction as the overall resolution of the kth grid unit. For example, in the above example, the maximum value of the solution resolutions Θ1, Θ2, and Θ3 is taken as the overall resolution of the kth grid unit.

[0085] It should be noted that the fluid state vector includes other parameters in addition to the fluid density. The present 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. For example, the gradient of the fluid pressure is calculated to obtain the solution resolution, etc., which are all calculated in the same way. The present application does not limit this.

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

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

[0088] Determine that the grid cells that share nodes with the kth grid cell and are directly adjacent to it are the first-level adjacent cells of the kth 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 nodes with the kth grid cell in the three fission directions η1, η2, and η3. These 6 grid cells are the first-level adjacent cells of the kth grid cell. Continue to extrapolate based on this principle. For any integer parameter 1≤g≤λ k , determine that the other directly adjacent grid cells that share 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 g-th level neighboring unit has not yet calculated the solution resolution, k -g is used as the solution resolution for the g-th level neighboring unit. When the g-th level neighboring unit already has a solution resolution, k The larger value of the solution resolution of -g and the g-th level adjacent unit is used as the updated solution resolution of the g-th level adjacent unit.

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

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

[0091] (1) In case 1, if the solution resolution of the grid cell at the current grid adaptive adjustment time is higher than the solution resolution of the grid cell in the corresponding area at the last grid adaptive adjustment time, 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 time 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 a grid unit is split into multiple grid units, it is also necessary to adjust the underwater explosion bubble motion forecast result at the grid unit at the t-th forecast time accordingly to obtain the fluid state vector of each grid unit 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. Including: assigning the fluid state vector of the grid unit before the fission contained in the underwater explosion bubble motion forecast at the t-th forecast time to each grid unit after the fission.

[0094] For example, in an example, taking the grid unit as a cube, one case is to convert 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 unit shown in (a) is divided into 4 levels in each fission direction, thereby obtaining 64 grid units 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 time is lower than the solution resolution of the grid cell in the corresponding area at the last grid adaptive adjustment time, it means that the solution resolution of the area where the grid cell is located has been reduced. In this case, the grid cell in the grid division result at the t-th forecast time and several of its 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 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.

[0097] Similarly, when merging multiple grid units into one grid unit, it is also necessary to adjust the underwater explosion bubble motion forecast results of the multiple grid units at the t-th forecast time accordingly to obtain the fluid state vector of the merged grid unit, 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 units before merging contained in the underwater explosion bubble motion forecast at the t-th forecast time to the merged grid unit.

[0098] For example, in one example, 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, a 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. Figure 4 As shown in (f) in FIG. 1 , in another example, a grid unit and an adjacent grid unit along each fission direction are merged into one grid unit, so that Figure 4 The 64 grid cells shown in (c) are merged into 8 grid cells. 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 motion 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 forecasting process, in the underwater explosion bubble motion forecasting process, a grid adaptive adjustment moment is set every several forecast moments, 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. And after the grid units are fissioned at a grid adaptive adjustment moment, they may continue to fission or may be merged again, and the dynamic adjustment is made according to the actual solution resolution of the forecasting process. For example, in Figure 4 In the example shown, Figure 4 After the grid unit shown in (a) 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 calculation domain is [0,1], the initial number of grids is 200, and the fixed grid discontinuous Ga Liu Jin method and the grid adaptive discontinuous Ga Liu Jin method of this application are used for calculation. The comparison results of the distribution curves of flow field density, flow field velocity, flow field pressure and specific internal energy 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 Jia Liu Jin method of this application, and the black curve is the forecast result obtained by using the fixed grid discontinuous Jia Liu Jin method. Figure 5-Figure 8 It can be seen that in the forecast results obtained by using the grid-adaptive discontinuous Galin method, the capture accuracy of density discontinuities, velocity and pressure discontinuities is significantly improved, the numerical dispersion is reduced, and the amount of calculation is not greatly increased. In the process of forecasting using the grid-adaptive discontinuous Galin method of this application, 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 encryption effect at these two locations is as follows: Fig. 9 Shown and Fig.10 shown.

[0102] 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 motion based on grid adaptive discontinuous galaxy, characterized in that: The underwater explosion bubble motion prediction method comprises: Based on the fluid Euler equation combined with the EOS state equation, the compressible fluid control equation is constructed, the computational domain is discretized into multiple grid cells to obtain the grid division result of the first forecast time and determine the initial conditions and boundary conditions of the computational domain; Initialize the parameter t=1, for any t-th forecast time, use the discontinuous Galician 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; wherein the underwater explosion bubble motion forecast result at any t-th forecast time includes the fluid state vector of each grid unit 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, and the solution of the t+1-th forecast time is entered according to the underwater explosion bubble movement forecast result of the t-th forecast time and the grid division result of the t+1-th forecast time; 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, characterized in that: The grid division result at the t-th forecast time is adjusted based on the adjustment target of improving the solution resolution to obtain the grid division result at the t+1-th forecast time, including any grid unit in the grid division result at the t-th forecast time: 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 last grid adaptive adjustment moment, the grid unit in the grid division result of the t-th forecast moment is split into multiple grid units; wherein the first forecast moment is taken as the last 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 last 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, characterized in that: The adjustment of the underwater explosion bubble motion forecast result at the t-th forecast time based on the grid division result at the t+1-th forecast time includes: When a grid cell in the grid division result at the t-th forecast time is fissioned 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, characterized in that: When a grid unit in the grid division result of 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, characterized in that: 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, characterized in that: The solution resolution of each grid unit in the grid division result at the t-th forecast time is determined according to the underwater explosion bubble motion forecast result at the t-th forecast time, including: Each grid cell in the grid division result at the t-th forecast time is traversed in turn. 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 , and 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 other grid cells; until all grid cells in the grid division result at the t-th forecast time are traversed, the solution resolution of all grid cells is obtained.

7. The underwater explosion bubble motion prediction method according to claim 6, characterized in that: 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. k include: 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 are calculated respectively, and the solution resolution in the corresponding direction is calculated 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 unit in each fission direction as the overall resolution of the k-th grid unit; 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 between 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, characterized in that: 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 Updating the solution resolution of each other mesh cell includes: Determine the grid cells that share nodes 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 1≤g≤λ k , 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 λ k -g 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 7, characterized in that: According to λ k -g updates the solution resolution of the g-th level neighboring cells of the k-th grid cell, including: 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, set λ k -g as the solution resolution of the g-th level neighboring unit; when the g-th level neighboring unit has a solution resolution, k -g and the solution resolution of the g-th level adjacent unit, the larger value is used 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: In the process of underwater explosion bubble movement prediction, a grid self-adaptive adjustment moment is set every several prediction moments, and the grid division result is dynamically adjusted at each grid self-adaptive adjustment moment.

Images