CN113591345A - Explosion reaction flow high-precision prediction method based on generalized Riemann solution method device - Google Patents

Explosion reaction flow high-precision prediction method based on generalized Riemann solution method device Download PDF

Info

Publication number
CN113591345A
CN113591345A CN202110771451.5A CN202110771451A CN113591345A CN 113591345 A CN113591345 A CN 113591345A CN 202110771451 A CN202110771451 A CN 202110771451A CN 113591345 A CN113591345 A CN 113591345A
Authority
CN
China
Prior art keywords
interface
detonation
wave
deflagration
conservation
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.)
Granted
Application number
CN202110771451.5A
Other languages
Chinese (zh)
Other versions
CN113591345B (en
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 CN202110771451.5A priority Critical patent/CN113591345B/en
Publication of CN113591345A publication Critical patent/CN113591345A/en
Application granted granted Critical
Publication of CN113591345B publication Critical patent/CN113591345B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • General Engineering & Computer Science (AREA)
  • Geometry (AREA)
  • Evolutionary Computation (AREA)
  • Fluid Mechanics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computing Systems (AREA)
  • Algebra (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

The invention discloses a high-precision prediction method for an explosion reaction flow based on a generalized Riemannian solution device, and belongs to the technical field of explosion mechanics. The method is realized by combining a two-step and four-step Lax-Wendroff method with a conservation interface method. The two-step and four-order Lax-Wendroff method integrates the simplicity characteristic of Runge-Kutta algorithm and the space-time coupling property of the generalized Riemann solver, can realize high-precision construction of numerical flux, and can reduce the time advancing step of the four-order precision algorithm by half. The conservation interface method is developed based on a combined algorithm of a level set method and a virtual fluid method. The level set method can quickly and accurately process the topological change of the interface and accurately capture the evolution process of the multi-medium interface; the virtual fluid method can skillfully convert the multi-medium problem into the single-medium problem, and effectively inhibit the non-physical oscillation of the multi-medium interface problem; the conservation interface method can improve the problem of non-conservation of the multi-medium interface and improve the reliability of a prediction result.

Description

Explosion reaction flow high-precision prediction method based on generalized Riemann solution method device
Technical Field
The invention relates to a high-precision prediction method for an explosion reaction flow based on a generalized Riemannian solution device, and belongs to the technical field of explosion mechanics.
Background
Explosion is accompanied by a violent chemical reaction, involves a fluid dynamic phenomenon under extreme conditions such as high speed, high temperature and high pressure, and involves a process of coupling various media such as gas, liquid and solid with each other. The research on the explosion phenomenon plays an important role in designing and damaging the weapons and ammunitions in the national defense industrial field, preventing and reducing the explosion accidents in the industrial field and the like. The experiment is an effective method for researching the explosion phenomenon, but has the defects of high danger, high cost, high equipment requirement and the like. Therefore, with the rapid development of computer science, especially after a large-scale parallel processing system is widely applied, a numerical prediction method with the advantages of safety, economy, easy operation and the like has been developed as an important means for researching explosion phenomena.
For fluidic problems involving the coupling of multiple media to each other, the location of the multi-media interface can be determined by a moving interface tracking technique. The Front tracking method can accurately capture the position of the interface, but the calculation process is very complicated. The VOF method can ensure the conservation of the consistency in the process of processing the topological change of the interface, but the high-order precision is difficult to realize. The level set method can accurately describe the topological change process of the interface and can quickly realize high-order precision, so that the method is widely applied to the fields of fluid calculation, image processing, computer vision and the like.
Level set methods are often combined with virtual fluid methods to convert multi-media problems to single-media problems through the construction of virtual fluids. The level set/virtual fluid method can effectively suppress non-physical oscillation of the predicted result, but has the obvious defect of lack of conservation. The conservation interface method improves the traditional level set/virtual fluid method, the level set method is applied to track the multi-medium interface, and the non-conservation problem can be effectively improved and the reliability of the prediction result can be improved through a modified finite volume method discrete control equation.
Solving numerical flux is a key step of a finite volume method, and a Generalized Riemannian Problem (GRP) method is used as a second-order expansion format of the Godunov method and has strong space-time coupling. The basic idea is to use a piece-wise linear function to represent the solution of the boundary of a computing unit and to construct numerical flux by solving a generalized Riemannian problem. The two-step and four-order Lax-Wendroff method based on the generalized Riemann solution device integrates the simplicity characteristic of Runge-Kutta algorithm and the space-time coupling property of the generalized Riemann solution device, four-order precision can be realized by only advancing two steps in the time direction, the calculation efficiency and compactness can be effectively improved, and the method is mainly applied to the single-medium fluid motion problem.
Disclosure of Invention
Aiming at the problems of low calculation efficiency, missing of conservation and the like of an explosion reaction flow prediction method in the prior art, the invention discloses a high-precision prediction method of an explosion reaction flow based on a generalized Riemannian solver, which aims to solve the technical problems that: the high-precision numerical flux construction method based on the generalized Riemann solution is realized, and the calculation efficiency of the high-precision prediction method is improved; meanwhile, the problem of non-conservation is effectively solved, the accurate prediction of detonation and deflagration processes is realized, and the technical problems of related engineering in the field of explosion mechanics are further solved. The related engineering technical problems in the field of explosion mechanics comprise weapon ammunition design and damage assessment and disaster prevention and reduction of explosion accidents.
The purpose of the invention is realized by the following technical scheme:
the invention discloses a high-precision prediction method of explosion reaction flow based on a generalized Riemannian solver, which is realized by combining a two-step four-order Lax-Wendroff method and a conservation interface method. And solving the topological changes of the detonation interface and the detonation interface by using a level set method, and quickly and accurately tracking the interface evolution process. And constructing virtual fluid in the computing units near the detonation interface and the detonation interface by using a virtual fluid method, and converting the multi-medium problem into the single-medium problem. An explosion reaction flow conservation type discrete format is constructed by utilizing a conservation interface method, and a conservation quantity exchange term along an interface is constructed by solving a Riemann problem containing detonation waves and deflagration waves. And combining a two-step and four-step Lax-Wendroff method and a conservation interface method to construct a high-precision numerical flux solving method based on a generalized Riemann solver. The method can improve the calculation efficiency of the high-precision prediction method of the explosion reaction flow; meanwhile, the problem of non-conservation of a multi-medium interface is effectively solved, and accurate prediction of detonation and deflagration processes of compressible reaction fluid is achieved.
The invention discloses a high-precision prediction method for explosion reaction flow based on a generalized Riemannian solution device, which comprises the following steps:
step 1: a rectangular coordinate system is established, space steps in the x direction and the y direction are defined, and (M +1) × (N +1) calculation units are arranged in a calculation region and are marked as units (i, j), wherein i is 0,1, …, and M, j is 0,1, …, N.
Step 2: the level set functions of all the calculation units in the calculation area at the initial time and the fluid physical quantities of unburned gas and burned gas, including density, moving speed in the x direction, moving speed in the y direction, and pressure, are defined according to the initial conditions.
In the level set method, the level set function is a symbol distance function, denoted by φ. Let regions where + 0 and + 0 denote distribution regions of two fluids, i.e., unburned gas and burned gas, respectively, and a position where + 0 denotes a detonation or deflagration interface. At the initial time, the value of the level set function φ is defined using the following equation:
Figure BDA0003153637430000021
wherein gamma is0A distribution curve representing the interface at the initial moment, d (x, y) representing the center (x, y) of the computing unit to the initial interface Γ0Distance of (d), omega1、Ω2Respectively, the areas of distribution of unburned gas and burned gas.
According to the initial conditions, the density, the moving speed in the x direction, the moving speed in the y direction and the pressure of unburned gas and burned gas on both sides of the initial interface are defined.
And step 3: a system of governing equations for the compressible reactive fluid is established and boundary conditions are set.
Step 3.1: a system of governing equations for compressible reactive fluids is established.
The CJ model is applied to establish a system of control equations for compressible fluids, i.e., without considering chemical reactivity, detonation and deflagration interfaces are considered as jump discontinuities, chemical reactants are immediately converted to products along the interfaces, and heat is released. The unburned gas and the burned gas satisfy the euler equation:
Figure BDA0003153637430000031
where ρ, u, v, p, E represent density, x-direction movement velocity, y-direction movement velocity, pressure, and total energy per unit mass of gas, respectively. Total energy represents the sum of internal energy and kinetic energy, i.e.
Figure BDA0003153637430000032
Where e represents the internal energy per mass of gas. The state equations of the unburned gas and the burned gas satisfy:
p＝(γ-1)ρ(e-q), (2)
where γ represents the specific heat ratio and q represents the heat parameter released due to the chemical reaction.
Step 3.2: setting the boundary conditions of the control equation set of the compressible reaction fluid.
In the construction process of the high-order precision numerical flux, physical quantities to surrounding calculation units need to be used. Therefore, in order to update the conservative value near the boundary of the calculation region, the calculation region needs to be extended by three layers.
Left and right borders:
for the incoming flow boundary condition, the conservative value of the continuation unit is set as:
Figure BDA0003153637430000033
where i is 1,2,3, j is 0,1,2, …, N, U is (ρ, ρ U, ρ v, ρ E)T。U-1,j、U-2,j、U-3,jRespectively representing the conservative quantities, U, of the first, second and third layers of the left continuation of the left boundary of the calculation regionM+1,j、UM+2,j、UM+3,jRespectively representing the conservative quantities of the first layer, the second layer and the third layer which extend to the right from the right boundary of the calculation region,
Figure BDA0003153637430000034
the specific numerical value of (a) needs to be set according to actual working conditions.
For the outflow boundary conditions, the conservative amount of the continuation unit is set as:
Figure BDA0003153637430000035
where i is 1,2,3, j is 0,1,2, …, N. That is, the layer of conservative value closest to the left boundary or the right boundary in the calculation region is directly assigned to the conservative value of the continuation unit.
For the fixed wall boundary condition, the conservation quantity of the continuation unit is set as follows:
Figure BDA0003153637430000041
where i is 1,2,3, j is 0,1,2, …, N. Respectively assigning the mass, the y-direction momentum and the total energy of the first layer, the second layer and the third layer which are closest to the left boundary or the right boundary in the calculation region to corresponding physical quantities of the first layer, the second layer and the third layer which extend outwards; and assigning opposite numbers of x-direction momentums of the first, second and third layers closest to the left or right boundary in the calculation region to corresponding physical quantities of the first, second and third layers extended outward, respectively.
Upper and lower boundaries:
for the incoming flow boundary condition, the conservative value of the continuation unit is set as:
Figure BDA0003153637430000042
where i is 0,1,2, …, M, j is 1,2, 3. U shapei,N+1、Ui,N+2、Ui,N+3Respectively representing the conservative values, U, of the first, second and third layers of the upper boundary continuation of the calculation regioni,-1、Ui,-2、Ui,-3Respectively representing the conservative quantities of the first layer, the second layer and the third layer which extend downwards from the lower boundary of the calculation region,
Figure BDA0003153637430000043
the specific numerical value of (a) needs to be set according to actual working conditions.
For the outflow boundary conditions, the conservative amount of the continuation unit is set as:
Figure BDA0003153637430000044
where i is 0,1,2, …, M, j is 1,2, 3. Namely, the layer of conservation quantity closest to the upper boundary or the lower boundary in the calculation region is directly assigned to the conservation quantity of the continuation unit.
For a solid wall boundary, the conservation quantity of the continuation unit is set as follows:
Figure BDA0003153637430000045
where i is 0,1,2, …, M, j is 1,2, 3. Respectively assigning the mass, the x-direction momentum and the total energy of the first layer, the second layer and the third layer which are closest to the upper boundary or the lower boundary in the calculation region to corresponding physical quantities of the first layer, the second layer and the third layer which extend outwards; and assigning opposite numbers of y-direction momentums of the first, second and third layers closest to the upper boundary or the lower boundary in the calculation region to corresponding physical quantities of the first, second and third layers extended toward the upper boundary or the lower boundary, respectively.
And 4, step 4: the time step for updating the unburned gas and burned gas conservation quantities is calculated.
Calculating the time step satisfying the CFL condition:
Figure BDA0003153637430000046
wherein CFL represents CFL coefficient, and the value range is (0, 1); Δ x and Δ y represent spatial step sizes in the x-direction and the y-direction, respectively; u. ofi,jAnd vi,jRespectively representing the moving speeds of the fluids in the computing units (i, j) along the x direction and the y direction; c. Ci,jRepresenting the speed of sound of the fluid in the calculation unit (i, j).
And 5: and updating the level set function according to the advection equation, performing reinitialization, and solving the evolution process of the detonation or deflagration interface.
The level set function is updated according to the advection equation as shown below:
φt+μφx+νφy＝0, (4)
where μ and ν represent the velocities of motion of the level set function along the x-direction and the y-direction, respectively, need to be defined by the following system of equations:
Figure BDA0003153637430000051
wherein D represents the velocity of movement of the detonation or deflagration interface, NxAnd NyRepresenting the components of the unit normal direction vector N along the x-direction and the y-direction, respectively.
The third-order TVD Runge-Kutta method is applied to carry out time dispersion on the equation (4), and partial derivatives phi of the level set function in the x direction and the y directionxAnd phiyIt needs to be obtained according to the modified Godunov method. The spatial derivative phi can be determined by the following formulaxAnd phiy：
Figure BDA0003153637430000052
Wherein
Figure BDA0003153637430000053
S (phi) represents the sign function of the level set function phi.
Figure BDA0003153637430000054
Wherein
Figure BDA0003153637430000055
To improve the computational efficiency, the level set function need only be updated by equation (4) in one to four layers of computational cells near the interface.
In order for the level set function to maintain the properties of the distance function, it needs to be reinitialized by the following equation:
Figure BDA0003153637430000056
step 6: the geometric parameters of the detonation or deflagration interface are defined in terms of the values of the level set function.
Let Γ denote the detonation or deflagration interface,
Figure BDA0003153637430000057
the left boundary, the right boundary, the lower boundary and the upper boundary of the calculation unit (i, j) are respectively represented by the proportion of the interface gamma cut, and specific numerical values can be obtained through operator operation calculation of the level set function. For fluids with φ > 0, A is defined by the following equation:
Figure BDA0003153637430000061
Figure BDA0003153637430000062
wherein
Figure BDA0003153637430000063
For fluids with phi < 0, the interfacial cutting ratio A is given by equation A-＝1-A+And obtaining the interface cutting ratio geometrical parameters corresponding to the fluid with phi less than 0 and phi more than 0 respectively by superscript-and-respectively.
The volume fraction α of each fluid can also be calculated by operator operation of the level set function, α being defined by the following equation for fluids with φ > 0:
Figure BDA0003153637430000064
wherein
Figure BDA0003153637430000065
For fluids with phi < 0, the volume fraction alpha is then passed through the equation alpha-＝1-α+The volume fraction geometric parameters corresponding to fluids where the superscripts-and + represent φ < 0 and φ > 0, respectively, are evaluated.
And 7: a virtual fluid is constructed in a computational cell near the interface according to a continuation equation.
And constructing the physical quantity state of the virtual fluid according to the continuation equation:
Vτ±N·▽V＝0, (12)
wherein V represents physical quantity needing continuation, including density, moving speed in x direction, moving speed in y direction and pressure; n represents a unit normal direction vector along the interface, by equation
Figure BDA0003153637430000071
Obtaining; τ denotes an artificial time step, defined as:
Figure BDA0003153637430000072
in equation (12), + sign is used to calculate the physical quantity of the fluid in the virtual cell with φ > 0; the sign is used to calculate the physical quantity of the fluid in the virtual cell for φ < 0.
Performing time dispersion on equation (12) by using a third-order TVD Runge-Kutta method, and calculating a spatial derivative by using a windward format
Figure BDA0003153637430000073
And
Figure BDA0003153637430000074
Figure BDA0003153637430000075
Figure BDA0003153637430000076
and 8: in a computing unit near a detonation or deflagration interface, a Riemann problem containing detonation waves or deflagration waves is constructed and solved, and physical quantity states of two sides of the detonation waves or deflagration waves, including density, movement speed, pressure and specific internal energy, are obtained through computing.
In a computing unit near the interface, initial conditions of the Riemann problem are constructed along the French direction according to the states of the real fluid and the virtual fluid. According to the mass, momentum and energy conservation equation, the fluid after the wave front and the wave back of the detonation wave or the deflagration wave satisfies the following relational expression:
Figure BDA0003153637430000077
where subscripts 0 and 1 denote the corresponding physical quantities of unburned gas and post-burned gas before or after the detonation wave or deflagration wave, respectively, and w denotes the velocity of movement of the fluid relative to the detonation wave or deflagration wave. The fluid state after the shock wave front and the shock wave back needs to satisfy a shock wave relational expression; the fluid state after the wave front and the wave back of the sparse wave needs to satisfy that the Riemann invariant is constant.
When CJ detonation occurs, Riemann's solution is composed of non-reaction wave (shock wave or rarefaction wave), contact discontinuity, rarefaction wave and CJ detonation wave; when strong detonation occurs, the riemann solution consists of an unreactive wave, a contact discontinuity and a strong detonation wave. When deflagration occurs, the Riemann solution consists of an unreactive wave, a contact discontinuity, a deflagration wave, and an unreactive wave. And solving the physical quantity states of two sides of each wave system, including density, motion speed, pressure and specific internal energy, according to the relation of shock waves, rarefaction waves, contact discontinuities and detonation waves or deflagration waves.
Physical quantity states of two sides of each wave system in Riemann solution containing detonation waves or deflagration waves are solved through an iterative method, preferably, a dichotomy is selected, and although the dichotomy solving process is relatively complicated, the dichotomy solving process has strong robustness.
And step 9: the conservative exchange terms of the unburned gas and the burned gas along the interface are constructed according to Riemann's solution containing the detonation wave and the deflagration wave.
For the detonation situation, calculating the wave velocity of the detonation wave according to Riemann's solution containing the detonation wave:
Figure BDA0003153637430000081
where ρ isrAnd urRespectively representing the density and the speed of the unburned gas before the detonation wave;
Figure BDA0003153637430000082
and
Figure BDA0003153637430000083
respectively, the density and velocity of the detonation wave post-burned gas.
For unburned gas, the interface term is calculated by:
Figure BDA0003153637430000084
where ρ is1、u1、v1、p1、E1Respectively representing the density of the unburned gas before the detonation wave, the moving speed in the x direction, the moving speed in the y direction, the pressure and the total energy of unit mass of the gas;
Figure BDA0003153637430000085
represents the length of the interface truncation in the computing unit (i, j),the following equation is used to obtain:
Figure BDA0003153637430000086
the interface term of the burnt gas directly takes the opposite value of the interface term of the unburnt gas, thereby ensuring that the conservation of mass, momentum and energy is satisfied along the detonation interface.
The detonation situation is slightly different from the detonation, and because the Riemannian solution containing the detonation wave consists of four waves and one more non-reaction wave is added to the side of the detonation wave, the detonation problem can be solved under the condition that the wave velocity of the detonation wave is known.
For unburned gas, the interface term can be calculated by:
Figure BDA0003153637430000087
where ρ is0、u0、p0、E0Respectively representing the density, the fluid motion speed, the pressure and the total energy of unit mass gas of the unburned gas in the deflagration wave front in Riemann solution; u. of1、v1、urRespectively representing the x-direction movement speed, the y-direction movement speed and the normal-direction movement speed of the unburned gas in the initial condition of the Riemann problem; d represents the velocity of the detonation wave.
The interface term of the burned gas takes the direct value of the inverse of the interface term of the unburned gas, thereby ensuring that conservation of mass, momentum and energy is satisfied along the deflagration interface.
Step 10: according to the modified finite volume method, a conservation-oriented discrete format of the control equation set is constructed.
In the cutting unit (i, j), the finite volume dispersion of the correction is carried out for the unburned gas and the burned gas respectively, and the discrete form of the control equation system is obtained:
Figure BDA0003153637430000091
wherein the geometric parameter alpha of the interfacei,j、
Figure BDA0003153637430000092
Obtaining according to the step 6; conservative exchange term X (gamma) of fluid along interfacei,j) Obtaining according to the step 7-9; F. g denotes the numerical flux of the calculation cell boundaries in the x-direction and y-direction, respectively.
Step 11: and combining a two-step four-step Lax-Wendroff method and a conservation interface method to realize the high-precision construction of numerical flux of detonation and deflagration problems, substituting a modified finite volume discrete format, and updating the conservation quantities of unburned gases and burned gases. The two-step four-order Lax-Wendroff method can reduce the time advancing step of the four-order precision algorithm by half, thereby obviously improving the calculation efficiency, but is mainly applied to the single medium fluid motion problem. There is therefore a need to apply the level set method and virtual fluid method, accurately capture the multimedia interface location, and translate the multimedia problem into a single media problem. The conservation interface method is developed based on a combined algorithm of a level set method and a virtual fluid method. The level set method can quickly and accurately process the topological change of the interface and accurately capture the evolution process of the multi-medium interface; the virtual fluid method can skillfully convert the multi-medium problem into the single-medium problem, thereby effectively inhibiting the non-physical oscillation of the multi-medium interface problem; by solving the Riemann problem containing detonation waves or deflagration waves, the conservation interface method can improve the non-conservation defect of the detonation reaction flow problem and improve the reliability of the prediction result of the detonation and deflagration processes. Therefore, the two-step four-order Lax-Wendroff method is required to be combined with the conservation interface method, so that the high-precision prediction of the detonation and deflagration process is realized, and meanwhile, the calculation efficiency of the high-precision prediction of the detonation and deflagration process is obviously improved.
Defining the areas where the unburnt gas and the burnt gas are located by a level set method by utilizing the step 5; with step 7, the multi-media detonation and detonation problems are converted into single-media problems with unburned gases and burned gases by a virtual fluid approach.
Step 11.1: according to calculationAverage conservation of energy within a cell
Figure BDA0003153637430000093
And calculating the conservation quantities at the cell boundaries
Figure BDA0003153637430000094
And
Figure BDA0003153637430000095
obtaining time derivatives of variables using HWENO reconstruction techniques and generalized Riemann solvers
Figure BDA0003153637430000096
And
Figure BDA0003153637430000097
step 11.2: calculating the average conservative value after half time step
Figure BDA0003153637430000101
And calculating the conservation quantities at the cell boundaries
Figure BDA0003153637430000102
And
Figure BDA0003153637430000103
obtaining time derivatives of the variables again using HWENO reconstruction techniques and generalized Riemann solvers
Figure BDA0003153637430000104
And
Figure BDA0003153637430000105
Figure BDA0003153637430000106
Figure BDA0003153637430000107
Figure BDA0003153637430000108
Figure BDA0003153637430000109
step 11.3: and constructing the numerical flux of the boundary of the calculation unit according to the time derivative of the numerical flux, and substituting the numerical flux into a finite volume discrete format to update the conservation quantities of the unburned gas and the burned gas.
Figure BDA00031536374300001010
Figure BDA00031536374300001011
Figure BDA00031536374300001012
Figure BDA00031536374300001013
Figure BDA00031536374300001014
Figure BDA00031536374300001015
Step 12: judging the calculated time t of the nth step according to the time step delta tnIn relation to the termination time T, if Tn+ Δ T ≦ T, then Tn+1＝tn+ Δ t, return to step 3, update the conservation quantity U at the next momentn+1(ii) a Else Δt＝T-tnAnd returning to the step 3, obtaining the conservative quantity of the termination time, calculating the mass of the termination time, the movement speed in the x direction, the movement speed in the y direction and the pressure according to the state equation, and ending the cycle.
Outputting a level set function, mass, moving speed in the x direction, moving speed in the y direction and pressure of all the computing units at the termination moment, predicting the position of a detonation or deflagration interface and the areas where unburnt gas and burnt gas are located, and computing the mass, the moving speed in the x direction, the moving speed in the y direction and the pressure of fluid in the areas.
Further comprising step 13: predicting detonation and deflagration problems by utilizing the steps 1 to 12, and improving the calculation efficiency of the high-precision prediction method; meanwhile, the problem of non-conservation is effectively solved, the detonation and deflagration processes are reliably predicted, and the technical problem of related engineering in the field of explosion mechanics is solved.
The related engineering technical problems in the field of explosion mechanics comprise weapon ammunition design and damage assessment and disaster prevention and reduction of explosion accidents.
Has the advantages that:
1. the invention discloses a high-precision prediction method for explosion reaction flow based on a generalized Riemannian solution device, which is used for solving the numerical flux of a calculation unit boundary by using a two-step and four-order Lax-Wendroff method based on the generalized Riemannian solution device, wherein the two-step and four-order Lax-Wendroff method integrates the characteristics of simplicity of Runge-Kutta algorithm and the space-time coupling property of the generalized Riemannian solution device, so that the prediction method can improve the calculation efficiency and stability of the prediction detonation and detonation process.
2. The invention discloses a high-precision prediction method for explosion reaction flow based on a generalized Riemannian solver, which is characterized in that a conservation type calculation format is obtained by applying a conservation interface method, a modified finite volume method is adopted to discretely control an equation, and conservation quantity exchange terms are solved by solving the Riemannian problem containing detonation waves or deflagration waves, so that the conservation of mass, momentum and energy along an interface can be met. Therefore, the prediction method can overcome the defect that the traditional level set/virtual fluid method is lack of conservation, and improve the accuracy of the detonation and deflagration process prediction.
3. The invention discloses a high-precision prediction method for an explosion reaction flow based on a generalized Riemannian solution device, which is characterized in that a level set method is applied to track a detonation or deflagration interface, the level set method can efficiently and quickly realize high-order precision, and meanwhile, the high-dimensional space can be simply and effectively expanded. Therefore, the prediction method can quickly and accurately capture the topological change process of the detonation or deflagration interface.
4. The invention discloses a high-precision prediction method for explosion reaction flow based on a generalized Riemannian solver, which is used for constructing states of unburned gas and burned gas in a virtual unit by using a virtual fluid method and can skillfully convert a multi-medium problem into a single-medium problem. Therefore, the invention can effectively inhibit non-physical oscillation near the multi-medium interface.
Drawings
FIG. 1 is a schematic diagram of a modified finite volume method;
FIG. 2 is a schematic view of a virtual fluid process;
FIG. 3 is a schematic diagram of Riemann's structure containing detonation waves (left diagram: CJ detonation, right diagram: strong detonation);
FIG. 4 is a schematic view of a Riemann solution structure including a deflagration wave;
FIG. 5 is a Riemannian solution profile (density, velocity, and pressure profiles, from left to right) of an example along a detonation interface at an initial time;
fig. 6 is a flowchart of a high-precision prediction method for an explosion reaction flow based on a generalized riemann solver disclosed in this embodiment.
Detailed Description
For better illustrating the objects and advantages of the present invention, the present invention will be further described with reference to the accompanying drawings and examples.
Example 1:
the embodiment is applied to the double explosion wave fusion algorithm.
As shown in fig. 6, the embodiment discloses a high-precision prediction method for an explosion reaction flow based on a generalized riemann solver, and the specific implementation method is as follows:
step 1: computing regionsSet to (0,2m) × (0,2m), and the grid size is set to
Figure BDA0003153637430000121
I.e. 400 x 400 grids are set within the calculation area. The calculation units are denoted as (i, j), where i is 0,1, …,399, j is 0,1, …, 399.
Step 2: the level set function at the initial time is taken as:
Figure BDA0003153637430000122
the densities, the moving speed in the x direction, the moving speed in the y direction, and the pressures of the unburned gas and the burned gas at the initial time are taken as follows:
ρ u v p
unburned gas 1.0kg/m3 0m/s 0m/s 1.0×105Pa
Burnt gas 0.142168kg/m3 0m/s 0m/s 9.45695×104Pa
Table 1
And step 3: both the unburned gas and the burned gas satisfy the euler equation, i.e., equation set (1). To close the system of control equations, the state equations for the unburned gases as well as the burned gases are given, equation (2). And the parameters of the state equation are taken as follows:
γ q
unburned gas 1.4 2.0×106J/kg
Burnt gas 1.4 0
Table 2
The left and right boundaries and the upper and lower boundaries are set as the fixed wall boundary conditions:
Figure BDA0003153637430000123
where i is 1,2,3, j is 0,1,2, …, 399.
Figure BDA0003153637430000131
Where i is 0,1,2, …,399, j is 1,2, 3.
And 4, step 4: the time step is calculated according to equation (3) where the CFL coefficient takes 0.4.
Then, step 5 to step 10 are performed, wherein the level set function updating process in step 5, the riemann solution solving process in step 7, and the conservative exchange term solving process in step 9 all require the known propagation velocity of the deflagration wave, and the deflagration wave velocity is calculated according to the following formula:
Figure BDA0003153637430000132
where V, p, ρ are the velocity, pressure and density of the unburned gas in the normal direction, respectively.
In step 8, a dichotomy is applied to iteratively solve the riemann problem containing the deflagration wave, and the density, the speed and the pressure of the deflagration wave front unburned gas in the initial time riemann solution can be calculated and obtained as follows:
ρ0 u0 p0
1.157kg/m3 55.594m/s 1.227×105Pa
table 3
Fig. 5 is a riemann solution distribution graph along a deflagration interface at an initial time, sequentially from left to right, density, velocity and pressure distribution graphs, and it can be seen that the riemann solution is composed of four waves, sequentially from left to right: shock, contact discontinuity, deflagration wave, and shock.
And then, step 11 is carried out, a two-step fourth-order Lax-Wendroff method is used for solving and calculating the unit boundary numerical flux, and the unit boundary numerical flux is substituted into a finite volume discrete format to update the conservation quantities of the unburned gas and the burned gas. And combining a two-step four-step Lax-Wendroff method and a conservation interface method to realize the high-precision construction of numerical flux of detonation and deflagration problems, substituting a modified finite volume discrete format, and updating the conservation quantities of unburned gases and burned gases. The two-step four-order Lax-Wendroff method can reduce the time advancing step of the four-order precision algorithm by half, thereby obviously improving the calculation efficiency, but is mainly applied to the single medium fluid motion problem. There is therefore a need to apply the level set method and virtual fluid method, accurately capture the multimedia interface location, and translate the multimedia problem into a single media problem. The conservation interface method is developed based on a combined algorithm of a level set method and a virtual fluid method. The level set method can quickly and accurately process the topological change of the interface and accurately capture the evolution process of the multi-medium interface; the virtual fluid method can skillfully convert the multi-medium problem into the single-medium problem, thereby effectively inhibiting the non-physical oscillation of the multi-medium interface problem; by solving the Riemann problem containing detonation waves or deflagration waves, the conservation interface method can improve the non-conservation defect of the detonation reaction flow problem and improve the reliability of the prediction result of the detonation and deflagration processes. Therefore, the two-step four-order Lax-Wendroff method is required to be combined with the conservation interface method, so that the high-precision prediction of the detonation and deflagration process is realized, and meanwhile, the calculation efficiency of the high-precision prediction of the detonation and deflagration process is obviously improved.
And then, step 12 to step 13 are carried out, and finally, the level set function, the mass, the movement speed in the x direction, the movement speed in the y direction and the pressure in the calculation area at different moments are output.
And (4) analyzing a calculation result:
the double deflagration waves initially propagate outward each other, and as time progresses, the two bubbles gradually begin to merge. The invention discloses a high-precision prediction method for explosion reaction flow based on a generalized Riemannian solver, which combines a two-step four-order Lax-Wendroff method and a conservation interface method, can accurately describe the topological change of a deflagration interface, can predict the fluid motion process of unburned gas and burned gas with high precision, and can effectively improve the calculation efficiency.
The above detailed description is intended to illustrate the objects, aspects and advantages of the present invention, and it should be understood that the above detailed description is only exemplary of the present invention and is not intended to limit the scope of the present invention, and any modifications, equivalents, improvements and the like made within the spirit and principle of the present invention should be included in the scope of the present invention.

Claims (10)

1. A high-precision prediction method for explosion reaction flow based on a generalized Riemannian solution device is characterized by comprising the following steps: comprises the following steps of (a) carrying out,
step 1: establishing a rectangular coordinate system, defining space steps in the x direction and the y direction, and arranging (M +1) × (N +1) calculation units in a calculation region, wherein i is 0,1, …, M, j is 0,1, …, N;
step 2: defining level set functions of all calculation units in a calculation area at the initial moment and fluid physical quantities of unburned gas and burned gas, including density, moving speed in the x direction, moving speed in the y direction and pressure according to initial conditions;
and step 3: establishing a system of control equations for the compressible reaction fluid and setting boundary conditions;
and 4, step 4: calculating and updating the time step of the unburnt gas and burnt gas conservation quantity;
and 5: updating a level set function according to an advection equation, performing reinitialization, and solving an evolution process of a detonation or deflagration interface;
step 6: defining the geometric parameters of the detonation or deflagration interface according to the numerical value of the level set function;
and 7: constructing a virtual fluid in a computing unit near the interface according to a continuation equation;
and 8: in a computing unit near a detonation or deflagration interface, constructing and solving a Riemann problem containing a detonation wave or a deflagration wave, and computing to obtain physical quantity states of two sides of the detonation wave or the deflagration wave, wherein the physical quantity states comprise density, motion speed, pressure and specific internal energy;
and step 9: constructing a conservation quantity exchange term of unburned gas and burned gas along an interface according to Riemannian decomposition containing detonation waves and deflagration waves;
step 10: constructing a conservation-type discrete format of a control equation set according to a modified finite volume method;
step 11: combining a two-step four-order Lax-Wendroff method and a conservation interface method, realizing the high-precision construction of numerical flux of detonation and deflagration problems, substituting a modified finite volume discrete format, and updating the conservation quantities of unburned gases and burned gases; the two-step four-order Lax-Wendroff method can reduce the time advancing step of the four-order precision algorithm by half, thereby obviously improving the calculation efficiency, but is mainly applied to the problem of single medium fluid movement; there is therefore a need to apply level set methods and virtual fluid methods to accurately capture multi-media interface locations and translate multi-media problems into single-media problems; the conservation interface method is developed based on a combined algorithm of a level set method and a virtual fluid method; the level set method can quickly and accurately process the topological change of the interface and accurately capture the evolution process of the multi-medium interface; the virtual fluid method can skillfully convert the multi-medium problem into the single-medium problem, thereby effectively inhibiting the non-physical oscillation of the multi-medium interface problem; by solving the Riemann problem containing detonation waves or deflagration waves, the conservation interface method can improve the non-conservation defect of the detonation reaction flow problem and improve the reliability of the prediction result of the detonation and deflagration processes; therefore, the two-step four-order Lax-Wendroff method and the conservation interface method are combined, so that the high-precision prediction of the detonation and the deflagration process is realized, and the calculation efficiency of the high-precision prediction of the detonation and the deflagration process is obviously improved;
step 12: judging the calculated time t of the nth step according to the time step delta tnIn relation to the termination time T, if Tn+ Δ T ≦ T, then Tn+1＝tn+ Δ t, return to step 3, update the conservation quantity U at the next momentn+1(ii) a Otherwise, T-TnReturning to the step 3, obtaining the conservation quantity of the termination time, calculating the mass of the termination time, the movement speed in the x direction, the movement speed in the y direction and the pressure according to the state equation, and ending the cycle;
outputting a level set function, mass, moving speed in the x direction, moving speed in the y direction and pressure of all the computing units at the termination moment, predicting the position of a detonation or deflagration interface and the areas where unburnt gas and burnt gas are located, and computing the mass, the moving speed in the x direction, the moving speed in the y direction and the pressure of fluid in the areas.
2. The method for predicting the explosion reaction flow with high precision based on the generalized Riemannian solver according to claim 1, wherein the method comprises the following steps: further comprising step 13: predicting detonation and deflagration problems by utilizing the steps 1 to 12, and improving the calculation efficiency of the high-precision prediction method; meanwhile, the problem of non-conservation is effectively solved, the reliable prediction of detonation and deflagration processes is realized, and the related engineering technical problem in the field of explosion mechanics is further solved;
the related engineering technical problems in the field of explosion mechanics comprise weapon ammunition design and damage assessment and disaster prevention and reduction of explosion accidents.
3. The method for predicting the explosion reaction flow with high precision based on the generalized Riemannian solver according to claim 1 or 2, wherein: the step 2 is realized by the method that,
in the level set method, the level set function is a symbol distance function, denoted by phi; let regions where + 0 and + 0 denote distribution regions of two fluids, i.e., unburned gas and burned gas, respectively, and a position where + 0 denotes a detonation or deflagration interface; at the initial time, the value of the level set function φ is defined using the following equation:
Figure FDA0003153637420000021
wherein gamma is0A distribution curve representing the interface at the initial moment, d (x, y) representing the center (x, y) of the computing unit to the initial interface Γ0Distance of (d), omega1、Ω2Areas showing distribution of unburned gas and burned gas, respectively;
according to the initial conditions, the density, the moving speed in the x direction, the moving speed in the y direction and the pressure of unburned gas and burned gas on both sides of the initial interface are defined.
4. The method for predicting the explosion reaction flow with high precision based on the generalized Riemannian solver according to claim 3, wherein the method comprises the following steps: the step 3 is realized by the method that,
step 3.1: establishing a system of control equations for the compressible reactant fluid;
establishing a control equation system of the compressible fluid by applying a CJ model, namely, considering chemical reaction rate and detonation and deflagration interfaces as jumping discontinuities, immediately converting chemical reactants into products along the interfaces and releasing heat; the unburned gas and the burned gas satisfy the euler equation:
Figure FDA0003153637420000031
wherein ρ, u, v, p, E represent density, x-direction movement velocity, y-direction movement velocity, pressure, and total energy per unit mass of gas, respectively; total energy represents the sum of internal energy and kinetic energy, i.e.
Figure FDA0003153637420000032
Wherein e represents the internal energy per unit mass of gas; the state equations of the unburned gas and the burned gas satisfy:
p＝(γ-1)ρ(e-q), (2)
wherein γ represents the specific heat ratio and q represents the heat released due to the chemical reaction parameter;
step 3.2: setting boundary conditions of a control equation set of the compressible reaction fluid;
in the construction process of the high-order precision numerical flux, physical quantities to surrounding computing units are required to be used; therefore, in order to update the conservative value near the boundary of the calculation region, the calculation region needs to be extended by three layers;
left and right borders:
for the incoming flow boundary condition, the conservative value of the continuation unit is set as:
Figure FDA0003153637420000033
where i is 1,2,3, j is 0,1,2, …, N, U is (ρ, ρ U, ρ v, ρ E)T；U-1,j、U-2,j、U-3,jRespectively representing the conservative quantities, U, of the first, second and third layers of the left continuation of the left boundary of the calculation regionM+1,j、UM+2,j、UM+3,jRespectively representing the conservative quantities of the first layer, the second layer and the third layer which extend to the right from the right boundary of the calculation region,
Figure FDA0003153637420000034
the specific numerical value of (2) needs to be set according to the actual working condition;
for the outflow boundary conditions, the conservative amount of the continuation unit is set as:
Figure FDA0003153637420000035
wherein i is 1,2,3, j is 0,1,2, …, N; directly assigning the layer of conservation quantity closest to the left boundary or the right boundary in the calculation region to the conservation quantity of the continuation unit;
for the fixed wall boundary condition, the conservation quantity of the continuation unit is set as follows:
Figure FDA0003153637420000036
wherein i is 1,2,3, j is 0,1,2, …, N; respectively assigning the mass, the y-direction momentum and the total energy of the first layer, the second layer and the third layer which are closest to the left boundary or the right boundary in the calculation region to corresponding physical quantities of the first layer, the second layer and the third layer which extend outwards; assigning the opposite numbers of the x-direction momentum of the first, second and third layers closest to the left boundary or the right boundary in the calculation region to the corresponding physical quantities of the first, second and third layers extended outwards respectively;
upper and lower boundaries:
for the incoming flow boundary condition, the conservative value of the continuation unit is set as:
Figure FDA0003153637420000041
wherein i is 0,1,2, …, M, j is 1,2, 3; u shapei,N+1、Ui,N+2、Ui,N+3Respectively representing the conservative values, U, of the first, second and third layers of the upper boundary continuation of the calculation regioni,-1、Ui,-2、Ui,-3Respectively representing the conservative quantities of the first layer, the second layer and the third layer which extend downwards from the lower boundary of the calculation region,
Figure FDA0003153637420000042
the specific numerical value of (2) needs to be set according to the actual working condition;
for the outflow boundary conditions, the conservative amount of the continuation unit is set as:
Figure FDA0003153637420000043
wherein i is 0,1,2, …, M, j is 1,2, 3; directly assigning the layer of conservation quantity closest to the upper boundary or the lower boundary in the calculation region to the conservation quantity of the continuation unit;
for a solid wall boundary, the conservation quantity of the continuation unit is set as follows:
Figure FDA0003153637420000044
wherein i is 0,1,2, …, M, j is 1,2, 3; respectively assigning the mass, the x-direction momentum and the total energy of the first layer, the second layer and the third layer which are closest to the upper boundary or the lower boundary in the calculation region to corresponding physical quantities of the first layer, the second layer and the third layer which extend outwards; and assigning opposite numbers of y-direction momentums of the first, second and third layers closest to the upper boundary or the lower boundary in the calculation region to corresponding physical quantities of the first, second and third layers extended toward the upper boundary or the lower boundary, respectively.
5. The method for predicting the explosion reaction flow with high precision based on the generalized Riemannian solver according to claim 4, wherein the method comprises the following steps: step 4, the method is realized by the following steps,
calculating the time step satisfying the CFL condition:
Figure FDA0003153637420000045
wherein CFL represents CFL coefficient, and the value range is (0, 1); Δ x and Δ y represent spatial step sizes in the x-direction and the y-direction, respectively; u. ofi,jAnd vi,jRespectively representing the moving speeds of the fluids in the computing units (i, j) along the x direction and the y direction; c. Ci,jRepresenting the speed of sound of the fluid in the calculation unit (i, j).
6. The method for predicting the explosion reaction flow with high precision based on the generalized Riemannian solver according to claim 5, wherein the method comprises the following steps: step 5 the method is realized by the following steps,
the level set function is updated according to the advection equation as shown below:
φt+μφx+νφy＝0, (4)
where μ and ν represent the velocities of motion of the level set function along the x-direction and the y-direction, respectively, need to be defined by the following system of equations:
Figure FDA0003153637420000051
wherein D represents the velocity of movement of the detonation or deflagration interface, NxAnd NyRepresenting the components of the unit normal direction vector N along the x-direction and the y-direction, respectively;
the third-order TVD Runge-Kutta method is applied to carry out time dispersion on the equation (4), and partial derivatives phi of the level set function in the x direction and the y directionxAnd phiyThe method needs to be solved according to a corrected Godunov method; the spatial derivative phi can be determined by the following formulaxAnd phiy：
Figure FDA0003153637420000052
Wherein
Figure FDA0003153637420000053
S (φ) represents a sign function of a level set function φ;
Figure FDA0003153637420000054
wherein
Figure FDA0003153637420000055
In order to improve the calculation efficiency, the level set function is updated only in one to four layers of calculation units near the interface through the equation (4);
in order for the level set function to maintain the properties of the distance function, it needs to be reinitialized by the following equation:
Figure FDA0003153637420000056
7. the method for predicting the explosion reaction flow with high precision based on the generalized Riemannian solver according to claim 6, wherein the method comprises the following steps: step 6 is realized by the method that,
let Γ denote the detonation or deflagration interface,
Figure FDA0003153637420000057
the left boundary, the right boundary, the lower boundary and the upper boundary of the computing unit (i, j) are respectively represented by the proportion of the interface gamma cut, and specific numerical values can be obtained through operator operation of a level set function; for fluids with φ > 0, A is defined by the following equation:
Figure FDA0003153637420000061
Figure FDA0003153637420000062
wherein
Figure FDA0003153637420000063
For fluids with phi < 0, the interfacial cutting ratio A is given by equation A-＝1-A+Obtaining the interface cutting ratio geometrical parameters corresponding to the fluid with phi less than 0 and phi more than 0 respectively by superscript-and-;
the volume fraction α of each fluid can also be calculated by operator operation of the level set function, α being defined by the following equation for fluids with φ > 0:
Figure FDA0003153637420000064
wherein
Figure FDA0003153637420000065
For fluids with phi < 0, the volume fraction alpha is then passed through the equation alpha-＝1-α+The volume fraction geometric parameters corresponding to fluids where the superscripts-and + represent φ < 0 and φ > 0, respectively, are evaluated.
8. The method for predicting the explosion reaction flow with high precision based on the generalized Riemannian solver according to claim 7, wherein the method comprises the following steps: step 7 is realized by the method that,
and constructing the physical quantity state of the virtual fluid according to the continuation equation:
Figure FDA0003153637420000066
wherein V represents physical quantity needing continuation, including density, moving speed in x direction, moving speed in y direction and pressure; n represents a unit normal direction vector along the interface, by equation
Figure FDA0003153637420000071
Obtaining; τ denotes an artificial time step, defined as:
Figure FDA0003153637420000072
in equation (12), + sign is used to calculate the physical quantity of the fluid in the virtual cell with φ > 0; -the number is used to calculate the physical quantity of the fluid in the virtual cell with phi < 0;
performing time dispersion on equation (12) by using a third-order TVD Runge-Kutta method, and calculating a spatial derivative by using a windward format
Figure FDA0003153637420000073
And
Figure FDA0003153637420000074
Figure FDA0003153637420000075
Figure FDA0003153637420000076
step 8 is realized by a method comprising the following steps of,
constructing initial conditions of the Riemann problem along a French direction according to states of the real fluid and the virtual fluid in a computing unit near the interface; according to the mass, momentum and energy conservation equation, the fluid after the wave front and the wave back of the detonation wave or the deflagration wave satisfies the following relational expression:
Figure FDA0003153637420000077
wherein subscripts 0 and 1 denote corresponding physical quantities of unburned gas and post-burned gas before or after a detonation wave or a deflagration wave, respectively, and w denotes a moving velocity of a fluid with respect to the detonation wave or the deflagration wave; the fluid state after the shock wave front and the shock wave back needs to satisfy a shock wave relational expression; the fluid state after the wave front and the wave back of the sparse wave needs to meet the condition that the Riemann invariant is a constant;
when CJ detonation occurs, Riemann's solution is composed of non-reaction wave, contact discontinuity, rarefaction wave and CJ detonation wave; when strong detonation occurs, Riemann's solution is composed of non-reaction wave, contact discontinuity and strong detonation wave; when deflagration occurs, Riemann's solution is composed of non-reaction wave, contact discontinuity, deflagration wave and non-reaction wave; according to the relation of shock waves, rarefaction waves, contact discontinuities and detonation waves or deflagration waves, solving the physical quantity states of two sides of each wave system, including density, movement speed, pressure and specific internal energy;
step 9 is implemented by a method comprising the following steps,
for the detonation situation, calculating the wave velocity of the detonation wave according to Riemann's solution containing the detonation wave:
Figure FDA0003153637420000081
where ρ isrAnd urRespectively representing the density and the speed of the unburned gas before the detonation wave;
Figure FDA0003153637420000082
and
Figure FDA0003153637420000083
respectively representing the density and the speed of the post-combustion gas after the detonation wave;
for unburned gas, the interface term is calculated by:
Figure FDA0003153637420000084
where ρ is1、u1、v1、p1、E1Respectively representing the density of the unburned gas before the detonation wave, the moving speed in the x direction, the moving speed in the y direction, the pressure and the total energy of unit mass of the gas;
Figure FDA0003153637420000085
the length of the interface cut in the calculation unit (i, j) is represented by the following equation:
Figure FDA0003153637420000086
the interface item of the burnt gas is directly taken as the opposite number of the interface item of the unburned gas, thereby ensuring that the conservation of mass, momentum and energy is met along the detonation interface;
the detonation situation is slightly different from the detonation, and because the Riemannian solution containing the detonation wave consists of four waves and one more non-reaction wave is added on the side of the detonation wave, the detonation problem can be solved under the condition that the wave velocity of the detonation wave is known;
for unburned gas, the interface term is calculated by:
Figure FDA0003153637420000087
where ρ is0、u0、p0、E0Respectively representing the density, the fluid motion speed, the pressure and the total energy of unit mass gas of the unburned gas in the deflagration wave front in Riemann solution; u. of1、v1、urRespectively representing the x-direction movement speed, the y-direction movement speed and the normal-direction movement speed of the unburned gas in the initial condition of the Riemann problem; d represents the wave velocity of the detonation wave;
the interface term of the burnt gas is directly taken as the inverse number of the interface term of the unburnt gas, thereby ensuring that the conservation of mass, momentum and energy is satisfied along the deflagration interface;
step 10 is implemented by a method comprising the steps of,
in the cutting unit (i, j), the finite volume dispersion of the correction is carried out for the unburned gas and the burned gas respectively, and the discrete form of the control equation system is obtained:
Figure FDA0003153637420000091
wherein the geometric parameter alpha of the interfacei,j、
Figure FDA0003153637420000092
Obtaining according to the step 6; conservative exchange term X (gamma) of fluid along interfacei,j) Obtaining according to the step 7-9; F. g denotes the numerical flux of the calculation cell boundaries in the x-direction and y-direction, respectively.
9. The method for predicting the explosion reaction flow with high precision based on the generalized Riemannian solver according to claim 8, wherein the method comprises the following steps: step 11 is implemented by a method comprising the following steps,
defining the areas where the unburnt gas and the burnt gas are located by a level set method by utilizing the step 5; converting the problems of multi-medium detonation and deflagration into the problems of single medium of unburnt gas and burnt gas through a virtual fluid method by utilizing the step 7;
step 11.1: according to average conservation quantity in computing unit
Figure FDA0003153637420000093
And calculating the conservation quantities at the cell boundaries
Figure FDA0003153637420000094
And
Figure FDA0003153637420000095
obtaining time derivatives of variables using HWENO reconstruction techniques and generalized Riemann solvers
Figure FDA0003153637420000096
And
Figure FDA0003153637420000097
step 11.2: calculating the average conservative value after half time step
Figure FDA0003153637420000098
And calculating the conservation quantities at the cell boundaries
Figure FDA0003153637420000099
And
Figure FDA00031536374200000910
obtaining time derivatives of the variables again using HWENO reconstruction techniques and generalized Riemann solvers
Figure FDA00031536374200000911
And
Figure FDA00031536374200000912
Figure FDA00031536374200000913
Figure FDA00031536374200000914
Figure FDA00031536374200000915
Figure FDA00031536374200000916
step 11.3: according to the time derivative of the numerical flux, the numerical flux of the boundary of the calculation unit is constructed and substituted into a finite volume discrete format, and the conservation quantities of the unburned gas and the burned gas are updated;
Figure FDA0003153637420000101
Figure FDA0003153637420000102
Figure FDA0003153637420000103
Figure FDA0003153637420000104
Figure FDA0003153637420000105
Figure FDA0003153637420000106
10. the method for predicting the explosion reaction flow with high precision based on the generalized Riemannian solver according to claim 9, wherein: physical quantity states of two sides of each wave system in Riemann solution containing detonation waves or deflagration waves are solved through an iterative method, a dichotomy is selected, and although the dichotomy solving process is relatively complicated, the method has strong robustness.
CN202110771451.5A 2021-07-08 2021-07-08 Explosion reaction flow high-precision prediction method based on generalized Riemann solver Active CN113591345B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110771451.5A CN113591345B (en) 2021-07-08 2021-07-08 Explosion reaction flow high-precision prediction method based on generalized Riemann solver

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110771451.5A CN113591345B (en) 2021-07-08 2021-07-08 Explosion reaction flow high-precision prediction method based on generalized Riemann solver

