Calculation method for the explosion process of combustible gas based on the reaction progress variable model

Abstract

The present invention provides a calculation method for the combustion process of combustible gas based on a reaction progress variable model, comprising the following steps: coupling auto-ignition and flame propagation through the reaction progress variable according to the construction of the source phase of the reaction progress variable equation; realizing smooth switching between a pressure-based solver and a density-based solver by judging the critical turbulent Mach number, thereby improving the calculation accuracy of the solver in the calculation of the turbulent flame acceleration process of deflagration-to-detonation transition; and improving the calculation accuracy of the explosion initiation process by establishing an ignition initiation model. The present invention constructs the source term of the reaction progress variable turbulent combustion model, uses the reaction progress variable to couple the flame surface density and the auto-ignition model; and improves the solution accuracy of the flow field variables in the control unit containing shock waves by the method of volume fraction weighting.

Description

Technical Field

[0001] The present invention relates to the technical field of computational fluid dynamics, and in particular, to a computational method for the deflagration-to-detonation turbulent flame acceleration process based on a reaction progress variable model. Background Art

[0002] With the development of the economy, problems such as environmental pollution, climate anomalies, and energy shortages have increasingly become the focus of public attention. Hydrogen energy is an important energy source for addressing global climate change and building a decarbonized society. One of the current energy security issues is to ensure explosion safety when using flammable gases. The most urgent problem is the safety of hydrogen, which has a wide explosion concentration range and a low ignition energy. Under obstacle conditions, it accelerates the flame transition from deflagration to detonation (DDT), which is the most destructive mode in gas explosions. As Figure 1 shown, the flame undergoes five stages from ignition to detonation: slow deflagration stage, fast deflagration stage, pre-detonation stage, over-driven detonation stage, and steady detonation stage. For a stoichiometric hydrogen-air mixture, the steady detonation stage can produce a flame propagation speed of about 1800 - 2000 m / s, with an overpressure as high as more than 20 atmospheres or higher. Numerical simulation can study and predict the explosion process and plays a crucial role in hydrogen safety engineering. Due to the highly nonlinear nature of the NS equations describing the hydrogen explosion process, the numerical solution process incurs a huge computational cost.

[0003] The turbulent combustion flamelet model was first proposed by the famous German scholar Norbert Peters based on the study of laminar diffusion flames. Under the flamelet model assumption, the internal structure of the flame is not affected by turbulent vortices and only undergoes distortion under the action of turbulence. Under this physical mechanism, the internal structure of the flamelet and the action of turbulence can be considered separately. In this way, decoupling between flow and chemical reaction is achieved to a certain extent, and various combustion phenomena can be predicted with a relatively small computational cost. Its high computational efficiency makes it very suitable for application in engineering. The flamelet model is restricted by specific combustion conditions, such as Figure 2 According to the degree of interaction between the flame and turbulence, the combustion process can be divided into: laminar flame zone, corrugated flame zone, wrinkled flame zone, thin flame zone, and broken-vortex reaction zone. The explosion acceleration of hydrogen passes through the "corrugated flame zone" and "thin flame zone" and may even reach the "broken-vortex reaction zone" mode. Quenching may occur with the increasing intensity of the interaction between turbulence and the flamelet. At the same time, due to the action of high-energy leading detonation waves, auto-ignition may occur in the unburned mixture region. According to the modal analysis of the flame, the original reaction progress variable model cannot be applied to the numerical calculation of the explosion kinetics of hydrogen, and its source term needs to be re-modeled to consider flame quenching and auto-ignition behaviors.

[0004] The most accurate method for predicting complex auto-ignition phenomena currently is to apply detailed chemical reaction models that include hundreds of intermediate components and reactions. These detailed reaction models can describe all auto-ignition phenomena across the entire range of temperature, pressure, and mixture fraction. However, with the existing CPU computing power, even for practical simplified mechanism models, it is impossible to directly solve the transport of these chemical reactions and intermediate components in industrial applications.

