阅读说明：本技术 一种基于广义黎曼解法器的爆炸反应流高精度预测方法 (Explosion reaction flow high-precision prediction method based on generalized Riemann solution method device ) 是由 王成 肖敏 昝文涛 杨同会 于 2021-07-08 设计创作，主要内容包括：本发明公开的一种基于广义黎曼解法器的爆炸反应流高精度预测方法,属于爆炸力学技术领域。本发明是通过两步四阶Lax-Wendroff方法以及守恒界面方法的结合实现的。两步四阶Lax-Wendroff方法融合了Runge-Kutta算法的简单性特点以及广义黎曼解法器的时空耦合性质,能够实现数值通量的高精度构造,并且能够将四阶精度算法的时间推进步骤减少一半。守恒界面方法是基于水平集方法和虚拟流体方法的联合算法进行发展的。水平集方法能够快捷且准确地处理界面的拓扑变化,精确地捕捉多介质界面的演化过程；虚拟流体方法能够巧妙地将多介质问题转化为单介质问题,有效抑制多介质界面问题的非物理振荡；守恒界面方法能够改善多介质界面非守恒性问题,提高预测结果的可靠性。(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.)

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 tⁿIn relation to the termination time T, if Tⁿ+ Δ T ≦ T, then Tⁿ⁺¹＝tⁿ+ Δ t, return to step 3, update the conservation quantity U at the next momentⁿ⁺¹(ii) a Otherwise, T-TⁿReturning 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:

wherein gamma is⁰A 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 Γ⁰Distance of (d), omega₁、Ω₂Areas 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:

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.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:

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 region_M+1,j、U_M+2,j、U_M+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,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:

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:

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:

wherein i is 0,1,2, …, M, j is 1,2, 3; u shape_i,N+1、U_i,N+2、U_i,N+3Respectively representing the conservative values, U, of the first, second and third layers of the upper boundary continuation of the calculation region_i,-1、U_i,-2、U_i,-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,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:

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:

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:

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. of_i,jAnd v_i,jRespectively representing the moving speeds of the fluids in the computing units (i, j) along the x direction and the y direction; c. C_i,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:

wherein D represents the velocity of movement of the detonation or deflagration interface, N_xAnd N_yRepresenting 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 direction_xAnd phi_yThe method needs to be solved according to a corrected Godunov method; the spatial derivative phi can be determined by the following formula_xAnd phi_y：

WhereinS (φ) represents a sign function of a level set function φ;

wherein

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:

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,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:

wherein

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:

wherein

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:

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 equationObtaining; τ denotes an artificial time step, defined as: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 formatAnd

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:

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:

where ρ is_rAnd u_rRespectively representing the density and the speed of the unburned gas before the detonation wave;andrespectively representing the density and the speed of the post-combustion gas after the detonation wave;

for unburned gas, the interface term is calculated by:

where ρ is₁、u₁、v₁、p₁、E₁Respectively 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;the length of the interface cut in the calculation unit (i, j) is represented by the following equation:

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:

where ρ is₀、u₀、p₀、E₀Respectively 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. of₁、v₁、u_rRespectively 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:

wherein the geometric parameter alpha of the interface_i,j、Obtaining according to the step 6; conservative exchange term X (gamma) of fluid along interface_i,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 unitAnd calculating the conservation quantities at the cell boundariesAndobtaining time derivatives of variables using HWENO reconstruction techniques and generalized Riemann solversAnd

step 11.2: calculating the average conservative value after half time stepAnd calculating the conservation quantities at the cell boundariesAndobtaining time derivatives of the variables again using HWENO reconstruction techniques and generalized Riemann solversAnd

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;

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.

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:

wherein gamma is⁰A 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 Γ⁰Distance of (d), omega₁、Ω₂Respectively, 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:

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.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:

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 region_M+1,j、U_M+2,j、U_M+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,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:

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:

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:

where i is 0,1,2, …, M, j is 1,2, 3. U shape_i,N+1、U_i,N+2、U_i,N+3Respectively representing the conservative values, U, of the first, second and third layers of the upper boundary continuation of the calculation region_i,-1、U_i,-2、U_i,-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,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:

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:

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:

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. of_i,jAnd v_i,jRespectively representing the moving speeds of the fluids in the computing units (i, j) along the x direction and the y direction; c. C_i,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:

wherein D represents the velocity of movement of the detonation or deflagration interface, N_xAnd N_yRepresenting 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 direction_xAnd phi_yIt needs to be obtained according to the modified Godunov method. The spatial derivative phi can be determined by the following formula_xAnd phi_y：

WhereinS (phi) represents the sign function of the level set function phi.

Wherein

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:

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,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:

wherein

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:

wherein

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 equationObtaining; τ denotes an artificial time step, defined as: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 formatAnd

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:

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:

where ρ is_rAnd u_rRespectively representing the density and the speed of the unburned gas before the detonation wave;andrespectively, the density and velocity of the detonation wave post-burned gas.

For unburned gas, the interface term is calculated by:

where ρ is₁、u₁、v₁、p₁、E₁Respectively 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;represents the length of the interface truncation in the computing unit (i, j),the following equation is used to obtain:

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:

where ρ is₀、u₀、p₀、E₀Respectively 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. of₁、v₁、u_rRespectively 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:

wherein the geometric parameter alpha of the interface_i,j、Obtaining according to the step 6; conservative exchange term X (gamma) of fluid along interface_i,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 cellAnd calculating the conservation quantities at the cell boundariesAndobtaining time derivatives of variables using HWENO reconstruction techniques and generalized Riemann solversAnd

step 11.2: calculating the average conservative value after half time stepAnd calculating the conservation quantities at the cell boundariesAndobtaining time derivatives of the variables again using HWENO reconstruction techniques and generalized Riemann solversAnd

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.

Step 12: judging the calculated time t of the nth step according to the time step delta tⁿIn relation to the termination time T, if Tⁿ+ Δ T ≦ T, then Tⁿ⁺¹＝tⁿ+ Δ t, return to step 3, update the conservation quantity U at the next momentⁿ⁺¹(ii) a Else Δt＝T-tⁿAnd 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 toI.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:

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:

where i is 1,2,3, j is 0,1,2, …, 399.

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:

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.