CN116663373B - High-precision numerical simulation method suitable for gas phase combustion and explosion - Google Patents
High-precision numerical simulation method suitable for gas phase combustion and explosion Download PDFInfo
- Publication number
- CN116663373B CN116663373B CN202310913470.6A CN202310913470A CN116663373B CN 116663373 B CN116663373 B CN 116663373B CN 202310913470 A CN202310913470 A CN 202310913470A CN 116663373 B CN116663373 B CN 116663373B
- Authority
- CN
- China
- Prior art keywords
- equation
- sub
- term
- grid cell
- diffusion
- 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.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 96
- 238000002485 combustion reaction Methods 0.000 title claims abstract description 31
- 238000004088 simulation Methods 0.000 title claims abstract description 21
- 238000004880 explosion Methods 0.000 title claims abstract description 20
- 238000004364 calculation method Methods 0.000 claims abstract description 83
- 238000009792 diffusion process Methods 0.000 claims abstract description 63
- 238000006243 chemical reaction Methods 0.000 claims abstract description 40
- 239000006185 dispersion Substances 0.000 claims abstract description 9
- 230000008878 coupling Effects 0.000 claims abstract description 5
- 238000010168 coupling process Methods 0.000 claims abstract description 5
- 238000005859 coupling reaction Methods 0.000 claims abstract description 5
- 230000008569 process Effects 0.000 claims description 41
- 230000004907 flux Effects 0.000 claims description 32
- 239000000376 reactant Substances 0.000 claims description 13
- 239000000126 substance Substances 0.000 claims description 8
- 238000012545 processing Methods 0.000 claims description 7
- 239000000446 fuel Substances 0.000 claims description 6
- 230000008859 change Effects 0.000 claims description 5
- 238000010276 construction Methods 0.000 claims description 5
- 230000004913 activation Effects 0.000 claims description 3
- 230000015572 biosynthetic process Effects 0.000 claims description 3
- 238000012546 transfer Methods 0.000 claims description 3
- 125000003275 alpha amino acid group Chemical group 0.000 claims description 2
- 150000001413 amino acids Chemical class 0.000 claims description 2
- 238000011438 discrete method Methods 0.000 claims description 2
- 238000005474 detonation Methods 0.000 abstract description 24
- 230000035939 shock Effects 0.000 abstract description 6
- 239000012530 fluid Substances 0.000 abstract description 5
- 238000004458 analytical method Methods 0.000 abstract description 3
- 238000004200 deflagration Methods 0.000 abstract description 3
- 238000005516 engineering process Methods 0.000 abstract description 3
- 238000001514 detection method Methods 0.000 abstract 1
- 239000007789 gas Substances 0.000 description 13
- VNWKTOKETHGBQD-UHFFFAOYSA-N methane Chemical compound C VNWKTOKETHGBQD-UHFFFAOYSA-N 0.000 description 8
- 241000722921 Tulipa gesneriana Species 0.000 description 5
- 238000011161 development Methods 0.000 description 5
- 239000000203 mixture Substances 0.000 description 5
- NJPPVKZQTLUDBO-UHFFFAOYSA-N novaluron Chemical compound C1=C(Cl)C(OC(F)(F)C(OC(F)(F)F)F)=CC=C1NC(=O)NC(=O)C1=C(F)C=CC=C1F NJPPVKZQTLUDBO-UHFFFAOYSA-N 0.000 description 5
- 238000010586 diagram Methods 0.000 description 3
- 238000002474 experimental method Methods 0.000 description 3
- 229910052739 hydrogen Inorganic materials 0.000 description 3
- 239000001257 hydrogen Substances 0.000 description 3
- 238000013473 artificial intelligence Methods 0.000 description 2
- 238000004422 calculation algorithm Methods 0.000 description 2
- 150000002431 hydrogen Chemical class 0.000 description 2
- 238000009413 insulation Methods 0.000 description 2
- 239000003345 natural gas Substances 0.000 description 2
- 230000002093 peripheral effect Effects 0.000 description 2
- UFHFLCQGNIYNRP-UHFFFAOYSA-N Hydrogen Chemical compound [H][H] UFHFLCQGNIYNRP-UHFFFAOYSA-N 0.000 description 1
- 230000001133 acceleration Effects 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 229910052799 carbon Inorganic materials 0.000 description 1
- 238000010835 comparative analysis Methods 0.000 description 1
- 230000001066 destructive effect Effects 0.000 description 1
- 230000004069 differentiation Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 239000002737 fuel gas Substances 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 230000007774 longterm Effects 0.000 description 1
- 230000010355 oscillation Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000002269 spontaneous effect Effects 0.000 description 1
- 230000001052 transient effect Effects 0.000 description 1
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
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
- G06F17/12—Simultaneous equations, e.g. systems of linear equations
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
- G16C20/00—Chemoinformatics, i.e. ICT specially adapted for the handling of physicochemical or structural data of chemical particles, elements, compounds or mixtures
- G16C20/10—Analysis or design of chemical reactions, syntheses or processes
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
-
- 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/08—Thermal analysis or thermal optimisation
-
- 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/12—Timing analysis or timing optimisation
-
- 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
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Physics (AREA)
- Pure & Applied Mathematics (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Data Mining & Analysis (AREA)
- Chemical & Material Sciences (AREA)
- General Engineering & Computer Science (AREA)
- Databases & Information Systems (AREA)
- Analytical Chemistry (AREA)
- Algebra (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Software Systems (AREA)
- Geometry (AREA)
- Operations Research (AREA)
- Chemical Kinetics & Catalysis (AREA)
- Crystallography & Structural Chemistry (AREA)
- Life Sciences & Earth Sciences (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Computing Systems (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
The invention provides a high-precision numerical simulation method suitable for gas phase combustion and explosion. Firstly, establishing a geometric model and dividing calculation grid units; dividing a control equation set into three sets of sub-equations comprising a convection term, a diffusion term and a source term by using an operator dividing method, performing space dispersion and time propulsion on each set of sub-equations respectively, and then coupling; describing a chemical reaction process by adopting a chemical-diffusion simplified reaction model; by adopting a high-order space differential format and combining with an optimized intermittent detection and capture technology, the method remarkably reduces the calculated amount while completing high-precision reconstruction of the cell interface value. The method can capture the flow characteristics of shock waves, detonation waves and the like with strong intermittent characteristics in a flow field formed by combustion, can be used for computational fluid dynamics simulation and theoretical analysis of phenomena such as flame, deflagration, detonation, deflagration-to-detonation and the like, and has the characteristics of high accuracy, wide application range, high computational efficiency and the like.
Description
Technical Field
The invention discloses a high-precision numerical simulation method suitable for gas-phase combustion and explosion, and belongs to the field of numerical solution of gas-phase combustion processes.
Background
In long-term development, the advantages of the carbon-free low-carbon clean fuel gas energy sources such as hydrogen, natural gas and the like are considered to occupy a place in a new energy system in the future. However, hydrogen, natural gas (methane), etc. all belong to combustible gases, and once mixed with a widely existing combustion supporting substance such as air and ignited, a deflagration process of continuous propagation in the form of combustion waves may occur. In particular, in the limited space of the pipeline, the storage tank and the like, the temperature and the pressure of the pipeline and the storage tank are both increased sharply, and the pipeline and the storage tank are very likely to be converted into detonation. Detonation is a supersonic combustion wave, the propagation speed and overpressure of the detonation can reach 2000 m/s and 20 bar respectively, the destructive power is extremely strong, and extremely bad consequences such as damage to surrounding facilities, house buildings, casualties and the like can be caused. In addition, other flammable gas explosion accidents are frequent, and huge life and property losses are usually caused. Therefore, studies on combustion and explosion of gaseous combustibles such as hydrogen and methane are necessary.
In recent years, with the development of computing technology and the improvement of computing power, research on these flow and combustion problems by using numerical simulation methods such as Computational Fluid Dynamics (CFD) has been increasingly paid attention to. Compared with direct experiments, the numerical simulation method has the advantages of short period, less investment and high safety, and can simulate working conditions which cannot be realized by a plurality of experimental techniques. Numerical calculation, experiment and theoretical analysis are matched with each other, so that the method is a good way for solving the problem of flow combustion. The numerical calculation for gas phase combustion and explosion is not separated from software and programs, and at present, some commercial fluid mechanics numerical simulation software (CFD software for short) such as Fluent, FDS, FLACS and the like is internationally available, so that the universality of the software is good, and various engineering flow problems can be solved. However, such CFD software generally adopts a numerical method with lower precision, generally adopts a first-order differential format, and at most does not exceed second-order precision, and adopts turbulence mode theory on turbulence treatment. Because the software is limited by numerical precision, the software can only be used for calculating the overall characteristic or average characteristic of the flow, but the complex structure and flow details of the flow are difficult to describe, the flow field formed by the flow and combustion is difficult to accurately predict, and the accurate result is difficult to be given to the strong intermittent characteristic formed in the flow field by the supersonic combustion phenomenon such as detonation.
Disclosure of Invention
The invention provides a high-precision numerical simulation method suitable for gas phase combustion and explosion, which aims to solve the problems of low precision, large uncertainty and the like of the numerical simulation result of the existing gas phase combustion and explosion process. The current method utilizes a high-order space differential format to complete high-precision reconstruction of the grid cell interface; and solving a chemical reaction source term by adopting a chemical-diffusion simplified chemical reaction model optimized based on an artificial intelligence technology, so as to realize the reproduction of flame and detonation target characteristics such as laminar flow combustion speed, laminar flow flame thickness, adiabatic flame temperature and the like. Aiming at the problems of large calculated amount and low efficiency of high-order reconstruction, a method for combining a smoothness detector and a total variation detector is provided to improve the reconstruction efficiency and reduce the calculated amount. The method is suitable for calculating most gas phase combustion explosion scenes, can be used for calculating fluid dynamics simulation and theoretical analysis of phenomena such as deflagration, detonation, deflagration-to-detonation and the like, and has the characteristics of high accuracy, wide application range, high calculation efficiency and the like.
In order to achieve the technical purpose, the technical scheme adopted by the invention comprises the following steps:
step 1, constructing a computational domain model and dividing grid cells according to the requirements of a computational scene to obtain grid cell node information including grid cell center coordinates, grid cell interface position coordinates and grid cell sizes;
Step 2, initializing physical variables including temperature, pressure, speed and reactant components in a calculation domain according to calculation requirements, establishing a Navier-Stokes control equation set of a coupling chemistry-diffusion simplified reaction model, and setting boundary conditions according to actual physical problems;
step 3, calculating a time step according to the set Brownian number according to the requirement of explicitly solving the stability condition;
step 4, splitting the control equation set into three sets of sub-equations comprising a convection term, a diffusion term and a source term by utilizing an operator splitting method, performing space dispersion and time propulsion on each set of sub-equations respectively, and finally decoupling the three sets of sub-equations to obtain a numerical solution of the control equation set;
step 5, differentiating the flow term sub-equation separated from the Navie-Stokes control equation set in the step 4 to obtain an algebraic equation set comprising a continuity equation, a momentum equation, an energy equation and a component equation;
step 6, utilizing a five-order precision weighted essential type concussion-free differential format, and reconstructing according to the average value of the conservation variable in the grid unit in the calculation domain to obtain the conservation variable values positioned at two sides of the grid unit interface;
step 7, substituting the conservation variable values at two sides of the grid cell interface obtained by reconstruction in the step 6 into an HLLC type Riemann solver, and calculating to obtain a convection item flux value of a convection item sub-equation on the grid cell interface;
Step 8, utilizing a high-order center differential format to complete the space dispersion of a diffusion term sub-equation on the grid unit to obtain a diffusion term flux value;
step 9, calculating a chemical reaction source term based on the chemical-diffusion simplified reaction model and completing integral time propulsion of a source term sub-equation;
step 10, adopting a stable-type Dragon-Kutta format, and respectively using the convection term flux values and diffusion term fluxes of the convection term sub-equations and the diffusion term sub-equations in the steps 7 and 8 to complete time-pushing updating of the conservation variable values so as to obtainPhysical state of each substance at the moment;
step 11, judging the current calculation time t and the target calculation timeAnd the current iteration number n and the target iteration number +.>Is a relationship of (2); if->And->Updating diffusion coefficients in the chemical-diffusion simplified reaction model according to the calculation result of the current calculation moment, updating grid cell node information according to preset boundary conditions, and returning to the step 3 to start a new cycle; if->Or->And converting the conservation variable values used in each grid cell in the calculation process into original variable values, outputting the original variable values to a result file, and ending the current calculation.
Further, the step 2 includes the following:
The set of Navie-Stokes control equations is as follows:
（1）
（2）
wherein U represents a conservation variable,，/>the convection term flux in the x and y axis directions is specifically expressed as:
（3）
wherein V represents a diffusion term, consisting of the first derivative of the diffusion term flux:
（4）
wherein,,/>the diffusion term flux in the x and y axis directions is specifically expressed as:
（5）
wherein,representing viscous stresses, the viscous stresses generated in different directions are specifically expressed as:
（6）
（7）
（8）
the chemical source term S represents the rate of formation and consumption of a substance during a chemical reaction and can be expressed as:
（9）
wherein t represents time, the x-axis represents horizontal direction, the y-axis represents vertical direction,represents density, u represents velocity in x-axis direction, v represents velocity in y-axis direction, p represents pressure, +.>Represents viscous stress, E represents specific energy, K represents thermal conductivity, T represents reactant temperature, q represents chemical reaction energy release amount,/and/or the like>Represents the chemical reaction rate, Y represents the mass fraction of the reactant, R represents the universal gas constant, M represents the molecular mass,/I>Represents dynamic viscosity, D represents component diffusion coefficient, < ->Represents the specific heat capacity ratio, which in the present calculation is +.>Assuming a constant;
describing the thermodynamic state change process of the gas-phase fuel by adopting an ideal gas state equation:
（10）
The specific total energy E is expressed as:
（11）
the chemical reaction process of fuel is described by a simplified chemical-diffusion reaction model, from which the chemical reaction rate is calculated：
（12）
Wherein,indicating the activation energy of the reactants, a indicating the pre-finger factor, e indicating the natural constant, the viscosity, the diffusion of the components and the heat transfer process being affected by temperature and density and the corresponding diffusion coefficient being described by an exponential relationship:
（13）
wherein,,/>,/>indicating the viscosity, component diffusion and heat conduction, respectively, with respect to the reference coefficients,/->Representing the specific heat capacity at constant pressure.
The current calculation method adopts a cell-centered finite volume method, and adopts a virtual grid method for auxiliary processing in the implementation process of boundary conditions. According to the method, the virtual grid cells used for auxiliary calculation are arranged outside the physical area, so that the grid cells near the boundary can be processed by adopting the same space discrete scheme as the internal grid cells, independent processing of each boundary position is not needed, vectorization and parallel processing in the calculation process are facilitated, and the calculation speed and efficiency are improved.
Further, the step 4 includes the following:
because of the obvious time scale difference between the various physical processes described by the set of the Navier-Stokes equations established in step 1, it is difficult to directly couple the processes together to solve by using a conservation format, so that an operator splitting method is required to split in space and time so as to realize a step-by-step propulsion solution. By the method, a proper numerical algorithm and a time advancing format can be adopted for solving each term. The current control equation set can be divided into a convection term sub-equation, a diffusion term sub-equation and a source term sub-equation according to different physical characteristics.
（14）
（15）
（16）
The operator splitting method is adopted to complete the propulsion calculation of the equation:
（17）
wherein,sub-equations representing the convection term, diffusion term, source term, respectively, < ->Representing a time step.
Further, the step 5 includes the following:
performing differential differentiation treatment on the conservation type flow term sub-equation obtained by operator splitting, wherein equation (14) can be written as:
（18）
equation (18) is a vector form of an algebraic system of equations comprising a continuity equation, a momentum equation, an energy equation, and a component equation. Wherein,representing the average value of the conservation variable in the grid cell with the center node of the grid cell at coordinates (i, j). />And->Flux values representing convection items in the x and y axis directions, subscript +.>Representing the coordinates of the grid cell interface at which the flux value is located. The flux value may be calculated by substituting the spatially differential reconstructed value into an HLLC type Riemann solver. />And->The grid cell sizes in the x and y directions are shown, respectively. Through the above processing, the semi-discrete of the fluid term sub-equation is completed, and the space derivative is converted into an algebraic equation set.
Further, the step 6 includes the following:
in order to acquire the conservation variable values at two sides of the grid cell interface, the space distribution of the conservation variable on the grid cell is firstly acquired, and the variable at two sides of the interface is further reconstructed in a mode of step 6.1. In the reconstruction process, a mode of step 6.2 is adopted to judge whether the reconstruction standard using the WENO5 format is met, so that a space discrete method is reasonably selected, and the high-precision reconstruction is realized and the calculated amount is reduced.
Step 6.1: adopting a WENO5 format, and reconstructing the average value of the conservation variable on the center of the grid unit to obtain the conservation variable values on two sides of the interface of the grid unit;
step 6.2: selecting a proper space discrete format in the calculation process by adding a smoothness detector and controlling a total variation detector;
the following is a specific implementation manner of the two steps:
step 6.1: reconstruction of WENO5 format to obtain conservation variable values at two sides of grid cell interface
For a grid cell with a center node at (i, j), the average value of the conservation variables in the grid cell isThe coordinate +.>Conservation variable values at left and right sides of grid cell interface/>And->. Likewise, the coordinate +.A. is obtained by means of the average reconstruction of the conservation variable at the center nodes on both sides of the grid cell in the y-axis direction>Values of upper and lower sides of the grid cell interface +.>And->The grid cell interface and related variable distribution are referred to as shown in fig. 2.
Furthermore, the WENO5 format is adopted, and the average value of the conservation variable at the center of the grid cell is utilized to reconstruct the conservation variable values at the two sides of the grid cell interface. The reconstruction process of the WENO5 format in the X-axis direction is as follows. For the grid cell with the center node at (i, j), taking two nodes on the left and right sides as base frame points and building three sub-templates as shown in fig. 3:
（19）
The conservation variable values at the grid cell interface side are written as weighted averages of the reconstructed values over three sub-templates:
（20）
here the number of the elements is the number,respectively represent by sub-templates->, />, />Reconstructing the obtained conservation variable value. The differential given by the conservation variable value is required to have third-order approximation precision, and the expression is given by using Taylor series expansion of the function at the point of coordinates (i, j):
（21）
weighting factorThe definition is given below with respect to the definition,
（22）
wherein,is a small amount to ensure that the distribution is not 0, here the value is +.>. p is an integer, here 2 is less for the purpose of calculating the dissipation. />The ideal weighting factors of the smooth area obtained according to Taylor series expansion are respectively taken as follows:
（23）
this value is small in the smooth region and relatively large in the discontinuous region, which is the smoothness metric factor of the template. The calculation method of the smoothness factor is defined as:
（24）
can be solved by using symmetryThus, the reconstruction of the conservation variable values at the two sides of the interface is completed, and the method is used for solving the flux value on the interface by substituting the value into the HLLC type Riemann solver.
Step 6.2: construction smoothness detector and total variation detector
The WENO5 format adopted in the step 6.1 reduces the influence of the region where the discontinuity is located on the reconstruction polynomial through the sub-template weighted combination, thereby ensuring that stable calculation can be performed when the discontinuity occurs in the flow field. However, in the smooth flow field region, the weights of the candidate templates are similar, and the continued adoption of the WENO5 format causes unnecessary calculation consumption. Accordingly, the present invention overcomes the above-described problems by incorporating a smoothness detector as well as a total variation detector. The area containing obvious interruption in the flow field adopts a WENO5 format, and the relatively smooth area in the flow field only adopts a linear weighted intrinsic concussion-free differential format (Linear Weighted Essentially Non-oscilation), hereinafter abbreviated as LWENO format, so that the calculation time consumption is greatly reduced. The smoothness detector is constructed as follows:
（25）
In the formula (25), the amino acid sequence of the amino acid,the calculation method of (2) can be referred to as formula (24).
The total variation detector is constructed as follows:
（26）
in the formula (26), the total variationThe definition is as follows:
（27）
where k represents the number of sub-templates in the WENO5 format, r represents the number of pedestal points in each template, and l represents the number of pedestal points in the sub-template.
When (when)Or->And->The region can be considered smooth at this point, < >>，/>. Grid cells located in the smooth area need only be reconstructed using a less computationally intensive LWENO format.
In addition, the WENO5 format may generate strong minor oscillations at strong discontinuities to cause subsequent computation to diverge, and in order to overcome this problem, a policy adjustment limiter needs to be added to the WENO5 format, and the construction method thereof may be expressed as:
（28）
wherein,representing the guard tone limiter, the guard tone limiter construct used in the present invention can be represented as:
（29）
wherein,the value of the conservation variable obtained by reconstruction of the WENO5 format is shown.
When the grid unit calculates the errorAnd when the formula (30) is satisfied, the policy adjustment limiter in the formula (28) is adopted to process the WENO5 format reconstruction result, so that the calculation stability is ensured.
（30）
Further, the step 9 includes the following:
Because an operator splitting method is adopted, the source term sub-equation is solved independently, and the source term sub-equation of the components and the energy to be solved is as follows:
（31）
（32）
solving the next moment by adopting an integral time propulsion formatThereby completing the updating of the conservation variable ρY in the component equation,
（33）
wherein the superscript indicates parameters at different moments, and C indicates a constant.
Similarly, the update of the energy equation is as follows:
（34）
（35）
the density at the same time t on two sides can finish updating the conservation variable ρE:
（36）
further, the step 10 includes the following:
the time variable adopts a stable-type Dragon lattice-Kutta format to convert the semi-discrete finite volume format into a space-time full-discrete finite volume format so as to finish time propulsion:
（37）
in the above formula, L represents residual error, subscript t represents current time and subscriptThe subscripts (1), (2), (3) in the other brackets represent the number of steps advanced in the stability-preserving Longgy-Kutta format, indicating the next moment. The method can calculate the conservation variable value of the next moment by utilizing the conservation variable values of all the nodes at the current moment, and has better stability, and the maximum number of the kura can be 2 during calculation, thereby increasing the time steps in the calculation process and accelerating the calculation process.
Compared with the prior art, the invention has the advantages that:
1. the high-precision numerical format is adopted, so that the problems of low precision, large uncertainty and the like of a traditional numerical simulation result aiming at a gas-phase combustion explosion process can be solved, the details of a flow field can be accurately analyzed, the flow field structures with strong intermittent characteristics such as shock waves, detonation waves and the like can be captured, and the accuracy of a calculation result can be ensured;
2. coupling a chemical-diffusion chemical reaction model optimized based on an artificial intelligence method in a Navier-Stokes control equation set, and greatly reducing the calculated amount while realizing accurate reproduction of flame and detonation target characteristics by using a numerical calculation method;
3. the numerical format adopted by the spatial dispersion is reasonably selected by adopting a mode of combining the smoothness detector and the total variation detector, so that the efficiency of reconstruction calculation is improved.
Drawings
FIG. 1 is a flow chart of a high-precision numerical simulation method suitable for gas phase combustion and explosion of the present invention;
FIG. 2 is a schematic diagram of a reconstruction variable structure at a two-dimensional grid cell number and grid cell interface;
FIG. 3 is a schematic diagram of the pedestal points and sub-templates employed in the construction of the WENO5 process in the x-axis direction;
FIG. 4 is a graph of a two-dimensional Riemann problem density contour calculated using the scheme of the present invention in accordance with an embodiment I, wherein the arrowed lines represent streamlines;
FIG. 5 is a temperature, composition and reaction rate profile for a one-dimensional premixed laminar flame structure calculated using the inventive approach in example two;
FIG. 6 shows the results of a two-dimensional deformed tulip flame face calculated using the inventive protocol in example three; FIG. 6 (a) is a schematic diagram of the calculation domain, and FIG. 6 (b) is a two-dimensional flame surface change with time;
FIG. 7 is a deflagration-to-detonation process calculated using aspects of the present invention in example four; fig. 7 (a) is a schematic view of a calculation domain, and fig. 7 (b) is a schlieren (left) and a temperature field (right) change with time during the development of deflagration-detonation.
Detailed Description
The present invention will be described more fully hereinafter with reference to the following examples, in order to more fully understand the technical aspects of the present invention, but the scope of the present invention is not limited to the following specific examples.
As shown in fig. 1, the high-precision numerical simulation method suitable for gas phase combustion and explosion of the present invention specifically comprises the following steps:
step 1, constructing a computational domain model and dividing grid cells according to the requirements of a computational scene to obtain grid node information including grid cell center coordinates, grid cell interface positions, grid cell sizes and the like;
Step 2, initializing physical parameters such as temperature, pressure, speed, reactant components and the like in a calculation domain according to calculation requirements, establishing a multidimensional, transient and compressible Navier-Stokes control equation set of a coupling chemistry-diffusion simplified reaction model, and setting boundary conditions according to actual physical problems;
the set of established control equations is as follows:
（1）
（2）
wherein U represents a conservation variable value,，/>the convection term flux in the x and y axis directions is specifically expressed as:
（3）
wherein V represents a diffusion term, consisting of the first derivative of the diffusion term flux:
（4）
wherein,,/>the diffusion term flux in the x and y axis directions is specifically expressed as:
（5）
wherein,representing viscous stresses, the viscous stresses generated in different directions are specifically expressed as:
（6）
（7）
（8）
the chemical reaction source term S represents the rate of formation and consumption of a substance during the reaction, and can be expressed as:
（9）
wherein t represents time, x represents horizontal direction, y represents vertical direction,represents density, u represents velocity in x-axis direction, v represents velocity in y-axis direction, p represents pressure, +.>Represents the viscous stress, E represents the specific energy, K represents the heat conductivity, T represents the temperature, q represents the chemical reaction energy release amount, +.>Represents the chemical reaction rate, Y represents the mass fraction of the reactant, R represents the universal gas constant, M represents the molecular mass,/I >Representing movementsForce viscosity, D denotes the component diffusion coefficient, +.>Represents the specific heat capacity ratio, which in the present calculation is +.>It is assumed that it is a constant and,
describing the thermodynamic state change process of the gas-phase fuel by adopting an ideal gas state equation:
（10）
the specific total energy E is expressed as:
（11）
the reaction process of the fuel is described by a chemical-diffusion simplified reaction model, from which the chemical reaction rate is calculated：
（12）
Wherein,indicating the activation energy of the reactants, a indicating the pre-finger factor, e indicating the natural constant, the viscosity, the diffusion of the components and the heat transfer process being affected by temperature and density and the corresponding diffusion coefficient being described by an exponential relationship:
（13）
wherein,d, K represents viscosity, component diffusion and heat transferThe diffusion coefficient of the conduction process affected by temperature and density,,/>,/>respectively represent corresponding reference coefficients->Representing the specific heat capacity at constant pressure.
Step 3, calculating a time step according to the set Brownian according to the requirement of explicitly solving the stability condition;
step 4, splitting the control equation set into three sub-equations comprising a convection term, a diffusion term and a source term by utilizing an operator splitting method, performing space dispersion and time propulsion on each set of equations respectively, and finally combining the solutions of the three sub-equations to obtain a numerical solution of the control equation;
Because of the obvious time scale difference between the various physical processes described by the set of the Navier-Stokes equations established in step 2, it is difficult to directly couple the processes together to solve by using a conservation format, so that an operator splitting method is required to split in space and time so as to realize a step-by-step propulsion solution. By the method, a proper numerical algorithm and a time advancing format can be adopted for solving each term. The current control equation set can be divided into a convection term sub-equation, a diffusion term sub-equation and a source term sub-equation according to different physical characteristics.
（14）
（15）
（16）
The propulsion calculation of the equation is completed by adopting a second order operator splitting method:
（17）
wherein,integrating operators respectively representing a stream item, a diffusion item and a source item, wherein the time step is that。
Step 5, differentiating the flow term sub-equation separated from the conservation type Navier-Stokes control equation set in the step 4 to obtain an algebraic equation set comprising a continuity equation, a momentum equation, an energy equation and a component equation;
（18）
wherein,representing the value of the conservation variable stored in the grid cell with the central node at coordinate (i, j). />And->Flux values representing convection items in the x-axis direction and the y-axis direction, subscript +. >Representing the coordinates of the grid cell interface at which the flux value is located. The flux values described above may be obtained by substituting the differentially reconstructed values into HLLC type Riemann solver calculations. />And->The grid dimensions in the x and y directions are shown, respectively. Through the above processing, the semi-discrete of the conservation type convection term sub-equation is completed, and the space derivative is converted into an algebraic equation set.
And 6, reconstructing the conservation variable mean value of the grid cells in the two-dimensional calculation domain by using a space high-order difference method to obtain the conservation variable values at two sides of the grid cell interface, wherein the related grid cell interface and related variable distribution conditions are shown in figure 2. The detector is added to realize the selection of a linear weighted intrinsic concussion-free differential format (Linear Weighted Essentially Non-oscilation) in a smooth area, hereinafter abbreviated as LWENO format, and the calculated amount is controlled; a five-order precision weighted intrinsic concussion-free differential format (Weighted Essentially Non-oscilation) is adopted in a strong discontinuous region where shock waves, detonation waves and the like exist, and is hereinafter abbreviated as WENO5 format, so that the calculation precision is improved;
first, we no5 space-value discrete format is used to reconstruct the values of the conservation variables at both sides of the grid cell interface based on the average value of the conservation variables at the center of the grid cell. Here, we no5 format reconstruction in the x-axis direction is described as an example. As shown in fig. 3, for a grid with a central node at coordinates (i, j), two nodes on each of the left and right sides are taken as base frame points and three sub-templates are created:
（19）
The conservation variable values on the grid cell interface side can then be written as a weighted average of the reconstructed values over three templates:
（20）
here the number of the elements is the number,respectively represent by sub-templates->,/>,/>And calculating the conservation variable value at one side of the obtained interface. The difference given by the conservation variable values is required to have third-order approximation precision, and the expression is easily given by using Taylor series expansion of the function at the point of coordinates (i, j):
（21）
selecting reasonable weighting factorsIs critical to the WENO5 format, where the weighting factors are defined as follows:
（22）
wherein,is a small amount to ensure that the distribution is not 0, here the value is +.>. p is an integer, here 2 is less for the purpose of calculating the dissipation. />The ideal weighting factors of the smooth area obtained according to Taylor series expansion are respectively taken as follows:
（23）
this value is small in the smooth region and relatively large in the discontinuous region, which is the smoothness metric factor of the template. The calculation method of the smoothness factor is defined as:
（24）
can be obtained by using symmetryAnd the reconstruction of the conservation variable values at the two sides of the grid cell interface is completed and is used for solving the flux value on the interface by substituting the conservation variable values into the HLLC type Riemann solver.
And selecting a proper space discrete format in the calculation process by defining a smoothness detector and controlling a total variation detector, wherein the smoothness detector is constructed in the following way:
（25）
In the expression (25), the calculation method of IS can refer to the expression (24).
The total variation detector is constructed as follows:
（26）
in the formula (26), the total variationThe definition is as follows:
（27）
where k represents the sub-template number in the WENO5 format, r represents the number of pedestal points in each template, and l represents the current pedestal point number.
When (when)Or->And->The region can be considered smooth at this point, < >>，/>. Grid cells located in the smooth area need only be reconstructed using a less computationally intensive LWENO format.
In addition, the WENO5 format may generate strong micro-oscillation at a strong break, which may cause subsequent calculation divergence, and in order to overcome this problem, a policy adjustment limiter needs to be added on the basis of the WENO5 format, and the construction method can be expressed as:
（28）
wherein,representing the guard tone limiter, the guard tone limiter construct used in the present invention can be represented as:
（29）
wherein,the value of the conservation variable obtained by reconstruction of the WENO5 format is shown.
When the grid unit calculates the errorWhen equation (30) is satisfied, the result of WENO5 reconstruction is performed using the policy adjustment limiter described in equation (29)Processing is performed, thereby ensuring computation stability.
（30）
Step 7, substituting the values at two sides of the unit obtained by reconstruction in the step 6 into an HLLC format Riemann solver, and calculating to obtain a convection item flux value of a convection item sub-equation on a grid unit interface;
Step 8, utilizing a high-order center differential format to complete the space dispersion of a diffusion term sub-equation on the grid unit to obtain a diffusion flux value;
step 9, calculating a chemical reaction source term based on the chemical-diffusion simplified reaction model and completing integral time propulsion of a source term sub-equation; because an operator splitting method is adopted for solving, the source term sub-equation is independently solved, and the components to be solved and the energy sub-equation are respectively:
（31）
（32）
solving the next moment by adopting an integral time propulsion formatAnd the component concentration of the (c) is changed, so that updating of a conservation variable rho Y in a component equation is completed.
（33）
Wherein, the superscript indicates parameters at different moments, and C indicates a constant;
similarly, the update of the energy equation is as follows:
（34）
（35）
the density at the same time t on two sides can finish updating the conservation variable ρE:
（36）
step 10, adopting a stable-type Dragon-Kutta format, and respectively using the convection term flux value and diffusion term flux of the convection term and diffusion term sub-equation in the steps 7 and 8 to complete time-pushing updating of the conservation variable value so as to obtainPhysical state of each substance at the moment:
（37）
in the above formula, L represents residual error, subscript t represents current time and subscriptThe subscripts (1), (2), (3) in the other brackets represent the number of steps advanced in the stability-preserving Longgy-Kutta format, indicating the next moment. The method can calculate the conservation variable value of the next moment by utilizing the conservation variable values of all grid cells at the current moment, and has better stability, and the maximum number of the kura can be 2 during calculation, thereby increasing the time steps in the calculation process and accelerating the calculation process.
Step 11, judging the current calculation time t and the target calculation timeAnd the current iteration number n and the target iteration number +.>Is a relationship of (2); if->And->Updating diffusion coefficients in the chemical-diffusion simplified reaction model according to the calculation result of the current calculation moment, updating grid cell information according to preset boundary conditions, and returning to the step 3 to start a new cycle; if->Or->And converting the conservation variable values used in each grid cell in the calculation process into original variable values, outputting the original variable values to a result file, and ending the current calculation.
The following are 3 examples given as specific embodiments of the disclosed method.
Embodiment one-dimensional and two-dimensional Riemann problem
This problem describes the interaction between two shock waves and two contact discontinuities, with the initial conditions:
,
wherein,the density is represented by u, the velocity along the x-axis, v, the velocity along the y-axis, and p, the pressure. The calculated domain size is +.>Grid cell resolution is +.>. This example, although being a numerical example only, includes a more complex flow structure that can be used to demonstrate the higher numerical accuracy of the present invention. Fig. 4 shows the flow field density distribution calculated using the inventive scheme at t=0.25. It can be seen that the scheme of the invention has a relatively high effect The high spatial resolution can accurately capture the shock wave surface position and the fine vortex structure in the interference zone wave system.
Example two, one-dimensional premixed laminar flame problem
The problem is a one-dimensional laminar flame problem, mainly describes a laminar flame surface development process formed after ignition of a one-dimensional premixed hydrogen-air mixture, and performs comparative analysis on numerical simulation and flame surface characteristic parameters obtained by theoretical calculation so as to verify the accuracy of the result of the current calculation method. The current calculation uses a calculated domain length of 1 cm, an initial reactant temperature of 293K, an initial pressure of 1 atm, and a stoichiometric ratio of 1 for the premixed reactants. A high temperature burnt zone is provided to the left of the zone to act as an ignition source. The temperature, composition, and reaction rate profile obtained in the calculation domain after ignition is formed is shown in fig. 5. The horizontal axis in the figure indicates the lateral distance, and fig. 5 shows only the structure near the flame face and redraws the horizontal axis range. The scattered points represent the results of numerical calculations, and the curves represent the exact solution results obtained from theoretical calculations. As can be seen from fig. 5, the calculated temperatureAbout 2141K, laminar flame speed +.>About 3.12 m/s, flame face thickness +. >Flame face thickness +.0.34 0.34 mm>Is defined as the distance between the component concentrations of 0.1 to 0.9. The result obtained by the calculation method has good matching degree with the accurate solution, and the error range is about 0.1 percent.
Three-dimensional deformation Tulip flame problem
This problem is a two-dimensional laminar flame problem and mainly describes the process of premixing the deformed tulip flame surface formed after ignition of the hydrogen-air mixture in a two-dimensional planar tube, as shown in fig. 6. The calculation domain of this problem is a two-dimensional planar pipeline with d=2 cm and l=14 cm, as shown in fig. 6 (a), and the initial temperature in the pipeline is 293K, and the initial pressure is 1 atm. The stoichiometric ratio of the premix reaction was 1. A semicircular high-temperature burnt area is arranged at the leftmost side of the pipeline and is used as an ignition source. The peripheral boundaries adopt the condition of heat insulation and no sliding wall surface. The situation in which the flame front evolves with time in the duct after ignition is shown in fig. 6 (b). It can be seen from the figure that different flame surface characteristics of fingertip flame, tulip flame, deformed tulip flame and the like can be successfully calculated by adopting the method of the invention, but commercial software adopting a low-order format cannot accurately capture the characteristics, so that the difference between the subsequent calculation result and the experimental result is obvious. The method has higher calculation precision, and is more suitable for capturing the detail characteristics of the flow field in the combustion process so as to ensure the reliability of the calculation result.
Fourth embodiment, deflagration-to-detonation problem in coarse and fine pipeline
This problem is a two-dimensional turbulent flame development problem and mainly describes the flame acceleration and deflagration-to-detonation process after ignition of a premixed hydrogen-air mixture in a two-dimensional coarse fine pipe, as shown in fig. 7. The calculation domain of this problem is a two-dimensional planar pipeline with a coarse element inside d=5 mm, as shown in fig. 7 (a), and the coarse element is an equilateral triangle with a height h=500 μm continuously distributed on the wall surface. The initial temperature in the pipeline was 293K and the initial pressure was 1 atm. The stoichiometric ratio of the premix reaction was 1. And arranging a planar high-temperature burnt area at the leftmost side of the pipeline as an ignition source. The peripheral boundaries adopt the condition of heat insulation and no sliding wall surface. The process of converting detonation inside the pipeline to detonation after ignition over time is shown in fig. 7 (b). The method can be used for successfully calculating the processes of spontaneous combustion of the boundary layer, generation of local detonation and global detonation caused by the local detonation wave. The results show that the method has high calculation accuracy and can capture the flow field characteristics with intermittent characteristics such as shock waves, detonation waves and the like generated in the detonation process.
While the foregoing describes illustrative embodiments of the present invention to facilitate an understanding of the present invention by those skilled in the art, it should be understood that the present invention is not limited to the scope of the embodiments, but is to be construed as protected by the accompanying claims insofar as various changes are within the spirit and scope of the present invention as defined and defined by the appended claims.
Claims (7)
1. A high-precision numerical simulation method suitable for gas phase combustion and explosion is characterized in that: the method comprises the following steps:
step 1, constructing a computational domain model and dividing grid cells according to the requirements of a computational scene to obtain grid cell node information including grid cell center coordinates, grid cell interface position coordinates and grid cell sizes;
step 2, initializing physical variables including temperature, pressure, speed and reactant components in a calculation domain according to calculation requirements, establishing a Navier-Stokes control equation set of a coupling chemistry-diffusion simplified reaction model, and setting boundary conditions according to actual physical problems;
step 3, calculating a time step according to the set Brownian number according to the requirement of explicitly solving the stability condition;
step 4, splitting the control equation set into three sets of sub-equations comprising a convection term, a diffusion term and a source term by utilizing an operator splitting method, performing space dispersion and time propulsion on each set of sub-equations respectively, and finally decoupling the three sets of sub-equations to obtain a numerical solution of the control equation set;
step 5, differentiating the flow term sub-equation separated from the Navie-Stokes control equation set in the step 4 to obtain an algebraic equation set comprising a continuity equation, a momentum equation, an energy equation and a component equation;
Step 6, reconstructing the conservation variable values positioned at two sides of the grid cell interface according to the average value of the conservation variable in the grid cell in the calculation domain by using a five-order precision weighted essential type concussion-free differential format;
step 7, substituting the conservation variable values at two sides of the grid cell interface obtained by reconstruction in the step 6 into an HLLC type Riemann solver, and calculating to obtain a convection item flux value of a convection item sub-equation on the grid cell interface;
step 8, completing the space dispersion of the sub-equation of the diffusion term on the grid unit by using a high-order center difference format and calculating the flux value of the diffusion term;
step 9, calculating a chemical reaction source term based on the chemical-diffusion simplified reaction model and completing integral time propulsion of a source term sub-equation;
step 10, adopting a stable-type Dragon-Kutta format, and respectively using the convection term sub-equation, the convection term flux value and the diffusion term flux of the diffusion term sub-equation in the steps 7 and 8 to complete the time pushing update of the conservation variable value so as to obtain the physical state of each substance at the moment of t+delta t;
step 11, judging the current calculation time t and the target calculation timeCurrent iteration number n and target iteration numberIs a relationship of (2); if->And->Updating diffusion coefficients in the chemical-diffusion simplified reaction model according to the calculation result of the current calculation moment, updating grid cell node information according to preset boundary conditions, and returning to the step 3 to start a new cycle; if- >Or->And converting the conservation variable values used in each grid cell in the calculation process into original variable values, outputting the original variable values to a result file, and ending the current calculation.
2. The high-precision numerical simulation method suitable for gas-phase combustion and explosion according to claim 1, wherein: the step 2 comprises the following steps:
the set of Navie-Stokes control equations is as follows:
（1）
（2）
wherein U represents a conservation variable,，/>the convection term flux in the x and y axis directions is specifically expressed as:
（3）
wherein V represents a diffusion term, consisting of the first derivative of the diffusion term flux:
（4）
wherein,,/>the diffusion term flux in the x and y axis directions is specifically expressed as:
（5）
wherein,representing viscous stresses, the viscous stresses generated in different directions are specifically expressed as:
（6）
（7）
（8）
the source term S represents the rate of formation and consumption of a substance during a chemical reaction, expressed as:
（9）
wherein t represents time, the x-axis represents horizontal direction, the y-axis represents vertical direction,represents density, u represents velocity in x-axis direction, v represents velocity in y-axis direction, p represents pressure, +.>Represents viscous stress, E represents specific energy, K represents thermal conductivity, T represents reactant temperature, q represents chemical reaction energy release amount,/and/or the like >Represents the chemical reaction rate, Y represents the mass fraction of the reactant, R represents the universal gas constant, M represents the molecular mass,/I>Represents dynamic viscosity, D represents component diffusion coefficient, < ->Represents the specific heat capacity ratio, which in the present calculation is +.>Assuming a constant;
describing the thermodynamic state change process of the gas-phase fuel by adopting an ideal gas state equation:
（10）
the specific total energy E is expressed as:
（11）
the chemical reaction process of fuel is described by a simplified chemical-diffusion reaction model, from which the chemical reaction rate is calculated：
（12）
Wherein,indicating the activation energy of the reactants, a indicating the pre-finger factor, e indicating the natural constant, the viscosity, the diffusion of the components and the heat transfer process being affected by temperature and density and the corresponding diffusion coefficient being described by an exponential relationship:
（13）
wherein,d, K represents the diffusion coefficient of viscosity, component diffusion and heat conduction process affected by temperature and density, +.>, ,/>Respectively represent corresponding reference coefficients->Representing the specific heat capacity at constant pressure.
3. A high-precision numerical simulation method suitable for gas-phase combustion and explosion according to claim 2, characterized in that: the step 4 comprises the following steps:
for the Navie-Stokes control equation set established in the step 2, an operator splitting method is adopted to split in space and time so as to realize step-by-step propulsion solution, the current Navie-Stokes control equation set is divided into a convection term sub-equation, a diffusion term sub-equation and a source term sub-equation according to different physical characteristics,
（14）
（15）
（16）
The operator splitting method is adopted to complete the propulsion calculation of the equation:
（17）
wherein,sub-equations representing the convection term, diffusion term, source term, respectively, < ->Representing a time step.
4. A high-precision numerical simulation method suitable for gas-phase combustion and explosion according to claim 3, wherein: the step 5 comprises the following steps:
differentiating the flow term sub-equation obtained by using an operator splitting method, and writing an equation (14):
（18）
equation (18) is a vector form of an algebraic system of equations comprising a continuity equation, a momentum equation, an energy equation, and a component equation, wherein,mean value of conservation variable in grid cell representing that center node of grid cell is located at coordinates (i, j),/>And->Flux values representing convection items in the x and y axis directions, subscript +.>Representing coordinates of the grid cell interface at which the flux value is located，/>And->The grid cell sizes in the x and y directions are expressed respectively, and through the above processing, the semi-dispersion of the flow term sub-equation is completed, and the spatial derivative is converted into an algebraic equation set.
5. The high-precision numerical simulation method suitable for gas-phase combustion and explosion according to claim 1, wherein: the step 6 comprises the following steps:
In order to acquire the conservation variable values at two sides of the grid cell interface, firstly, acquiring the space distribution of the conservation variable values on the grid cell, further reconstructing the conservation variable values at two sides of the grid cell interface in a mode of step 6.1, completing the reconstruction of the step by adopting a five-order precision weighted intrinsic concussion-free differential format, and judging whether the reconstruction standard by adopting the five-order precision weighted intrinsic concussion-free differential format is met or not in the reconstruction process in a mode of step 6.2, so that a space discrete method is reasonably selected, and the calculation amount is reduced while the high-precision reconstruction is realized; the method comprises the following specific steps:
step 6.1: adopting a five-order precision weighted intrinsic concussion-free differential format, and reconstructing the average value of the conservation variable of the grid cell on the central node of the grid cell to obtain the conservation variable values at two sides of the grid cell interface;
step 6.2: selecting a proper space discrete format in the calculation process by adding a smoothness detector and controlling a total variation detector; the specific implementation manner of the step 6.1 is as follows:
for a grid cell with its center node at coordinates (i, j), the average value of the conservation variables in the grid cell is Then utilize the gate keeper at the center nodes of the grid unit along the two sides of the x-axis directionAverage reconstruction of constant variable value to obtain coordinatesConservation variable value of left and right sides at grid cell interface +.>And->Likewise, the coordinate ++is obtained by using the average value reconstruction of the conservation variable at the center nodes on both sides of the grid cell in the y-axis direction>Conservation variable values of upper side and lower side of grid interface +.>And->；
Further, a five-order precision weighted intrinsic concussion-free differential format is adopted, the conservation variable values at two sides of the grid unit interface are reconstructed by using the conservation variable average value at the center node of the grid unit, and the reconstruction process of the five-order precision weighted intrinsic concussion-free differential format in the x-axis direction is as follows:
for the grid unit with the central node at the coordinates (i, j), taking two nodes on the left and right sides of the grid unit as base frame points and establishing three sub-templates:
（19）
the conservation variable values at the grid cell interface side are written as weighted averages of the reconstructed values over three sub-templates:
（20）
here the number of the elements is the number,respectively represent by sub-templates->, />,/>Reconstructing the obtained conservation variable value, wherein the difference given by the conservation variable value is required to have third-order approximation precision, and the expression is given by using Taylor series expansion of a function at a coordinate (i, j) point:
（21）
Weighting factorThe definition is given below with respect to the definition,
（22）
wherein,is a small amount to ensure that the distribution is not 0, here the value is +.>P is an integer, where 2, ++for the purpose of calculating the dissipation is smaller>The ideal weighting factors of the smooth area obtained according to Taylor series expansion are respectively taken as follows:
（23）
as the smoothness factor of the template, the calculation mode of the smoothness factor is defined as:
（24）
obtained by using symmetry principleThereby completing the reconstruction of the conservation variable values at the two sides of the interface, and being used for solving the flux value on the interface by substituting the HLLC type Riemann solver subsequently;
the method of constructing the smoothness detector and the total variation detector in step 6.2 is as follows,
the smoothness detector is constructed as follows:
（25）
in the formula (25), the amino acid sequence of the amino acid,with reference to the calculation method of (2),
the total variation detector is constructed as follows:
（26）
in the formula (26), the total variationThe definition is as follows:
（27）
wherein k represents the number of sub-templates in the five-order precision weighted essential type concussion-free differential format, r represents the number of base frame points in each sub-template, and l represents the number of base frame points in the sub-template;
when (when)Or->And->The area is considered smooth, where,，/>the grid unit positioned in the smooth area adopts a linear weighted intrinsic concussion-free differential format with smaller calculated amount during reconstruction;
A policy adjustment limiter is added on the basis of a five-order precision weighted intrinsic concussion-free differential format, and the construction method is as follows:
（28）
wherein,representing a policy adjustment limiter, the configuration of which is represented as:
（29）
wherein,representing conservation variable values obtained by reconstruction of a five-order precision weighted intrinsic concussion-free differential format;
when the grid unit calculates the errorWhen the formula (30) is satisfied, the policy adjustment limiter in the formula (28) is adopted to process the five-order precision weighted intrinsic concussion-free differential format result, so as to ensure stable calculation,
（30）。
6. a high-precision numerical simulation method suitable for gas-phase combustion and explosion according to claim 2, characterized in that: the step 9 comprises the following steps:
because an operator splitting method is adopted, the source term sub-equation is solved independently, and the source term sub-equation of the components and the energy to be solved is as follows:
（31）
（32）
solving the next moment by adopting an integral time propulsion formatThe component concentration of the component is changed, thereby completing the updating of the conservation variable ρY in the source item sub-equation of the component,
（33）
wherein, the superscript t is marked up,parameters representing different moments, C representing a constant;
similarly, the update of the source sub-equation for energy is:
（34）
（35）
The density at the same time t on two sides completes updating the conservation variable ρE:
（36）。
7. the high-precision numerical simulation method suitable for gas-phase combustion and explosion according to claim 1, wherein: the step 10 includes the following:
the time variable adopts a stable-type Dragon lattice-Kutta format to convert the semi-discrete finite volume format into a space-time full-discrete finite volume format so as to finish time propulsion:
（37）
in the above formula, L represents residual error, subscript t represents current time and subscriptThe next time, the other subscripts (1), (2), (3) in brackets represent the number of time-pushed steps in the stability-preserving Longgy-Kutta format.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310913470.6A CN116663373B (en) | 2023-07-25 | 2023-07-25 | High-precision numerical simulation method suitable for gas phase combustion and explosion |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310913470.6A CN116663373B (en) | 2023-07-25 | 2023-07-25 | High-precision numerical simulation method suitable for gas phase combustion and explosion |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116663373A CN116663373A (en) | 2023-08-29 |
| CN116663373B true CN116663373B (en) | 2023-11-17 |
Family
ID=87715600
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310913470.6A Active CN116663373B (en) | 2023-07-25 | 2023-07-25 | High-precision numerical simulation method suitable for gas phase combustion and explosion |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116663373B (en) |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| FR3031210A1 (en) * | 2014-12-31 | 2016-07-01 | Landmark Graphics Corp | SEISMIC ELASTIC WAVE SIMULATION FOR TRANSVERSALLY INCLINED ISOTROPIC ENVIRONMENT USING DECADED LEBEDEV GRID |
| EP3499161A1 (en) * | 2017-12-13 | 2019-06-19 | Novamet Sàrl | Method of controlling a gas furnace for melting metal and control system therefor |
| CN112082798A (en) * | 2020-09-14 | 2020-12-15 | 中国科学技术大学 | Visual test device for accurately testing unsteady detonation flame arrester effect of combustible gas |
| CN112163312A (en) * | 2020-08-17 | 2021-01-01 | 空气动力学国家重点实验室 | Method for carrying out numerical simulation on compressible flow problem through high-order WENO format reduction |
| CN114169184A (en) * | 2021-10-19 | 2022-03-11 | 北京理工大学 | High-precision self-adaptive finite volume-finite difference coupling numerical simulation method |
| CN115935841A (en) * | 2022-12-02 | 2023-04-07 | 山东科技大学 | Numerical simulation method for explosion propagation characteristic of multi-element gas in semi-closed space |
| CN116305789A (en) * | 2023-02-06 | 2023-06-23 | 北京理工大学 | Multiphase explosion flow field temperature simulation method based on saturated vapor pressure model |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20160004802A1 (en) * | 2014-07-03 | 2016-01-07 | Arizona Board Of Regents On Behalf Of Arizona State University | Multiscale Modelling of Growth and Deposition Processes in Fluid Flow |
-
2023
- 2023-07-25 CN CN202310913470.6A patent/CN116663373B/en active Active
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| FR3031210A1 (en) * | 2014-12-31 | 2016-07-01 | Landmark Graphics Corp | SEISMIC ELASTIC WAVE SIMULATION FOR TRANSVERSALLY INCLINED ISOTROPIC ENVIRONMENT USING DECADED LEBEDEV GRID |
| EP3499161A1 (en) * | 2017-12-13 | 2019-06-19 | Novamet Sàrl | Method of controlling a gas furnace for melting metal and control system therefor |
| CN112163312A (en) * | 2020-08-17 | 2021-01-01 | 空气动力学国家重点实验室 | Method for carrying out numerical simulation on compressible flow problem through high-order WENO format reduction |
| CN112082798A (en) * | 2020-09-14 | 2020-12-15 | 中国科学技术大学 | Visual test device for accurately testing unsteady detonation flame arrester effect of combustible gas |
| CN114169184A (en) * | 2021-10-19 | 2022-03-11 | 北京理工大学 | High-precision self-adaptive finite volume-finite difference coupling numerical simulation method |
| CN115935841A (en) * | 2022-12-02 | 2023-04-07 | 山东科技大学 | Numerical simulation method for explosion propagation characteristic of multi-element gas in semi-closed space |
| CN116305789A (en) * | 2023-02-06 | 2023-06-23 | 北京理工大学 | Multiphase explosion flow field temperature simulation method based on saturated vapor pressure model |
Non-Patent Citations (1)
| Title |
|---|
| 基于燃烧与水动力耦合模型的锅炉蒸汽管超温特性研究;喻聪，等;《热能动力工程》;第92-98页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116663373A (en) | 2023-08-29 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Ellail et al. | Description and validation of a three-dimensional procedure for combustion chamber flows | |
| Xiao et al. | Dynamics of premixed hydrogen/air flame in a closed combustion vessel | |
| Liu et al. | Numerical investigation of hydrogen detonation suppression with inert particle in pipelines | |
| Li et al. | Numerical application of additive Runge-Kutta methods on detonation interaction with pipe bends | |
| CN116663373B (en) | High-precision numerical simulation method suitable for gas phase combustion and explosion | |
| Kobiera et al. | A new phenomenological model of gas explosion based on characteristics of flame surface | |
| CN109977431B (en) | Smoke modeling method in large-scene environment | |
| Kuzuu et al. | Numerical simulation of premixed flame propagation in a closed tube | |
| Shang et al. | Unstructured adaptive grid method for reacting flow computation | |
| Nicolás-Pérez et al. | Evaluation of different models for turbulent combustion of hydrogen-air mixtures. Large Eddy Simulation of a LOVA sequence with hydrogen deflagration in ITER Vacuum Vessel | |
| Yao et al. | A pressure-stabilized Lagrange-Galerkin method in a parallel domain decomposition system | |
| Taylor et al. | Mesh-independent large-eddy simulation with anisotropic adaptive mesh refinement for hydrogen deflagration prediction in closed vessels | |
| CN113591345A (en) | Explosion reaction flow high-precision prediction method based on generalized Riemann solution method device | |
| Cercignani | The Boltzmann equation and fluid dynamics | |
| Shen et al. | The meshless direct simulation Monte Carlo method | |
| Zhang et al. | Development of a Multi-Phase Flamelet Generated Manifold for Spray Combustion Simulations | |
| Heidarinejad et al. | Simulation of premixed combustion flow around circular cylinder using hybrid random vortex | |
| Qizhong et al. | Numerical simulation study of gas explosion in confined space based on deep learning algorithm | |
| Liu et al. | Detonation propagation characteristics of H2/O2/AR in multi-angle bifurcation tubes | |
| CN114021499A (en) | Method for calculating heat conduction of aircraft thermal protection structure based on FVM-TLBFS method | |
| VANKA | Block-implicit computation of viscous internal flows-recent results | |
| Kumar et al. | Study of Leakage and Explosion of Hydrogen and Blast Wall Failures in an Offshore Platform | |
| Jeong et al. | A study on smoke movement in room fires with various pool fire location | |
| Ingram et al. | On the role of turbulence in detonation induced by Mach stem reflection | |
| Li et al. | Numerical simulations of transient flows in ram accelerators |
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 |