Publications (2)

Publication Number Publication Date
CN113591345A true CN113591345A (en) 2021-11-02
CN113591345B CN113591345B (en) 2024-01-23

Family

ID=78246415

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110771451.5A Active CN113591345B (en) 2021-07-08 2021-07-08 Explosion reaction flow high-precision prediction method based on generalized Riemann solver

Country Status (1)

Country Link
CN (1) CN113591345B (en)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114254573A (en) * 2021-12-17 2022-03-29 北京航空航天大学 High-order processing method for multi-medium interface with symmetric source items
CN116595746A (en) * 2023-05-11 2023-08-15 中国船舶科学研究中心 Method for determining material interface state of strong impact compressible multiphase fluid

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103823916A (en) * 2013-10-23 2014-05-28 沈智军 Arbitrary Lagrange Euler method based on multi-dimensional Riemann solution
CN103971007A (en) * 2014-05-16 2014-08-06 北京航空航天大学 Interface decoupling technology used for processing coupling between compressible fluid and ideal elastic-plastic solid in compression state
US20160038788A1 (en) * 2013-02-06 2016-02-11 Blur Sports Inc. Performance monitoring systems and methods for edging sports
CN106503379A (en) * 2016-10-28 2017-03-15 京工博创（北京）科技有限公司 A kind of gas burst emulation mode that is reacted based on adaptive simplifying with grid subdivision
US20180357347A1 (en) * 2013-02-04 2018-12-13 Comsol Ab Apparatus and method for defining coupled systems on spatial dimensions and extra dimensions
CN109214082A (en) * 2018-09-03 2019-01-15 北京理工大学 A kind of high resolution numerical simulation method of near field underwater blast wave load
CN110750933A (en) * 2019-11-19 2020-02-04 北京理工大学 Accurate interface tracking processing method for coupling Lagrange particles and Euler method

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20180357347A1 (en) * 2013-02-04 2018-12-13 Comsol Ab Apparatus and method for defining coupled systems on spatial dimensions and extra dimensions
US20160038788A1 (en) * 2013-02-06 2016-02-11 Blur Sports Inc. Performance monitoring systems and methods for edging sports
CN103823916A (en) * 2013-10-23 2014-05-28 沈智军 Arbitrary Lagrange Euler method based on multi-dimensional Riemann solution
CN103971007A (en) * 2014-05-16 2014-08-06 北京航空航天大学 Interface decoupling technology used for processing coupling between compressible fluid and ideal elastic-plastic solid in compression state
CN106503379A (en) * 2016-10-28 2017-03-15 京工博创（北京）科技有限公司 A kind of gas burst emulation mode that is reacted based on adaptive simplifying with grid subdivision
CN109214082A (en) * 2018-09-03 2019-01-15 北京理工大学 A kind of high resolution numerical simulation method of near field underwater blast wave load
CN110750933A (en) * 2019-11-19 2020-02-04 北京理工大学 Accurate interface tracking processing method for coupling Lagrange particles and Euler method

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
丁建旭: "爆炸与冲击问题的高精度界面处理及数值研究", 基础科学 *

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114254573A (en) * 2021-12-17 2022-03-29 北京航空航天大学 High-order processing method for multi-medium interface with symmetric source items
CN114254573B (en) * 2021-12-17 2024-05-24 北京航空航天大学 Multi-medium interface high-order processing method with symmetry source item
CN116595746A (en) * 2023-05-11 2023-08-15 中国船舶科学研究中心 Method for determining material interface state of strong impact compressible multiphase fluid

