CN109214082A - A kind of high resolution numerical simulation method of near field underwater blast wave load - Google Patents

A kind of high resolution numerical simulation method of near field underwater blast wave load Download PDF

Info

Publication number
CN109214082A
CN109214082A CN201811018574.6A CN201811018574A CN109214082A CN 109214082 A CN109214082 A CN 109214082A CN 201811018574 A CN201811018574 A CN 201811018574A CN 109214082 A CN109214082 A CN 109214082A
Authority
CN
China
Prior art keywords
grid
zoning
boundary
formula
explosive
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN201811018574.6A
Other languages
Chinese (zh)
Inventor
张之凡
王成
胡皓亮
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beijing Institute of Technology BIT
Original Assignee
Beijing Institute of Technology BIT
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Beijing Institute of Technology BIT filed Critical Beijing Institute of Technology BIT
Priority to CN201811018574.6A priority Critical patent/CN109214082A/en
Publication of CN109214082A publication Critical patent/CN109214082A/en
Pending legal-status Critical Current

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

This invention proposes a kind of high resolution numerical simulation method of near field underwater blast wave load, belongs to underwater explosion numerical simulation technology field.The present invention carries out high-precision numerical discretization to the non-linear Eulerian equation of meter and strong compressibility using the golden method (RKDG) of interruption gal the Liao Dynasty and three rank TVD Runge-Kutta methods;The physical quantity of gas-liquid interface two sides is pre-processed using ghost fluid method (GFM), effectively cuts down non-physical concussion caused by being interrupted by physical quantity;Using the complex topology structure variation occurred on Level Set Method (LSM) treated substance interface, the accurate capture to moving interface is realized.Method proposed by the present invention can effectively simulate the generation and communication process of near field underwater blast wave load;Numerical simulation result and classical empirical equation, experimental result can quite well, demonstrate this method and handling strong compressible, strong discontinuity, there is some superiority in transient state strong nonlinearity problem.

Description

A kind of high resolution numerical simulation method of near field underwater blast wave load
Technical field
This invention proposes a kind of high resolution numerical simulation method of near field underwater blast wave load, belongs to underwater Explosion numerical simulation technology field.
Background technique
It is increasingly prominent with marine safety problem, countries in the world competitively develop high speed, big dose, precise guidance elder generation Into Underwater Battery, great threat is constituted to large surface warships such as submarine and aircraft carriers.The shock wave that Underwater Battery explosion generates Load, pressure peak may be up to GPa magnitude, globality will be caused to injure destruction to naval ship structure.Therefore, effectively simulation is underwater Explosion wave load character is of great significance for naval vessels safety and Protection.In addition, Underwater Battery Performance Evaluation, ocean Also there is great demand in the fields such as military installations construction to this technology.To sum up, method proposed by the present invention has important skill Art background.
Traditional underwater explosion research mode mainly includes that theory analysis and experiment are tested.
Theory analysis is by establishing and solving the Closure equation including fluid motion equation, state of matter equation etc. Group obtains the analytic solutions of related physical quantity.This method in the research of early stage have certain guidance meaning, but simultaneously there is also More defect: (1) theoretical model is relatively simple, and research achievement can not solve actual engineering problem.(2) draw in derivation process The artificial it is assumed that there are certain errors between actual condition of some simplicity is entered.(3) it is related to the mathematic(al) manipulation of many complexity, Research difficulty is big, the period is long.
Experiment test is can the most directly to provide reliably for other correlative studys with accurate research method, result Data basis, the empirical equation of many classics are also derived from the summary to abundant experimental results.The deficiency of experiment test is: (1) duration of underwater explosion process, pressure peak reached as high as GPa magnitude in ms magnitude, and experiment measurement difficulty is very big. (2) explosion test has greatly destructive, and the loss of measurement of correlation instrument is very big, and experimental cost is higher.Meanwhile testing participant There is a certain security risk by member.(3) often there is uncertain factor in experiment, repeatability is poor.
With the continuous development of computer technology, numerical simulation has become the research essential method of underwater explosion problem One of.The advantage of the technology is: (1) solve the problems, such as the high-efficient of Practical Project, it is adaptable.(2) compared with experiment, numerical value That simulates is low in cost, and without security risk.(3) easy to operation.Currently, numerical simulation field of exploding under water, mainly There are two kinds of research modes: first is that using existing business software, such as LS-DYNA, AUTODYN, ABAQUS, though this mode Right simple possible, but its numerical result is often not accurate enough.Second is that the certain methods proposed using domestic and foreign scholars, such as finite element Grid model, mesh free SPH method, ALE algorithm etc., achieve many positive achievements.However, there are still two at present for these methods Item is insufficient: (1) the case where applicable far field is exploded mostly more lacks the research of near field underwater explosion;(2) many methods are adopted With low order precision format, the accuracy of calculated result is affected.Precision guided weapon high speed development, naval vessels security protection demand under water Under the overall background of surge, the above two o'clock deficiency is urgently made up.
For these deficiencies, the invention proposes the high resolution numerical simulation methods of near field underwater blast wave load. Specifically, this method is a kind of based on the golden method (RKDG) of interruption gal the Liao Dynasty, while Level Set Method (LSM) and virtually is coupled The unified algorithm of stream method (GFM).RKDG method can fluid motion equation carry out high-precision it is discrete;LSM method can have Imitate the complex topology structure variation at treated substance interface;GFM can effectively near inhibiting substances interface Nonphysical Oscillation.Pass through This method can carry out high-precision analog near field underwater blast wave load.
Summary of the invention
The invention proposes a kind of high resolution numerical simulation methods of near field underwater blast wave load.Its object is to The relative absence of current near field underwater blast wave load value simulation field is made up, and improves the precision of calculated result.
Above-mentioned purpose that the invention is realized by the following technical scheme:
Step 1: determining the zoning of underwater explosion, establish rectangular coordinate system, be m × n grid by the region division (wherein, m indicates that the direction x number of grid, n indicate the direction y number of grid).
Step 2: each physical quantity, correlation function in definition zoning under original state, and establish fluid governing equation And state of matter equation, comprising:
(a) interface function of each substance, i.e. Level Set function:
Wherein, Γ (0) indicates that zero contour surface, i.e. material interface, d (x, y, Γ (0)) indicate that a bit (x, y) is arrived in computational domain The distance of material interface, Ω1、Ω2The respectively region of material interface two sides.
(b) physical parameter of each substance: density, speed, the pressure of explosive, original state are lauched close under original state Degree, speed, pressure.
(c) according to the physical characteristic of underwater explosion field, viscous adaptable nothing, irrotationality, meter and can fluid governing equation: be established The fluid governing equation of compressibility:
Ut+ ▽ F (U)=S (U) (2)
Wherein,
(d) state equation (Tait equation) and relevant parameter of water:
Wherein, N, AW、BWIt is constant relevant to the physical property of water, ρW0Indicate the density of water in the initial state, pW、 ρWRespectively indicate the pressure and density of water.
(e) detonation products state equation (JWL equation) and relevant parameter:
Wherein, AE、BE、R1、R2, ω be constant relevant to explosive property, ρE0Indicate explosive in the initial state close Degree, ρEIndicate the density of detonation product, PEIndicate the density of detonation product, E0Indicate the specific internal energy of detonation product.
Step 3: expanding zoning, boundary condition is set according to different operating conditions: since boundary lacks difference scheme institute Required mesh point, therefore zoning is expanded to N layers outward.According to the format that this method uses, at least need to open up outward Open up 3 layers of grid (i.e. N >=3).
For free boundary condition, have: UΓ=U-1,U-i=U-i+1, i=2N, wherein UΓIndicate zoning side Random physical quantity on boundary's grid, U-iIndicate the random physical quantity on virtual grid.Specifically, first outside the boundary of zoning Physical quantity on layer virtual grid is equal to the physical quantity on the boundary mesh of zoning;Second layer virtual net outside the boundary of zoning Physical quantity on lattice is equal to the physical quantity on first layer virtual grid;Physics outside the boundary of zoning on third layer virtual grid Amount is equal to the physical quantity on second layer virtual grid, and so on.
For wall boundary condition, have:Wherein, UΓIndicate zoning Random physical quantity on boundary mesh in addition to speed, U-iIndicate the random physical quantity on virtual grid in addition to speed, VΓTable Show the speed on the boundary mesh of zoning, V-iIndicate the speed on virtual grid.Specifically, first outside the boundary of zoning Layer virtual grid on perpendicular in the speed and boundary mesh on boundary perpendicular to boundary velocity vector and be zero, other physical quantitys It is equal;Outside the boundary of zoning on second layer virtual grid perpendicular to boundary speed and first layer virtual grid on perpendicular to side The velocity vector on boundary and be zero, other physical quantitys are equal;Perpendicular to boundary on third layer virtual grid outside the boundary of zoning On speed and second layer virtual grid perpendicular to boundary velocity vector and be zero, other physical quantitys are equal, and so on.
Step 4: choose CFL parameter and calculate time step: CFL parameter is the constant between 0 to 1;It is current to calculate It is indicated with t=n Δ t total time used, n is to calculate step number, and time step Δ t can be obtained by following formula:
Δ x is the mesh spacing in the direction x in zoning;Δ y is the mesh spacing in the direction y in zoning;U is to calculate Velocity component of the current grid in the direction x in region;Velocity component of the v for current grid in zoning in the direction y;C is meter Calculate the velocity of sound in region at current grid.
Step 5: for the parameter assignment on zoning and virtual grid: as t=0, the calculating that will be provided in step 1 The physical parameter of each substance and interface function value in region, the parameter initial value as grid each in zoning;It, will as t > 0 The physical parameter of each substance and interface function value in the zoning being calculated in step 13, as net each in zoning Initial value of the lattice parameter when next step calculates.For virtual grid, then according to the method established in step 3, pass through interpolation and extrapolation Assignment is carried out to it.
Step 6: solving material interface normal direction: according to Level Set equation, to the normal of material interface each in zoning Direction is solved, solution formula are as follows:
Step 7: solving Riemannian problem and to grid assignment near material interface: establishing explosive along material interface normal direction Riemannian problem between water simultaneously solves the physical quantity between explosive and water at material interface.It is obtained with solution Riemannian problem Physical quantity carries out assignment close to one layer of grid of material interface in explosive region.Meanwhile by extrapolate to explosive region other than Grid carry out assignment, method particularly includes: at the flow field of treated substance interface explosive region side, by the material interface other side Grid be set as Level Set method grid, the one of the pressure on these Level Set method grids and normal velocity and the side real units It causes, and across the material interface interpolation and extrapolation of real fluid that tangential velocity and density then pass through explosive region acquires.
Step 8: by the discrete Eulerian equation of RKDG method, the numerical flux of positive pressure of shock wave is solved, specifically: to (2) Formula both sides are obtained with multiplied by tentative function Φ (x, y), and in grid cell K upper integral:
After carrying out integration by parts to (7) formula, the line integral and Line Integral in Gauss quadrature formula solution formula are utilized.Due to net Positive pressure of shock wave flux is discontinuous at lattice elementary boundary, therefore uses numerical fluxInstead of real fluxesUse approximate solutionInstead of accurately solving U, with the tentative function after numerical approximationInstead of tentative function Φ, It can thus be concluded that:
Wherein, ω indicates weight coefficient, | e | indicate the length of integral domain grid cell boundary, | K | indicate integral domain The area of grid cell.
For the numerical flux in (8) formulaIt is solved using (9) formula:
Wherein, α indicates Jacobi matrixThe maximum maximum value of characteristic value, U+(xel,yel, t) and it indicates The U value of grid cell boundary, U are approached from positive direction-(xel,yel, t) and indicate the U that grid cell boundary is approached from negative direction Value.For the pressure of explosive region detonation product, obtained by (4) formula in step 1.
For the approximate solution in (8) formulaIt is solved using (10) formula:
Wherein, for grid cell [xi-1/2,xi+1/2]×[yj-1/2,yj+1/2],Δxi=xi+1/2-xi-1/2, Δ yj=yj+1/2-yj-1/2
The semi-discrete scheme that (9), (10) two formulas are brought into available (2) formula of (8) formula, is denoted as:
Ut=R U (11)
Wherein, the Discrete Operator of R U representation space derivative term.
Step 9: it is discrete by time term progress of the three rank TVD Runge-Kutta methods to (11) formula equation left side, it obtains (2) the complete discrete scheme of formula, and then solve each physical quantity of detonation products in subsequent time zoning.Specifically Formula are as follows:
Step 10: similar operations are carried out with step 7, solve Riemannian problem again and to grid assignment near material interface: The Riemannian problem established between explosive and water along material interface normal direction simultaneously solves between explosive and water at material interface Physical quantity.The physical quantity obtained with solution Riemannian problem carries out assignment close to one layer of grid of material interface in air section. Meanwhile assignment is carried out to the grid other than air section by extrapolating, method particularly includes: in treated substance interfacial air region one When the flow field of side, Level Set method grid is set by the grid of the material interface other side, the pressure on these Level Set method grids It is consistent with the side real units with normal velocity, and tangential velocity and density then pass through the real fluid in explosive region across object The extrapolation of matter interface interpolation acquires.
Step 11: on the basis of step 10, carrying out the operation similar with step 8, step 9.But the difference is that: step In rapid 8, for the pressure of water area, obtained by (3) formula in step 1;What solution obtained in step 9 is that subsequent time calculates Each physical quantity of water in region.
Step 12: the space derivation item of the Level Set equation of motion is carried out discrete:
For the lateral numerical flux of the interface function of explosive and water, can be solved by (14) formula:
Wherein,Indicate grid cell in zoningIn the interface function transverse direction of upper explosive and water Numerical flux,Indicate the interface function numerical flux in lateral positive direction on the grid cell,Indicate the net Interface function numerical flux in lateral negative direction on lattice unit;Indicate grid cell (i-r+s, j) in zoning The flux of the lateral positive direction of the interface function of upper explosive and water;It indicates in zoning on grid cell (i-r+s, j) The flux of the lateral negative direction of the interface function of explosive and water.Design parameter can be solved by (15), (16) formula:
Wherein,Indicate the interface function value of explosive and water on grid cell (i, j) in zoning.
For the Vertical Numerical flux of the interface function of explosive and water, method for solving is identical.
By step 12, the semi-discrete scheme of available (13) formula is denoted as:
Wherein,The Discrete Operator of representation space derivative term:
Step 13: by the time-derivative item of the three discrete Level Set equations of motion of rank TVD Runge-Kutta method, The complete discrete scheme of (13) formula is obtained, and then solves the interface function value of each grid cell in subsequent time zoning.
Step 14: on the basis of step before, solve and export the pressure value of all grid cells in zoning: if Grid cell is in detonation products region, then is acquired by (4) the formula calculating in step 1;If grid cell is in water In region, then acquired by (3) formula in step 1.
Step 15: the current calculating total time t for calculating total time t=n Δ t and setting used of judgement*Relationship: if t > t*, then terminate to calculate, export the pressure value of grid cell in zoning;If t > t*, then it is back to step 3, continues operation.
Beneficial effect
Underwater blast wave load high resolution numerical simulation method near field proposed by the present invention and traditional underwater explosion are rushed It hits wave load method for numerical simulation to compare, have the advantage that
(1) present invention carries out Eulerian equation space derivation item using RKDG method discrete, using three rank TVD Runge- Kutta method is discrete to the progress of Eulerian equation time-derivative item, has higher precision compared to conventional method;
(2) physical quantity discontinuous problem of the present invention using the method treated substance interface the DGFM two sides in GFM method group, height Maintain to resolution ratio the discontinuous nature of material interface;DGFM method can effectively inhibit to solve compared to original GFM method The non-physical reforming phenomena that interface generates when Eulerian equation, ensures the stable operation of algorithm;
(3) the Level Set method that the present invention uses can accurately capture moving interface, with traditional Lagrange circle Face method for tracing is compared, and Level Set method more stable can effectively handle complex material interface and its topological structure change Change, when the Nonlinear Large Deformations such as broken, tearing, fusion for especially occurring during solving the problems, such as underwater explosion, has obviously Advantage.
(4) method proposed by the present invention can effectively simulate near field underwater blast wave load, accurately reflect underwater quick-fried For fried shock wave in the propagation law near field, numerical result and experimental result meet ground preferably;With the tradition mainly for middle far field Method for numerical simulation is compared, and the scope of application is bigger, more adaptable.
Detailed description of the invention
Fig. 1 is Riemannian problem schematic diagram
Fig. 2 is to construct local Riemannian problem schematic diagram
Fig. 3 is Level Set method schematic diagram
Fig. 4 is DGFM method schematic diagram
Fig. 5 is the initial geometric model figure of example
Fig. 6 is the pressure time-process curve graph in example at measuring point B
Fig. 7 is underwater blast wave simulation effect picture
Specific embodiment
The present invention is specifically described below with reference to example, the initial geometric model in example is presented in figure 5, is seen Pressure time-process curve at measuring point B is presented in figure 6, and underwater blast wave simulation effect is presented in fig. 7.
Information needed for steps 1 and 2 is provided first: as shown in Fig. 5, using O point as coordinate origin, zoning size A=1000mm, b=1800mm, transverse grid number are 500, and longitudinal grid number is 900, mesh spacing 2mm.The total time of setting For 0.6ms, CFL number takes 0.2.Powder charge model is spherical charge, and centre coordinate is (0,0.9), radius r=11.15mm.Take measuring point The coordinate of B is (0.5,0.9).Explosive area equation uses JWL state equation, and water area uses Tait state equation, specific object The selection of reason parameter see the table below.
The JWL state equation parameter of 1 TNT explosive of table
AE/Pa BE/Pa R1 R2 ω ρE0/(kg/m3) E0/(J/m3)
3.712×1011 3.231×1109 4.15 0.95 0.3 1630 4.29×106
2 Tait state equation parameter of table
N Aw/Pa Bw/Pa ρW0/(kg/m3)
7.15 1.0×105 3.31×108 1000
Steps 1 and 2 are completed as a result,.
Then, step 3 is carried out to step 15, obtains the numerical simulation result of underwater blast wave load.
Calculated result analysis:
Attached drawing 6 is the pressure time-process curve at observation point B, and solid line is the experiment value that pressure sensor measures, dotted line in figure The pressure value calculated for numerical value.It can be seen that shock motion is to observation point B, and pressure increases suddenly at this time when t ≈ 0.30ms Greatly to about 19.8MPa, the attenuation law of exponential form is presented in shock wave in water later.By correlation curve it is found that by we The pressure trend and peak value that method is calculated match with experimental result, error is less than 5%, it was demonstrated that the numerical simulation side Method is feasible effective.
Attached drawing 7 is that the underwater blast wave that this method is calculated simulates effect picture, is clearly labelled with impact in figure Wave surface, air detonation product interface and cohesion rarefaction wave region.The biography as can be seen that wave surface and material interface explode under water Good geometrical symmetry is maintained during broadcasting, interface svelteness meets objective physical rule, illustrates that this method can be quasi- Really effectively capture underwater explosion material interface motion conditions.
The drawings and the specific embodiments two parts are combined above, and the present invention is elaborated.But it needs to state Be: these pictures and text are merely illustrative, and any restrictions are not constituted to protection scope of the present invention, without departing from this hair In the case where bright spirit and scope, those skilled in the art carry out a variety of improvement, equivalencing to the present invention and embodiments thereof Or modification, these fall within the protection scope of the present invention.Protection scope of the present invention is subject to the appended claims.

Claims (4)

1. a kind of high resolution numerical simulation method of near field underwater explosion loading, it is characterised in that: its operating process are as follows:
Step 1: determining the zoning of underwater explosion, establish rectangular coordinate system, be m × n grid (its by the region division In, m indicates that the direction x number of grid, n indicate the direction y number of grid);
Step 2: each physical quantity, correlation function in definition zoning under original state, and establish fluid governing equation and object Matter state equation, comprising:
(a) interface function of each substance, i.e. Level Set function:
Wherein, Γ (0) indicates that zero contour surface, i.e. material interface, d (x, y, Γ (0)) indicate that a bit (x, y) arrives substance in computational domain The distance at interface, Ω1、Ω2The respectively region of material interface two sides;
(b) physical parameter of each substance: density, speed, the pressure of explosive under original state, density that original state is lauched, speed Degree, pressure;
(c) according to the physical characteristic of underwater explosion field, viscous adaptable nothing, irrotationality, meter and compressible fluid governing equation: are established The fluid governing equation of property:
Ut+ ▽ F (U)=S (U) (2)
Wherein,F (U)=[f (U), g (U)],
(d) state equation (Tait equation) and relevant parameter of water:
Wherein, N, AW、BWIt is constant relevant to the physical property of water, ρW0Indicate the density of water in the initial state, pW、ρWPoint Not Biao Shi water pressure and density;
(e) detonation products state equation (JWL equation) and relevant parameter:
Wherein, AE、BE、R1、R2, ω be constant relevant to explosive property, ρE0Indicate the density of explosive in the initial state, ρETable Show the density of detonation product, PEIndicate the density of detonation product, E0Indicate the specific internal energy of detonation product;
Step 3: expand zoning, according to different operating conditions be arranged boundary condition: by boundary lack difference scheme institute it is required Mesh point, therefore zoning is expanded to N layers outward;According to the format that this method uses, at least need to expand 3 layers outward Grid (i.e. N >=3);
For free boundary condition, have: UΓ=U-1,U-i=U-i+1, i=2N, wherein UΓIndicate zoning boundary net Random physical quantity on lattice, U-iIndicate the random physical quantity on virtual grid;Specifically, first layer is empty outside the boundary of zoning Physical quantity on quasi- grid is equal to the physical quantity on the boundary mesh of zoning;Outside the boundary of zoning on second layer virtual grid Physical quantity be equal to first layer virtual grid on physical quantity;Physical quantity etc. outside the boundary of zoning on third layer virtual grid Physical quantity on second layer virtual grid, and so on;
For wall boundary condition, have:Wherein, UΓIndicate zoning boundary net Random physical quantity on lattice in addition to speed, U-iIndicate the random physical quantity on virtual grid in addition to speed, VΓIt indicates to calculate Speed on the grid of zone boundary, V-iIndicate the speed on virtual grid;Specifically, first layer is virtual outside the boundary of zoning On grid perpendicular in the speed and boundary mesh on boundary perpendicular to boundary velocity vector and be zero, other physical quantitys are equal; Outside the boundary of zoning on second layer virtual grid perpendicular to boundary speed and first layer virtual grid on perpendicular to boundary Velocity vector and be zero, other physical quantitys are equal;Speed outside the boundary of zoning on third layer virtual grid perpendicular to boundary With on second layer virtual grid perpendicular to boundary velocity vector and be zero, other physical quantitys are equal, and so on;
Step 4: choose CFL parameter and calculate time step: CFL parameter is the constant between 0 to 1;Used in current calculating Indicated with t=n Δ t total time, n is to calculate step number, and time step Δ t can obtain by following formula:
Δ x is the mesh spacing in the direction x in zoning;Δ y is the mesh spacing in the direction y in zoning;U is zoning Velocity component of the middle current grid in the direction x;Velocity component of the v for current grid in zoning in the direction y;C is to calculate area Velocity of sound in domain at current grid;
Step 5: for the parameter assignment on zoning and virtual grid: as t=0, the zoning that will be provided in step 1 The physical parameter and interface function value of interior each substance, the parameter initial value as grid each in zoning;As t > 0, by step The physical parameter of each substance and interface function value in the zoning being calculated in 13, as grid parameter each in zoning Initial value when next step calculates;For virtual grid, then according in step 3 establish method, by interpolation and extrapolation to its into Row assignment;
Step 6: solving material interface normal direction: according to Level Set equation, to the normal direction of material interface each in zoning It is solved, solution formula are as follows:
Step 7: solving Riemannian problem and to grid assignment near material interface: establishing explosive and water along material interface normal direction Between Riemannian problem and solve the physical quantity between explosive and water at material interface;The physics obtained with solution Riemannian problem Amount carries out assignment close to one layer of grid of material interface in explosive region;Meanwhile by extrapolating to the net other than explosive region Lattice carry out assignment, method particularly includes: at the flow field of treated substance interface explosive region side, by the net of the material interface other side Lattice are set as Level Set method grid, and the pressure on these Level Set method grids is consistent with the side real units with normal velocity, And across the material interface interpolation and extrapolation of real fluid that tangential velocity and density then pass through explosive region acquires;
Step 8: by the discrete Eulerian equation of RKDG method, the numerical flux of positive pressure of shock wave is solved, specifically: to (2) formula two While being obtained with multiplied by tentative function Φ (x, y), and in grid cell K upper integral:
After carrying out integration by parts to (7) formula, the line integral and Line Integral in Gauss quadrature formula solution formula are utilized;Due to grid list First boundary positive pressure of shock wave flux is discontinuous, therefore uses numerical fluxInstead of real fluxesUse approximate solutionInstead of accurately solving U, with the tentative function after numerical approximationInstead of tentative function Φ, It can thus be concluded that:
Wherein, ω indicates weight coefficient, | e | indicate the length of integral domain grid cell boundary, | K | indicate integral domain grid The area of unit;
For the numerical flux in (8) formulaIt is solved using (9) formula:
Wherein, α indicates Jacobi matrixThe maximum maximum value of characteristic value, U+(xel,yel, t) and it indicates from just The U value of grid cell boundary, U are approached on direction-(xel,yel, t) and indicate the U value that grid cell boundary is approached from negative direction;It is right In the pressure of explosive region detonation product, obtained by (4) formula in step 1;
For the approximate solution in (8) formulaIt is solved using (10) formula:
Wherein, for grid cell [xi-1/2,xi+1/2]×[yj-1/2,yj+1/2],Δ xi=xi+1/2-xi-1/2, Δ yj=yj+1/2-yj-1/2
The semi-discrete scheme that (9), (10) two formulas are brought into available (2) formula of (8) formula, is denoted as:
Ut=R U (11)
Wherein, the Discrete Operator of R U representation space derivative term;
Step 9: it is discrete by time term progress of the three rank TVD Runge-Kutta methods to (11) formula equation left side, obtain (2) formula Complete discrete scheme, and then solve each physical quantity of detonation products in subsequent time zoning;Specific formula are as follows:
Step 10: carrying out similar operations with step 7, solve Riemannian problem again and to grid assignment near material interface: along object Riemannian problem that matter interface normal direction is established between explosive and water simultaneously solves the physics between explosive and water at material interface Amount;The physical quantity obtained with solution Riemannian problem carries out assignment close to one layer of grid of material interface in air section;Together When, assignment is carried out to the grid other than air section by extrapolating, method particularly includes: in treated substance interfacial air region side Flow field when, set Level Set method grid for the grid of the material interface other side, the pressure on these Level Set method grids with Normal velocity is consistent with the side real units, and tangential velocity and density then pass through the real fluid in explosive region across substance Interface interpolation extrapolation acquires;
Step 11: on the basis of step 10, carrying out the operation similar with step 8, step 9;But the difference is that: step 8 In, for the pressure of water area, obtained by (3) formula in step 1;What solution obtained in step 9 is subsequent time zoning Each physical quantity of interior water;
Step 12: it is discrete to the space derivation item progress of the Level Set equation of motion, it can obtain:
For the lateral numerical flux of the interface function of explosive and water, can be solved by (14) formula:
Wherein,Indicate grid cell in zoningNumerical value in the interface function transverse direction of upper explosive and water Flux,Indicate the interface function numerical flux in lateral positive direction on the grid cell,Indicate the grid list Interface function numerical flux in lateral negative direction in member;It indicates in zoning on grid cell (i-r+s, j) The flux of the lateral positive direction of the interface function of explosive and water;It indicates in zoning on grid cell (i-r+s, j) The flux of the lateral negative direction of the interface function of explosive and water;Design parameter can be solved by (15), (16) formula:
Wherein,Indicate the interface function value of explosive and water on grid cell (i, j) in zoning;
For the Vertical Numerical flux of the interface function of explosive and water, method for solving is identical;
By step 12, the semi-discrete scheme of available (13) formula is denoted as:
Wherein,The Discrete Operator of representation space derivative term:
Step 13: by the time-derivative item of the three discrete Level Set equations of motion of rank TVD Runge-Kutta method, obtaining (13) the complete discrete scheme of formula, and then solve the interface function value of each grid cell in subsequent time zoning;
Step 14: on the basis of step before, solving and export the pressure value of all grid cells in zoning: if grid Unit is in detonation products region, then is acquired by (4) the formula calculating in step 1;If grid cell is in water area It is interior, then it is acquired by (3) formula in step 1;
Step 15: the current calculating total time t for calculating total time t=n Δ t and setting used of judgement*Relationship: if t > t*, Then terminate to calculate, exports the pressure value of grid cell in zoning;If t > t*, then it is back to step 3, continues operation.
By the operation of above-mentioned steps, the shock wave in underwater explosion zoning at each moment each grid cell can be obtained Pressure value.
2. a kind of high resolution numerical simulation method of near field underwater blast wave load as described in claim 1, feature It is: extrapolation assignment is carried out to grid cell near material interface using Riemann problem solution in its step 7, method particularly includes:
At the flow field of treated substance interface explosive region side, Level Set method net is set by the grid of the material interface other side Lattice, the pressure and normal velocity on these Level Set method grids are consistent with the side real units, and tangential velocity and density Then acquired by across the material interface interpolation and extrapolation of real fluid in explosive region.
3. a kind of high resolution numerical simulation method of near field underwater blast wave load as described in one of claims 1 or 2, It is characterized by: passing through the joint of step 7, step 11 and step 12, have the two or more medium interaction problems of processing Function.
4. a kind of high resolution numerical simulation method of near field underwater blast wave load as described in claims 1 to 3, special Sign is: the discrete Eulerian equation of RKDG method is used in its step 8, it is only necessary to read grid cell self information and participate in calculating i.e. The above precision of three ranks can be obtained.
CN201811018574.6A 2018-09-03 2018-09-03 A kind of high resolution numerical simulation method of near field underwater blast wave load Pending CN109214082A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811018574.6A CN109214082A (en) 2018-09-03 2018-09-03 A kind of high resolution numerical simulation method of near field underwater blast wave load

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811018574.6A CN109214082A (en) 2018-09-03 2018-09-03 A kind of high resolution numerical simulation method of near field underwater blast wave load

Publications (1)

Publication Number Publication Date
CN109214082A true CN109214082A (en) 2019-01-15

Family

ID=64987167

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811018574.6A Pending CN109214082A (en) 2018-09-03 2018-09-03 A kind of high resolution numerical simulation method of near field underwater blast wave load

Country Status (1)

Country Link
CN (1) CN109214082A (en)

Cited By (15)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109902376A (en) * 2019-02-25 2019-06-18 北京理工大学 A kind of fluid structurecoupling high resolution numerical simulation method based on Continuum Mechanics
CN110879919A (en) * 2019-11-18 2020-03-13 中国人民解放军陆军防化学院 Sectional type simulation method for poison diffusion under explosion action
CN111008492A (en) * 2019-11-22 2020-04-14 电子科技大学 High-order unit Euler equation numerical simulation method based on Jacobian-free matrix
CN111783365A (en) * 2020-06-04 2020-10-16 三多(杭州)科技有限公司 Virtual medium method, device and equipment applied to low-voltage interface processing
CN112395722A (en) * 2019-08-15 2021-02-23 武汉理工大学 Method for acquiring motion response of hull beam under action of underwater explosion and wave load
CN112861263A (en) * 2021-02-22 2021-05-28 西北工业大学 Calculation simulation method suitable for compressible two-phase flow
CN112861445A (en) * 2020-12-29 2021-05-28 中国船舶重工集团公司第七研究院 Simulation method of grid-free numerical model
CN113408168A (en) * 2021-06-18 2021-09-17 北京理工大学 High-precision numerical simulation method based on Riemann problem accurate solution
CN113515900A (en) * 2021-04-07 2021-10-19 南京航空航天大学 Numerical simulation method for simulating wake flow of continuous rotation detonation engine
CN113553739A (en) * 2021-07-06 2021-10-26 西安近代化学研究所 Method for calculating explosion output characteristics of mixed explosive
CN113591345A (en) * 2021-07-08 2021-11-02 北京理工大学 Explosion reaction flow high-precision prediction method based on generalized Riemann solution method device
CN114595603A (en) * 2022-03-01 2022-06-07 哈尔滨工程大学 Underwater structure impact resistance calculation method considering underwater explosion bubble jet flow slamming
CN114647990A (en) * 2022-02-22 2022-06-21 哈尔滨工程大学 Method for efficiently and highly accurately calculating pressure wavelet of air gun seismic source
CN116306364A (en) * 2023-03-13 2023-06-23 武汉理工大学 Wave-absorbing simulation method for explosion water mist in cabin
CN116992796A (en) * 2023-09-27 2023-11-03 中国科学技术大学 Self-adaptive low-dissipation SPH-HLLC Riemann solver coupling algorithm

Cited By (24)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109902376A (en) * 2019-02-25 2019-06-18 北京理工大学 A kind of fluid structurecoupling high resolution numerical simulation method based on Continuum Mechanics
CN112395722B (en) * 2019-08-15 2022-05-31 武汉理工大学 Method for acquiring motion response of hull beam under action of underwater explosion and wave load
CN112395722A (en) * 2019-08-15 2021-02-23 武汉理工大学 Method for acquiring motion response of hull beam under action of underwater explosion and wave load
CN110879919A (en) * 2019-11-18 2020-03-13 中国人民解放军陆军防化学院 Sectional type simulation method for poison diffusion under explosion action
CN110879919B (en) * 2019-11-18 2023-08-18 中国人民解放军陆军防化学院 Sectional simulation method for poison diffusion under explosion effect
CN111008492A (en) * 2019-11-22 2020-04-14 电子科技大学 High-order unit Euler equation numerical simulation method based on Jacobian-free matrix
CN111008492B (en) * 2019-11-22 2022-10-14 电子科技大学 High-order unit Euler equation numerical simulation method based on Jacobian-free matrix
CN111783365A (en) * 2020-06-04 2020-10-16 三多(杭州)科技有限公司 Virtual medium method, device and equipment applied to low-voltage interface processing
CN111783365B (en) * 2020-06-04 2023-09-22 三多(杭州)科技有限公司 Virtual medium method, device and equipment applied to low-voltage interface processing
CN112861445A (en) * 2020-12-29 2021-05-28 中国船舶重工集团公司第七研究院 Simulation method of grid-free numerical model
CN112861263A (en) * 2021-02-22 2021-05-28 西北工业大学 Calculation simulation method suitable for compressible two-phase flow
CN112861263B (en) * 2021-02-22 2024-02-13 西北工业大学 Calculation simulation method suitable for compressible two-phase flow
CN113515900A (en) * 2021-04-07 2021-10-19 南京航空航天大学 Numerical simulation method for simulating wake flow of continuous rotation detonation engine
CN113408168A (en) * 2021-06-18 2021-09-17 北京理工大学 High-precision numerical simulation method based on Riemann problem accurate solution
CN113553739A (en) * 2021-07-06 2021-10-26 西安近代化学研究所 Method for calculating explosion output characteristics of mixed explosive
CN113591345B (en) * 2021-07-08 2024-01-23 北京理工大学 Explosion reaction flow high-precision prediction method based on generalized Riemann solver
CN113591345A (en) * 2021-07-08 2021-11-02 北京理工大学 Explosion reaction flow high-precision prediction method based on generalized Riemann solution method device
CN114647990A (en) * 2022-02-22 2022-06-21 哈尔滨工程大学 Method for efficiently and highly accurately calculating pressure wavelet of air gun seismic source
CN114595603A (en) * 2022-03-01 2022-06-07 哈尔滨工程大学 Underwater structure impact resistance calculation method considering underwater explosion bubble jet flow slamming
CN114595603B (en) * 2022-03-01 2023-06-13 哈尔滨工程大学 Underwater structure impact resistance calculation method considering underwater explosion bubble jet slamming
CN116306364A (en) * 2023-03-13 2023-06-23 武汉理工大学 Wave-absorbing simulation method for explosion water mist in cabin
CN116306364B (en) * 2023-03-13 2024-03-22 武汉理工大学 Wave-absorbing simulation method for explosion water mist in cabin
CN116992796B (en) * 2023-09-27 2023-12-22 中国科学技术大学 Self-adaptive low-dissipation SPH-HLLC Riemann solver coupling method
CN116992796A (en) * 2023-09-27 2023-11-03 中国科学技术大学 Self-adaptive low-dissipation SPH-HLLC Riemann solver coupling algorithm

Similar Documents

Publication Publication Date Title
CN109214082A (en) A kind of high resolution numerical simulation method of near field underwater blast wave load
Yan et al. Experimental and numerical investigation of water impact on air-launched AUVs
Cheng et al. Towards the modeling of the ditching of a ground-effect wing ship within the framework of the SPH method
Zhang et al. SPH-BEM simulation of underwater explosion and bubble dynamics near rigid wall
CN106124349A (en) Determine the critical erosion rate of gas well, the method and device of flow
Nguyen et al. Multiphase flow simulation of water-entry and-exit of axisymmetric bodies
Liang et al. Numerical study of spherical blast-wave propagation and reflection
Faucher et al. High resolution adaptive framework for fast transient fluid-structure interaction with interfaces and structural failure–Application to failing tanks under impact
Rebelo et al. A Comparison between three air blast simulation techniques in LS-DYNA
Ning et al. Numerical analysis of the shaped charged jet with large cone angle
Shenshen et al. Two-phase flow characteristics of different charging structure in ignition and flame-spreading
Ren et al. Numerical studies of penetration problems by an improved particle method
Mu et al. Numerical simulation on the cavitation flow of high speed oblique water entry of revolution body
Yang et al. Characteristics and load reduction method of the deep-sea implosion of the ellipsoidal and egg-shaped pressure hulls
Georgievskiy et al. Shock focusing upon interaction of a shock with a cylindrical dust cloud
Grządziela Ship impact modeling of underwater explosion
Hsiao et al. Development of compressible-incompressible link to efficiently model bubble dynamics near floating body
Zhu et al. Modeling and analysis of numerical wave tank based on the ALE algorithm
Pappalardo et al. Comparison of CFD numerical approaches for the simulation of accidental gas release in energy applications
Chen et al. Transient Response of the Steel Plate to Underwater Explosion Bubble
Lu et al. Risk assessment method and protection goals of high concrete gravity dam subjected to far-field underwater nuclear explosion
Jie et al. Numerical simulation on underwater explosion in small-sized containers
Charles et al. A numerical study on cavity expansion in water: hydraulic ram under ballistic impacts
Wei Development of a new steam explosion model for the MC3D software
Duan et al. Design and Implementation of a Damage Assessment System for Large-Scale Surface Warships Based on Deep Learning

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
WD01 Invention patent application deemed withdrawn after publication

Application publication date: 20190115

WD01 Invention patent application deemed withdrawn after publication