[0005] Analyzing the chemical reactions within the flame front (~0.3 - 0.5 mm) requires a grid size of (~10 - 125 μm). The extremely high propagation speed (~1800 - 2000 m / s) requires even smaller time steps, consuming huge computing resources. According to the literature (Heidari, 2014), simulating the DDT process within a 20 ms pipeline typically takes 2 - 4 months (on a 64 CPU workstation).

[0006] Meanwhile, the basic reaction progress variable solver cannot simulate the auto-ignition and quenching phenomena during the explosion process. The general pressure-based solver cannot accurately capture the shock discontinuity surface (~10 nm), and the premixed combustion model cannot be used in the density-based solver in commercial software (FLUENT, CFX). Summary of the Invention

[0007] In view of the technical problems mentioned in the above background art, a calculation method for the deflagration-to-detonation turbulent flame acceleration process based on the reaction progress variable model is provided. The present invention mainly utilizes

[0008] Compared with the prior art, the present invention has the following advantages:

[0009] The present invention discloses a calculation method for the deflagration-to-detonation turbulent flame acceleration process based on the reaction progress variable model, belonging to the field of computational fluid dynamics (CFD). The present invention constructs the source term of the reaction progress variable turbulent combustion model, uses the reaction progress variable to couple the flame surface density and the auto-ignition model; and improves the solution accuracy of the flow field variables within the control unit containing shock waves through the method of volume fraction weighting. By judging the critical turbulent Mach number, a smooth switch between the pressure-based solver and the density-based solver is achieved, improving the calculation accuracy of the solver in the turbulent flame acceleration and deflagration-to-detonation processes. An ignition model is developed to achieve precise control of the explosion startup process. The present invention can accurately and quickly predict the explosion process, providing key guidance for safety planning and explosion protection engineering.

[0010] Meanwhile, the present invention can also be applied to the numerical calculation of the explosion of other flammable gases besides hydrogen. Brief Description of the Drawings

[0011] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required for the description of the embodiments or the prior art. Obviously, the drawings in the following description are some embodiments of the present invention. For those of ordinary skill in the art, without creative efforts, other drawings can also be obtained based on these drawings.

[0012] Figure 1 It is a schematic diagram of the evolution of the combustion system from deflagration to detonation in the background art of the present invention;

[0013] Figure 2 It is the Borghi diagram of the characteristic state of turbulent combustion of the present invention;

[0014] Figure 3 It is the curve of the laminar flame propagation speed, expansion coefficient and total thermodynamic index of premixed hydrogen-air varying with the hydrogen concentration of the present invention;

[0015] Figure 4 It is a schematic diagram of the flow field distribution within the shock wave-containing grid unit of the present invention;

[0016] Figure 5 It is a schematic diagram of the Sod shock tube of the present invention;

[0017] Figure 6 It is the comparison between the calculation results of the pressure-based PISO scheme and the density-based HLLC scheme of the present invention and the exact solution of the shock tube;

[0018] Figure 7 It is a geometric schematic diagram of the closed placement channel of the present invention;

[0019] Figure 8 It is the temperature cloud diagram of the explosion flame propagation of the present invention;

[0020] Figure 9 It is the comparison and verification of the propagation position of the flame front of the present invention. Detailed implementation manners

[0021] To enable those skilled in the art to better understand the solution of the present invention, the following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

[0022] It should be noted that the terms "first", "second", etc. in the description, claims and above-mentioned drawings of the present invention are used to distinguish similar objects, and do not necessarily describe a specific order or sequence. It should be understood that the data used in this way can be interchanged under appropriate circumstances, so that the embodiments of the present invention described here can be implemented in an order other than those illustrated or described here. In addition, the terms "comprising" and "having" and any variations thereof are intended to cover non-exclusive inclusion. For example, a process, method, system, product or device comprising a series of steps or units does not necessarily limit to those steps or units clearly listed, but may include other steps or units not clearly listed or inherent to these processes, methods, products or devices.

[0023] As Figures 1-9 shown, the present invention provides a calculation method for the combustible gas explosion process based on the reaction progress variable model. By constructing the source term of the reaction progress variable equation, the reaction progress variable is used to couple auto-ignition and flame propagation; by judging the critical turbulent Mach number, a smooth switch between the pressure-based solver and the density-based solver is realized, improving the calculation accuracy of the solver in the turbulent flame acceleration process of deflagration-to-detonation transition. The calculation accuracy of the explosion initiation process is improved through an ignition initiation model. Finally, actual explosion test cases are given to verify the applicability of the model. The specific implementation is as follows:

[0024] In premixed combustion, the flame propagation (controlled by transport) and auto-ignition (controlled by chemical reaction mechanism) phenomena are two of the most basic combustion mechanisms, which coexist and interact in the actual explosion process. Research shows that the combustible premixed gas needs to experience a slow reaction process at low temperature before auto-ignition. The auto-ignition model describing this process can obtain the auto-ignition delay time and reaction rate, so as to calculate the local temperature of the mixture. When the temperature is higher than 1000K, the reaction progress variable equation describing premixed flame propagation can be used to track the development of the flame front, thereby realizing the coupling between the auto-ignition model and the flame front model. The propagation of the flame front is simulated by solving the transport equation of the density-weighted average reaction progress variable c. Combustion is described by the reaction progress variable, c = 0 corresponds to the unburned mixture, and c = 1 corresponds to the fully burned mixture. The present invention uses a method for coupling the auto-ignition and flame propagation models based on the reaction progress variable, and this method describes the flame propagation process using the global reaction progress variable c. The reaction progress variable equation is:

[0025]

[0026] In the formula, P is the pressure; ρ is the density, R mix is the mixture gas constant; t is the time; x jis the spatial coordinate component in the j direction, and u j is the velocity component in the j direction, where j = 1, 2, 3 are spatial coordinate indices; the reaction progress variable c can be defined by any parameter (such as temperature T or reactant mass fraction Y f ). The superscripts "-" and "~" represent the physically filtered and mass-weighted Favre-averaged filtered quantities respectively. The source term of this equation consists of two subterms:

[0027]

[0028] In the formula, the two subterms respectively represent: ① generated by flame propagation modelled by the product of the flame propagation speed and the flame surface density; ② generated by auto-ignition calculated by the look-up table method.

[0029] In summary, the present invention combines the flame surface propagation model and the auto-ignition model through the reaction progress variable to successfully predict the hybrid combustion mode in which both flame surface propagation and auto-ignition coexist.

[0030] Generated by flame propagation modelled by the product of the flame surface propagation speed and the flame surface density, and considering the quenching effect of turbulent flames, the flame surface density equation is as follows:

[0031]

[0032] In the formula, 〈ρω〉 S is the average reaction rate per unit flame surface area; Σ is the flame surface density; ρ u is the density of the unburned mixture; S t is the turbulent combustion speed, which is modelled as the product of the laminar combustion speed S l and the flame wrinkling factor ξ Δ . This model is also called the flame thickening technique, which can accurately capture the flame surface on a low-resolution grid, reducing the computational cost due to the number of grids.

[0033] The quenching effect of turbulent flames is evaluated by the stretching factor G ∈ [0, 1], which represents the probability that the flame stretching will not extinguish; if there is no stretching (G = 1), the probability that the flame will not extinguish is 100%.

[0034]

[0035] In the formula, erfc(x) is the complementary error function, σ = 0.26·L / η is the standard deviation of turbulent dissipation, L is the turbulent integral length scale, η is the Kolmogorov microscale; ε is the turbulent dissipation rate; is the critical dissipation rate, and the critical strain rate is g cr can be obtained through experiments or through the model g cr = 0.5·S l 2 / α, where ν is the kinematic viscosity of the unburned mixture and α is the thermal diffusivity of the unburned mixture.

[0036] The sub-grid wrinkling factor model is a simplified method for evaluating the flame propagation state of premixed gases, developed using different semi-empirical correlations and the basic equations of gas dynamics. This method is based on the exponential fitting of the basic laminar burning velocity S l , the normalized dimensionless turbulent flame speed S t , the main dimensionless characteristics of the system, the reactivity and stability of the mixture.

[0037]

[0038] In the formula, u′ Δ is the sub-grid pulsation velocity; Δ is the sub-grid length; δ l is the laminar flame thickness; Le is the Lewis number; Re Δ is the sub-grid Reynolds number; Ka Δ is the sub-grid Karlovitz number. In actual modeling, several semi-empirical correlations are generally selected to cover a wider range of variables in Equation (5). Table 1 lists several representative wrinkling factor models.

[0039] Table 1 Several representative wrinkling factor models

[0040]

[0041] The premixed flame propagation speed S l can be obtained through the Metghalchi-Keck formula. The general expression between the laminar flame propagation speed S l in a closed space and the temperature T and pressure P is as follows (Molkov, 2000):

[0042]

[0043] where S l0 is the laminar flame propagation speed at normal temperature and pressure (NTP); α is the temperature exponent; β is the pressure exponent; in an adiabatic or isentropic system S l = S l0 (P / P0) ε , and the total thermokinetic exponent can be expressed as ε = α + β - α / γ, where γ is the adiabatic exponent of the unburned mixture. The above parameters can be obtained by referring to manuals or literature.

[0044] Taking premixed hydrogen-air as an example, according to the literature (Wiki and Jordan, 2007), the overall thermokinetic index (ε), laminar flame speed (S l0 ), and expansion coefficient (θ) are fitted to polynomial functions of the hydrogen volume fraction (x). As Figure 3 shown in their function images, the fitting curve expressions are as follows:

[0045]

[0046] As a preferred embodiment, in the actual combustion process, auto-ignition can be regarded as the spontaneous oxidation process of the fuel-air mixture. The challenge in modeling it lies in the prediction of the auto-ignition occurrence time and location. Currently, the most accurate method for predicting complex auto-ignition phenomena is to apply detailed chemical reaction models that include hundreds of intermediate components and reactions. These detailed reaction models can describe all auto-ignition phenomena within the entire range of temperature, pressure, and mixture fraction. However, with the existing CPU computing power, it is impossible to directly solve the transport of these chemical reactions and intermediate components in industrial applications. The look-up table method for chemical reactions solves this difficulty. The look-up table method uses a dedicated program to solve complex chemical reactions and stores their results on the computer. In the flow field solution, there is no need to recalculate the chemical reactions, and as long as the results are looked up in a certain way. This not only reduces the amount of computation but also ensures the accuracy of the chemical reactions used.

[0047] The chemical reaction table is based on the reaction mechanism of a homogeneous mixture of fuel and oxidizer under atmospheric pressure, and the reaction progress variable realizes the connection between the flow field and the look-up table. Once the local state is obtained in the flow field, its auto-ignition process can be obtained according to the look-up table. The early reaction process of auto-ignition is mainly the decomposition of the fuel without heat release. There is a time delay between the mixing of the fuel and oxidizer and the heat release, which must be considered in the model. To solve this problem, the auto-ignition delay is also stored in the calculated table, which avoids the problem of premature heat release in chemical reactions. To make the look-up table as simple as possible, the concept of ignition delay time is used in the present invention. The physical ignition delay time describing the auto-ignition process:

[0048] t ign = t ign (T, P, f) (12)

[0049] When the ignition delay time has passed, the ignition time variable τ = 1 (τ = t / t ign is the dimensionless time). As long as the ignition delay time has not been reached, it will not affect the flow characteristics. Only when τ ≥ 1 is reached, the mixture will be ignited. The auto-ignition source term can be expressed as:

[0050]

[0051] where Δt represents the current time step, and H(x) is the Heaviside function:

[0052]

[0053] The Heaviside function only activates the part of the computational cell where the local ignition delay time has expired, which can be achieved by weighting the autoignition source term with the volume fraction.

[0054] As a preferred embodiment, since autoignition in the DDT environment can be triggered by heating caused by shocks, a sub-model is introduced to improve the accuracy of autoignition modeling on coarse grids. Although the average temperature in the computational cell can be high enough to trigger autoignition, the detonation simulation must consider the fact that the shock that causes the temperature rise may not have passed through the entire computational cell. Therefore, a model that predicts the autoignition of the computational cell based on the average temperature and average pressure will lead to incorrect results. As Figure 4 shown, due to the huge difference in scales, the shock discontinuity cannot be resolved on the computational grid. Therefore, each computational cell is divided into two parts: in the part with volume fraction α, the temperature and pressure rise to T high and P high ; in the part with volume fraction 1 - α, the temperature and pressure remain at T low and P low . P high is defined by the highest value, and P low is defined by the lowest value that can be found in the surrounding computational cells. According to the consistency with the average pressure P in the computational cell, the volume fraction α can be determined:

[0055]

[0056] The flow field parameters before and after the shock are calculated from the dynamic shock relations of the ideal gas:

[0057]

[0058] When P high and P low are known, the shock Mach number Ma can be solved, and then the corresponding temperature ratio T high / T low can be obtained. This model is applicable not only to the computational cells with shocks but also to the entire computational domain. It should be noted that although the presence of the shock changes the temperature and pressure, the mixture fraction f(x) is not affected. The ignition delay time is evaluated separately on both sides of the shock wave:

[0059]

[0060] The volume-weighted autoignition source term can be expressed as:

[0061]

[0062] As a preferred embodiment, the wrinkling factor and ignition delay time in the flow field can be obtained through the following transport equation:

[0063]

[0064] In a non-uniform mixture, it is very important to correctly obtain the flame speed in the flow field. By evaluating the molar fraction x of the combustible gas based on the mixture fraction f(x) (volume fraction of fuel), the laminar flame propagation speed can be obtained. The spatial distribution of the mixture volume fraction is described by the transport equation:

[0065]

[0066] As a preferred embodiment, the flame acceleration process of deflagration-to-detonation belongs to a transonic reactive flow. In the initial stage of the explosion ignition start, the flame propagation speed is low and it belongs to an incompressible flow. At this time, using a pressure-based solver will obtain accurate calculation results. As the hydrodynamic instability induces flame acceleration, the flame changes from an incompressible flow to a compressible flow and shock waves are generated. At this time, it is more reasonable to use a density-based solver.

[0067] The evolution problem of the initial discontinuity of one-dimensional inviscid flow (Sod shock tube problem) belongs to a typical discontinuous problem and there is an exact solution, so it is widely used to compare and verify the accuracy of discrete formats and numerical methods in CFD. As a test case, the present invention compares the calculation results of the pressure-based PISO solver and the density-based HLLC solver in the 1D Sod shock tube problem. As Figure 5 shown, the medium inside the shock wave is an ideal gas with a molecular weight M = 28.85 g / mol and a specific heat ratio γ = 1.4. The initial conditions on the left side (x ≤ 0) of the shock tube are: P L = 1E6 Pa, T L = 1000K, U L = 0m / s; the initial conditions on the right side (x > 0) of the shock tube are: P R = 1E5 Pa, T R = 300K, U R = 0m / s. Figure 6 The results shown were obtained on a uniform grid with a spacing of 1 mm. It can be seen that the pressure-based PISO solver not only predicts a wrong shock position, but also the flow field variables obtained are very dissipative and show overshoots at the discontinuities. The density-based solver, especially the HLLC solver, shows better performance, producing accurate shock propagation speeds and less dissipation at the discontinuities.

[0068] It can be concluded that the standard pressure-based solver is not suitable for simulating rapid deflagration and explosion. The disadvantage of the density-based solver is that it cannot work in very low Mach number flows. Therefore, in the present invention, a pressure-based solver is used to start the calculation of the incompressible flow at the initial stage of the explosion, and then it is switched to the density-based solver once the combustion-induced flow develops. The turbulent critical Mach number Ma t = 0.1 is used to achieve a smooth transition between the two solvers; that is, the pressure-based calculator is used when Ma t ≤ 0.1, and the program switches to the density-based solver when Ma t > 0.1.

[0069] Preferably, the combustion is started at a set time and position in the combustion chamber, and this goal can be achieved by a high-pressure pulsed spark igniter. The spark ignition event occurs very briefly relative to the combustion process in the combustion chamber. The physical description of this simple event is very complex, so it is difficult to accurately simulate the ignition process. In addition, the energy generated by the spark event is several orders of magnitude smaller than the chemical energy released by the fuel combustion. Although a large amount of research has been devoted to the spark ignition theory and ignition devices, the ignition in the computational domain depends more on the mixture composition than on the spark energy (Heywood, 1988). Therefore, it is not necessary to model the spark event in great detail, and it can simply be set as the start of combustion within the continuous ignition time.

[0070] Generally, due to the relatively small initial ignition radius, its resolution on the CFD grid is insufficient, and ignition may not be initiated. Igniting by burning several grid cells around the ignition position will result in strong sensitivity to the grid and time step, leading to incorrect simulation results. Therefore, a suitable ignition start model is very important for the accurate simulation of the explosion process. To reduce the sensitivity to the ignition initial conditions, in the present invention, within the short ignition period t in the ignition cell, the progress variable of the unburned mixture in the ignition cell linearly increases from 0 to 1, and at the same time, the temperature in this region linearly increases from the initial temperature to the adiabatic flame temperature to initiate ignition. The ignition time t is defined as the time required for the laminar flame to propagate to half of the equivalent length of the ignition cell at the initial temperature and pressure.

[0071]

[0072] In the formula, θ = ρ u / ρ b is the expansion coefficient, where ρ u is the density of the unburned gas, ρ b is the density of the burned gas; S l is the laminar flame propagation speed. The equivalent length Δ of the fire cell in the computational domain can be expressed as a function of the control volume V c for evaluation:

[0073]

[0074] The present invention verifies the numerical algorithm by obtaining verification data through the DDT open online platform (http: / / www.td.mw.tum.de / ddt). The Gravent test bench (Boeck, 2016) is a closed channel with a rectangular cross-section of 5400 mm × 300 mm × 60 mm. Figure 7 shows the configuration and cross-sectional view of the channel. The channel consists of two parts: an obstacle area and a smooth area. The obstacle area is 2050 mm long, and the smooth area is 3350 mm long. Seven obstacles with a thickness of 12 mm and a blockage ratio of BR = 60% are installed in the channel. The first obstacle is located 250 mm from the ignition source, and the obstacle spacing is 300 mm. The channel is filled with a hydrogen-air mixture with a stoichiometric concentration of 25% and a vertical concentration gradient (diffusion time 3 s). The initial pressure and temperature of the mixture are 101325 Pa and 293 K respectively. As Figure 8 and 9 shown, it can be seen that the flame front position obtained by using the numerical algorithm of the present invention is in good agreement with the experimental results, and the maximum error does not exceed 5%.

[0075] The serial numbers of the above embodiments of the present invention are only for description and do not represent the advantages or disadvantages of the embodiments. In the above embodiments of the present invention, the descriptions of each embodiment have their own emphases. For the parts not detailed in a certain embodiment, reference can be made to the relevant descriptions of other embodiments. In the several embodiments provided in this application, it should be understood that the disclosed technical content can be implemented in other ways.

[0076] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those of ordinary skill in the art should understand that they can still modify the technical solutions recorded in the foregoing embodiments, or perform equivalent replacements on some or all of the technical features; and these modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A calculation method for the explosion process of combustible gas based on a reaction progress variable model, characterized in that, Including the following steps: S1: According to the construction of the source phase of the reaction progress variable equation, coupling auto-ignition and flame propagation through the reaction progress variable; S2: Achieving smooth switching between the pressure-based solver and the density-based solver by judging the critical turbulent Mach number, thereby improving the calculation accuracy of the solver in the turbulent flame acceleration process of deflagration-to-detonation transition; S3: Improving the calculation accuracy of the explosion initiation process by establishing an ignition initiation model; The ignition initiation model initiates ignition by linearly increasing the progress variable of the unburned mixture in the ignition unit from 0 to 1 within a short ignition period t in the ignition unit, and simultaneously linearly increasing the temperature in this region from the initial temperature to the adiabatic flame temperature; the ignition time t is defined as the time required for the laminar flame to propagate to half of the equivalent length of the ignition unit at the initial temperature and pressure; where θ = ρ u / ρ b is the expansion coefficient, where ρ u is the density of unburned gas, and ρ b is the density of burned gas; S l is the laminar flame speed; the equivalent length Δ of the fire cell in the computational domain is expressed as a function of the control volume V c for evaluation:

2. The calculation method of the combustible gas explosion process based on the reaction progress variable model according to claim 1, characterized in that During the coupling process, the flame propagation process is described by the global reaction progress variable c; The transport equation of the reaction progress variable is: wherein, P represents pressure; ρ represents density, R mix represents the mixture gas constant; t represents time; x j represents the spatial coordinate component in the j direction, u j represents the velocity component in the j direction, and j = 1, 2, 3 are spatial coordinate indices; the reaction progress variable c is defined according to the temperature T or the mass fraction Y of the reactant f for definition, The superscripts "-" and "~" respectively represent the physical space filtering and the mass-weighted Favre averaging filtered quantity; the source term of the transport equation includes: Among them, represents the modeling of the product of flame propagation speed and flame surface density, which is generated by flame propagation; Generated by auto-ignition and calculated by the look-up table method.

3. The calculation method of the combustible gas explosion process based on the reaction progress variable model according to claim 1, characterized in that, Combining the flame surface propagation model and the auto-ignition model through the reaction progress variable, Generated by flame propagation Modelled by the product of the flame surface propagation speed and the flame surface density, and considering the quenching effect of the turbulent flame, the flame surface density equation is as follows: where <ρω> S represents the average reaction rate per unit flame surface area; Σ represents the flame surface density; ρ u represents the density of the unburned mixture; S t represents the turbulent burning velocity, which is modeled as the product of the laminar burning velocity S l and the flame wrinkling factor ξ Δ ; The turbulent flame quenching effect is evaluated by the stretch factor G ∈ [0, 1], and the stretch factor represents the probability that the flame stretch will not extinguish; if there is no stretch factor, i.e., G = 1, the probability that the flame will not extinguish is 100%; wherein, represents the complementary error function; σ = 0.26·Lη is the standard deviation of turbulent dissipation, L is the integral length scale of turbulence, η is the Kolmogorov microscale; ε is the turbulent dissipation rate; is the critical dissipation rate, and the critical strain rate g cr is obtained through experiments or through the model g cr = 0.5·S l 2 / α, where ν is the kinematic viscosity of the unburned mixture and α is the thermal diffusivity of the unburned mixture; The sub-grid wrinkling factor model is a simplified method for evaluating the propagation state of premixed gas flames, developed using different semi-empirical correlations and the basic equations of gas dynamics; based on the fundamental laminar burning velocity S l , the normalized dimensionless turbulent flame velocity S t , the dimensionless characteristics of the system, and the exponential fitting of the mixture reactivity and stability; where, u′ Δ represents the sub-grid fluctuating velocity; Δ represents the sub-grid length; δ l represents the laminar flame thickness; Le is the Lewis number; Re Δ is the sub-grid Reynolds number; Ka Δ is the sub-grid Karlovitz number; Furthermore, the premixed flame propagation speed S is obtained through the Metghalchi-Keck formula l ; the expression between the laminar flame propagation speed S l in a closed space and the temperature T and pressure P is as follows: where S l0 is the laminar flame speed under normal temperature and pressure (NTP); α is the temperature exponent; β is the pressure exponent; in an adiabatic or isentropic system, S l = S l0 (P / P0) ε , and the total thermokinetic exponent is expressed as ε = α + β - αγ, where γ is the adiabatic exponent of the unburned mixture.

Images