Also Published As

Publication number Publication date
CN113591345B (en) 2024-01-23

Similar Documents

Publication Publication Date Title
CN113591345A (en) Explosion reaction flow high-precision prediction method based on generalized Riemann solution method device
CN111241742A (en) Multiphase flow calculation method
Ha et al. A compressive interface-capturing scheme for computation of compressible multi-fluid flows
CN105512363A (en) Simulation method of water-column separation of pressure pipelines based on Godunov format
CN113792432A (en) Flow field calculation method based on improved FVM-LBFS method
CN115587551B (en) Multi-scale prediction method for ablation behavior of heat-proof structure of hypersonic aircraft
Rokhy et al. Fluid structure interaction with a finite rate chemistry model for simulation of gaseous detonation metal-forming
Guilkey et al. An Eulerian-Lagrangian approach for large deformation fluid structure interaction problems, Part 1: algorithm development
CN116663373A (en) High-precision numerical simulation method suitable for gas phase combustion and explosion
Niu Advection upwinding splitting method to solve a compressible two‐fluid model
Garoosi et al. Presenting a novel higher-order bounded convection scheme for simulation of multiphase flows and convection heat transfer
Wu et al. A new approach to aircraft ditching analysis by coupling free surface lattice Boltzmann and immersed boundary method incorporating surface tension effects
Sharif et al. Numerical study of turbulent natural convection in a side-heated square cavityat various angles of inclination
CN114021499B (en) Aircraft heat protection structure heat conduction calculation method based on FVM-TLBFS method
Witt et al. A study in multiphase modeling of fluidized beds
Narayanan et al. Numerical predictions of the flow and heat transfer characteristics in the film boiling regime during tube quenching
Straka Numerical simulation of heat transfer in packed beds by two population thermal lattice Boltzmann method
Zhang et al. Coupling of level set and volume of fluid methods for simulations of transient internal flow field in solid rocket motors
Balam et al. Numerical solution of natural-convection flow in enclosures: An implicit vorticity boundary condition type method
Irikura et al. Numerical simulation of slugging of stagnant liquid at a V-shaped elbow in a pipeline
Sarmiento et al. A numerical method for shell and thermosyphon heat exchanger analysis
Hao et al. Investigation of oblique water entry of high-speed supercavitating projectiles using transient fluid-structure interaction simulation
CN110705185A (en) Method for predicting pipeline air hammer
Pericleous et al. Free surface Nvier-Stokes flows with simultaneous heat transfer and solidification/melting
Lu et al. Dynamical and thermal performance of molten salt pipe during filling process

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant