Turbulent flow-chemical reaction decoupling combustion modeling method

Abstract

The invention provides a turbulent flow-chemical reaction decoupling modeling method. The turbulent flow-chemical reaction decoupling modeling method comprises the following steps: step 1, basic premise and assumption; step 2, treatment of different types of reactions; according to the law of mass action, only single-molecule and double-molecule reactions are considered; step 3, correcting temperature distribution of the reaction rate in the grid; and non-uniform temperature distribution correction is introduced to obtain an average rate constant # imgabs0 #, and # imgabs1 # is used for chemical source item calculation. The method has the advantages that when a high-resolution multi-step skeleton chemical mechanism is applied, the method has the advantage of predicting reaction source items under all combustion modes and working conditions, and the calculation cost can be greatly reduced. The turbulence effect is considered in source item calculation of an existing turbulence combustion model, and the existing turbulence combustion model is limited to a turbulence model, so that the existing turbulence combustion model has relatively good prediction precision only for partial combustion modes.

Description

Technical Field

[0001] The invention relates to the technical field of turbulent flow numerical simulation, and in particular to a combustion modeling method for turbulent flow-chemical reaction decoupling. Background Art

[0002] In the context of Large Eddy Simulation (LES) and Finite Rate Chemistry (FRC) modeling, the chemical reaction rate of species i in turbulent combustion is a function of temperature T, pressure p and species X. i The mass fraction y i A strong nonlinear function, using To express. Among them, is the number of components X i The net reaction rate of all elementary reaction steps, y i represents the mass fraction of all species. In actual combustion numerical simulation calculations, all variables for solving conservation equations exist in filtered form. Therefore, the filtered chemical reaction source term is not equivalent to the net reaction rate calculated based on the filtered temperature, pressure and filtered species mass fraction, that is, To solve this problem, relevant researchers have proposed a series of turbulent chemical interaction models to predict the source terms of chemical reactions after filtration.

[0003] The Eddy Break-Up Model (EBU) assumes that the combustion reaction occurs in small-scale vortex structures, which are constantly broken and mixed during turbulent motion. When the characteristic time scale of the flow is comparable to the characteristic time scale of the chemical reaction, the combustion chemical reaction will be triggered. Therefore, the EBU determines the chemical reaction rate by comparing these two time scales. The FRC direct modeling method assumes that the reaction is carried out in a laminar state and ignores the turbulent chemical interaction (TCI). The chemical reaction rate is directly calculated using the Arrhenius equation. The Eddy-Dissipation Concept Model (EDC) assumes that the flow within a grid can be divided into large-scale turbulent structures and small-scale fine structures, and that chemical reactions occur in small-scale turbulent structures, which are constantly generated and dissipated during the cascade transfer of turbulent energy. The Partially Stirred Reactor (PaSR) assumes that the fluid in the reaction area is not completely uniformly mixed, but is partially stirred. This means that in the reactor, there are both areas with sufficient stirring, so that the reactants can contact and react well; and areas with insufficient stirring, where the degree of mixing of the reactants is relatively low. The reaction rate is determined by the chemical reaction kinetics and turbulent mixing. The flame surface model (Flamelet-like Model) decomposes the complex turbulent combustion problem into two relatively simple sub-problems: the turbulent flow problem and the laminar flame structure problem. It assumes that the turbulent flame is composed of a series of laminar flame surfaces (Flamelet, also known as small flames). These laminar flame surfaces are constantly stretched, twisted and wrinkled under the action of turbulence, but at each local position, the combustion process can still be described by the characteristics of laminar flames. The probability density function (PDF) method is a method for dealing with the probability distribution of physical quantities in turbulent combustion. In the turbulent combustion process, due to the randomness and complexity of turbulence, the physical quantities in the flow field have the characteristics of random fluctuations. The purpose of the PDF method is to simulate the turbulent combustion process more accurately by describing the probability distribution of these physical quantities.

[0004] Current turbulent combustion interaction models are often similar to sub-grid modeling of turbulence. People try to introduce turbulent characteristic parameters to consider the influence of small vortices on reactions within the grid. Solving the control equation can only give the average species concentration of the grid, and cannot analyze the characteristics of small vortices within the grid. In fact, even for small vortices, the vortex scale is much larger than the molecular scale of chemical reactions, and the chemical reaction within the vortex depends on molecular diffusion and collision. What is needed to determine the chemical reaction rate within the grid, that is, the chemical source term, is the spatial distribution of temperature and concentration at the sub-grid scale. Introducing pulsating values ​​to improve the modeling of chemical source terms may lead to the complexity of mathematical modeling and the introduction of unnecessary turbulence parameters, leading to the complication of the problem. The present invention uses the average species concentration and temperature of the flow field obtained by the species transport equation to establish a chemical source term calculation model to achieve decoupled calculation of flow and chemical kinetics. In the case of non-uniform mixing, the non-uniformity of the spatial distribution of concentration is introduced in the concentration product calculation of the bimolecular reaction to overcome the error caused by the uniform mixing assumption. According to the average species concentration and the "elementary" steps of the mechanism and the number of reaction molecules, the concentration product spatial integral calculation method under the assumption of linear spatial distribution of concentration is established (this method is also applicable to other concentration spatial distribution situations), and then the chemical source term is obtained. On the other hand, when the temperature is non-uniformly distributed inside the grid, the average temperature will underestimate the reaction rate. When the activation energy is high, it may lead to a large error in the calculation of the rate constant. Therefore, the weighted average is used in combination with the temperature of the adjacent grid to obtain the rate constant value rate constant correction. This treatment adopts the spatial distribution modeling of species concentration and temperature in the grid to avoid the introduction of turbulence parameters in the source term calculation process, thereby realizing the decoupled calculation of chemical kinetics and turbulence.

[0005] A commonly used model is the PaSR model, which is characterized by dividing the inner part of each control volume into a reaction zone and a mixing zone, and assuming that the chemical components in the reaction zone have been completely mixed and chemical reactions have occurred therein, while no chemical reactions have occurred in the mixing zone. The volume fraction κ of the control volume occupied by the reaction zone is used to approximate the left side of the right side of the formula. The existing PaSR model uses the characteristic time scale of the chemical reaction of a component and the characteristic time scale of the mixing of a component to calculate the volume fraction κ as follows:

[0006]

[0007] Among them, τ c is the characteristic time scale of the reaction, τ m is the characteristic time scale.

[0008] In the current PaSR model, the calculation of turbulent combustion often only uses a simple total package reaction mechanism, that is, the existing PaSR model has good accuracy when facing simple reaction processes. However, when facing turbulent combustion calculations involving complex chemical reactions, the accuracy is not high enough, that is, when the existing PaSR model needs to use a high-resolution skeleton chemical mechanism, the chemical characteristic time is difficult to determine and the accuracy is difficult to guarantee.

[0009] Another commonly used turbulent combustion model is the EDC model, which divides the flow field into two regions: large-scale turbulent eddies and small-scale micro-mixing zones. In large-scale turbulent eddies, the motion and properties of the fluid are described by traditional turbulence models (such as the k-∈ model). Chemical reactions occur in small-scale turbulent structures, which are continuously generated and dissipated during the cascade transfer of turbulent energy. The source term of the component conservation equation is expressed as:

[0010]

[0011] Among them, the length scale fraction Characteristic time scale is the mass fraction of species α in the fine structure, is the mass fraction of the fluid species α outside the fine structure.

[0012] Current EDC models often require the use of simplified chemical reaction mechanisms because the models assume that chemical reactions reach equilibrium in a very short time, which may not be applicable to all types of combustion reactions. The accuracy of EDC models is highly dependent on the description of turbulent mixing. If the mixing model is not accurate enough, the simulation results of the chemical reactions will also be affected. EDC models are suitable for fast chemical reactions and may not be accurate enough for some slow reactions or when detailed chemical kinetics need to be considered. Although EDC models are computationally cheaper than some more complex combustion models (such as PDF transport models), they still require high computing resources, especially when dealing with large-scale or complex geometries.

[0013] There is another model called CMC model, the core idea of ​​which is to condition the concentrations of reactants and products in the flow field. This means that the model does not directly deal with the statistical characteristics of the entire flow field, but focuses on the local flow field characteristics under given conditions (such as under specific mixing fractions or temperature conditions). In the CMC model, the main thing to solve is the conditional moment equations, which describe the statistical characteristics of thermodynamic variables such as species concentration and temperature under given conditions. The coupling of turbulence and chemical reactions is achieved in the CMC model through conditional expectation values. The model assumes that under given conditions, the chemical reaction rate can be obtained by averaging, while the influence of turbulence is described by a turbulence model.

[0014] Although the CMC model can provide relatively accurate simulation results, its computational cost is relatively high. This is mainly because it is necessary to solve multiple conditional moment equations, each of which may involve complex chemical reactions and turbulence models. The accuracy of the CMC model depends largely on the accuracy and complexity of the chemical mechanism. For complex chemical mechanisms with a large number of reactants and reaction pathways, the computational burden of the CMC model will increase significantly. The CMC model is very sensitive to model parameters (such as mixing time scale, initial conditions, etc.), which may lead to uncertainty in model predictions. The performance of the CMC model depends largely on the selected turbulence model. If the turbulence model itself is not accurate, the prediction accuracy of the entire CMC model will also be affected. Summary of the invention

[0015] The present invention aims to solve the technical problem of decoupling calculation of chemical kinetics and turbulence in the prior art when facing turbulent combustion calculation involving complex chemical reactions, and to provide a calculation method and model for decoupling turbulence and reaction in the context of high-resolution chemical mechanism application.

[0016] In order to solve the above technical problems, the turbulent flow-chemical reaction decoupling modeling method proposed in the present invention comprises the following steps:

[0017] Step 1: Basic premises and assumptions;

[0018] (1) The chemical reaction rate depends on the species concentration and temperature distribution within the grid;

[0019] (2) The turbulence in the grid cannot be resolved, and the use of a method similar to small eddy modeling cannot solve the problem that the chemical source term cannot be calculated using the average concentration;

[0020] (3) Sub-grid vortex modeling provides the reaction temperature and pressure of grid G;

[0021] (4) Obtain the grid average concentration by solving ODE based on CKL or other chemical reaction mechanisms;

[0022] (5) Accurate calculation of single-molecule reaction rate equations;

[0023] (6) The bimolecular reaction concentration product of the reaction rate equation adopts the linear species concentration assumption. The concentration distribution curve of the grid G ​​is determined according to the average concentration of the adjacent grid points G-1 and G+1. The concentration product in three directions is spatially integrated to obtain the reaction rate equation.

[0024] (7) The trimolecular reaction and Lindemann unimolecular reaction are degenerated into bimolecular and unimolecular reactions respectively by assuming the three-body concentration is constant;

[0025] (8) For the total package reaction with non-integer order, the reaction rate is calculated by using the double midpoint concentration approximation of the linear spatial distribution approximation;

[0026] (9) Directly obtain species X in grid G ​​by summing up the multi-step reactions of the reaction rate equation i of chemical sources.

[0027] The rate of change of species concentration is as follows:

[0028]

[0029] Among them, y i For species X i The mass fraction, c i is the molar concentration, M i is the molecular weight and ρ is the density.

[0030] For the reaction of N species L, the reaction mechanism of the elementary step list is:

[0031]

[0032] According to the law of mass action, species X i The rate of generation in the rth reaction is but

[0033]

[0034] Step 2: Processing of different types of reactions;

[0035] According to the law of mass action, only unimolecular and bimolecular reactions are considered.

[0036] (1) Single molecule reaction:

[0037] Cracking and isomerization reactions fall into this category:

[0038]

[0039] Where P is the product. The rate of change of the number of molecules of A in the grid is only related to the total number of A and has nothing to do with the spatial distribution of concentration. The average concentration obtained by solving the ODE is The reaction rate is calculated accurately and has nothing to do with the flow, that is:

[0040]

[0041] (2) Lindemann unimolecular reaction mechanism,

[0042] A * →P

[0043] The true reverse reaction rate constants of the first reversible reaction are k1 and k -1The rate constant of the second one-way reaction is k2. According to the first step fast equilibrium approximation, the quasi-steady-state species A can be obtained. * The concentration of P can be expressed as

[0044]

[0045] Where c M is the three-body concentration, that is, the sum of the concentrations of all species, which can be regarded as a constant. The above formula degenerates into the case of a single molecule reaction, that is,

[0046]

[0047] (3) Bimolecular reaction:

[0048] Shape

[0049] The reaction rate of a bimolecular reaction is written as

[0050]

[0051] Such reactions include oxidation, free radical recombination reactions, etc. The contribution of species concentration to the chemical source term depends on the concentration product factor c A c B ;

[0052] When A and B are uniformly distributed in the grid volume V, that is, at any position,

[0053]

[0054] but

[0055]

[0056] Of the two species included in the concentration product, if one has a uniform spatial distribution of concentration within the grid G, equation (13) holds true;

[0057] The reaction rate is calculated by setting the spatial distribution of the concentration of each species in the grid. When the concentration distribution of A or B in the grid is uneven, it is necessary to set the spatial distribution of the concentration. In the actual flow field calculation, the average concentration of grid G ​​and the adjacent grids G-1 and G+1 is obtained. First, a one-dimensional grid in the x direction is discussed. The concentration is assumed to change linearly along the x direction. The average concentration of the left and right adjacent grids G-1 and G+1 is taken as the actual concentration of the corresponding boundary of grid G. In order to ensure the conservation of the number of species A particles in grid G, that is, the number of moles of species A in the grid is unchanged, where V is the volume of the grid G; two methods are used to represent the concentration distribution in the grid. First, consider the linear concentration distribution case, let the one-dimensional concentration distribution, straight line L1 is the straight line connecting the left and right boundaries of the grid G, and the concentrations at the left and right boundaries are set to the average concentrations of the left adjacent grid G-1 and the right adjacent grid G+1 respectively and In order to ensure the conservation of the number of moles of species A in the grid, the straight line is translated, and the concentration at x = 0.5 is In this way, the linear distribution of the spatial concentration of species A in grid G ​​is represented by straight line L2.

[0058]

[0059] Consider two chemical species A and B.

[0060]

[0061] Make the same assumptions in the y and z directions, then:

[0062]

[0063] Calculation of reaction rate when considering the concentration of two adjacent grids under one-dimensional conditions. Integrate along the x-axis within the grid, and integrate the concentration product of the bimolecular reaction rate along the x-direction within the grid G, that is,

[0064]

[0065] Substituting equations (12) and (13) into equation (18), we get the reaction rate under one-dimensional conditions:

[0066]

[0067] exist or When , the concentration product is reduced to the case of uniform concentration distribution or unimolecular reaction, that is,

[0068]

[0069] The reaction rate is calculated by averaging the integration results of the three one-dimensional concentration products, and the expression of the chemical source term is:

[0070]

[0071] The summation term is the summation in the x, y, and z directions, and m is the number of summation terms;

[0072] (4) Trimolecular reaction and overall mechanism

[0073] The trimolecular reaction contains one species M, whose concentration is c M is the trisomic concentration, cM is the sum of the concentrations of all species and can be assumed to be constant. The rate equation of the trimolecular reaction degenerates into a bimolecular reaction, and the rate constant is treated as a bimolecular reaction, i.e., equation (21).

[0074] Complex reaction mechanisms also have a general package mechanism, and the reaction rate is as follows

[0075]

[0076] The total package reaction is a form of reaction rate equation with relatively low accuracy. The reaction rate can be calculated by taking the average of multiple average concentrations to simplify the calculation. Considering the linear concentration distribution assumption, the average concentration of x = 0 to 0.5 is taken in the grid G. and the average value of x = 0.5 to 1.0 The simple average of c A，x ,Right now

[0077]

[0078] A similar treatment to equation (21) can yield the reaction rate equation when considering the m-dimensional concentration spatial distribution.

[0079] The chemical source term is obtained by summing the reaction rates of all steps r in the generation and consumption of species A in the reaction mechanism, that is:

[0080]

[0081] Step 3, temperature distribution correction of reaction rate within the grid;

[0082] Introducing non-uniform temperature distribution correction to obtain the average rate constant and will Used for chemical source term calculations. For bimolecular reactions, the reaction rates are calculated as

[0083]

[0084] The rate constant k is expressed by the Arrhenius equation. a When the average temperature is Substituting the Arrhenius equation will increase the deviation Right now

[0085]

[0086] For the convenience of mathematical processing, according to Arrhenius's differential form:

[0087]

[0088] When the activation energy E a≈0, the effect of T on the rate constant k value weakens. At this time, the average temperature Substituting the reaction rate equation into the calculation of the reaction rate will not lead to a large error;

[0089] In E a >20kcal / mol, the effect of uneven temperature distribution on the reaction rate calculation needs to be considered, and the second-order derivative

[0090]

[0091] The inflection point temperature of the k~T curve is If the activation energy E a =20 kcal·mol -1 ,T inf is 5000K, and the visible grid G The value is much lower than this temperature. It is easy to find that the average temperature below the inflection point temperature Chemical reaction rates will be underestimated. The present invention considers the method of weighted average, assuming The corresponding equivalent temperature is T b ,Due to the exponential dependence of k on T, T is obtained by weighted average in the following way b ,Right now:

[0092]

[0093] Where T x- and T x+ are the left and right boundary temperatures of grid G, respectively. The average temperature of grids G-1 and G+1 can be taken. x+ The weighting factor is taken as Logarithmic term Middle molecular To ensure E a =0 o'clock No divergence and weighting factor Get T b To ensure

[0094]

[0095] Weighted average to obtain equivalent temperature T b , we can get the rate constant term. In addition to considering the average temperature of the two adjacent grids G-1 and G+1, we also need to consider the average temperature of the grid G ​​itself. Take T b and The chemical reaction rate is calculated by taking the arithmetic mean of

[0096]

[0097] Where c A c B It is obtained by the concentration product calculation method. a ≈0, simply take

[0098] When considering a two-dimensional or three-dimensional grid, take: In order to consider the dimension of the actual temperature distribution in space, the chemical source term for the correction of temperature inhomogeneity is obtained by substituting it into equation (35).

[0099] The advantage of the present invention is that it has an advantage in predicting the reaction source terms under all combustion modes and conditions when applying high-resolution multi-step skeleton chemical mechanisms. Existing turbulent combustion models all take turbulence effects into account in the source term calculations, which are limited to turbulence models and therefore only have good prediction accuracy for some combustion modes. BRIEF DESCRIPTION OF THE DRAWINGS

[0100] Figure 1 It is a schematic diagram of one-dimensional linear concentration distribution of the present invention;

[0101] Figure 2 Calculate the spatial distribution of linear concentration within the region;

[0102] Figure 3 The average concentration of each grid point when the grid is divided into 10 equal parts. DETAILED DESCRIPTION

[0103] The specific technical solution of the present invention is described in conjunction with the accompanying drawings:

[0104] A turbulent flow-chemical reaction decoupling modeling method includes the following steps.

[0105] Step 1: Basic premises and assumptions;

[0106] (1) The chemical reaction rate depends on the species concentration and temperature distribution within the grid;

[0107] (2) The turbulence within the grid cannot be resolved, and the use of methods such as small eddy modeling cannot solve the problem of the spatial distribution of species concentration within the grid;

[0108] (3) Sub-grid vortex modeling provides the reaction temperature and pressure of grid G;

[0109] (4) Obtain the grid average concentration based on the ODE solution of the chemical reaction mechanism;

[0110] (5) Accurate calculation of single-molecule reaction rate equations;

[0111] (6) The bimolecular reaction concentration product of the reaction rate equation uses a linear concentration approximation (other concentration spatial distribution forms can also be used), and the concentration distribution curve of the grid G ​​is determined based on the average concentration of adjacent grid points; the concentration product in three directions is spatially integrated to obtain the reaction rate equation and chemical source term;

[0112] (7) The trimolecular reaction and the Lindemann unimolecular reaction are degenerated into bimolecular and unimolecular reactions respectively by assuming the three-body concentration is constant;

[0113] (8) For the total package reaction with non-integer order, the double midpoint concentration approximation is used to calculate the reaction rate equation of concentration change;

[0114] (9) The chemical source terms of species in the grid G ​​are directly obtained by summing up the multi-step reactions of the reaction rate equation.

[0115] Step 2: Processing of different types of reactions;

[0116] According to the law of mass action, the present invention only considers unimolecular and bimolecular reactions.

[0117] (1) Single molecule reaction

[0118] Many cracking and isomerization reactions belong to this category.

[0119]

[0120] Where P is the product. The rate of change of the number of molecules of A in the grid is only related to the total number of A and has nothing to do with the spatial distribution of concentration. The average concentration obtained by solving the ODE is The reaction rate is calculated exactly and is independent of the flow, i.e.

[0121]

[0122] (2) Lindemann single molecule reaction mechanism

[0123]

[0124] c M c A Calculation of concentration product, c M is the trisomic concentration, i.e., the sum of the concentrations of all species, which can be regarded as c M is a constant, so the above equation degenerates into the case of a single molecule reaction.

[0125]

[0126] (3) Bimolecular reaction

[0127] Such as: The reaction rate is generally written as:

[0128] Such reactions include oxidation, free radical recombination, etc. The contribution of species concentration to the chemical source term depends on the concentration product factor c A c B .

[0129] When A and B are uniformly distributed in the grid volume V, that is, at any position:

[0130]

[0131] but

[0132]

[0133] If one of the two species included in the concentration product is uniformly distributed in the grid G, that is, satisfies equation (4), then equation (5) holds.

[0134] In general, bringing the average concentration into the calculation may overestimate the reaction rate. In order to achieve the decoupled calculation of the chemical source term and the flow, the present invention sets the spatial distribution of the concentration of each species in the grid to calculate the reaction rate. When the concentration distribution of A or B in the grid is uneven, it is necessary to set the spatial distribution of the concentration. In the actual flow field calculation, the average concentration of the grid G ​​and the adjacent grids G-1 and G+1 can be obtained. First, take the one-dimensional grid in the x direction for discussion, assume that the concentration changes linearly along the x direction, and take the average concentration of the left and right adjacent grids G-1 and G+1 as the actual concentration of the corresponding boundary of the grid G. In order to ensure that the species X of the grid G i The number of particles is conserved, that is, the species X in the grid i The number of moles is unchanged, where V is the volume of the grid G ​​(or the length in the one-dimensional case). Assume that all species X i The spatial concentration distribution is linear, and the one-dimensional concentration distribution is as follows: Figure 1 As shown, Figure 1 The middle straight line L1 is a straight line connecting the left and right boundaries of the grid G. The concentrations at the left and right boundaries are set to the average concentrations of the left neighboring grid G-1 and the right neighboring grid G+1 respectively. and To ensure that the X i The number of moles of the species is conserved, so the straight line is shifted and the concentration at x = 0.5 is In this way, the straight line L2 represents the spatial concentration distribution of species in the grid G. Considering two species A and B,

[0135]

[0136] The same assumption can be made in the y and z directions. Here, x, y, and z do not have a strict rectangular coordinate system concept. The following discusses the reaction rate calculation when considering the concentration of two adjacent grids under one-dimensional conditions. Integrate along the x-axis within the grid and get:

[0137]

[0138] Substituting equations (6) and (7) into equation (8), we can obtain the reaction rate under one-dimensional conditions:

[0139]

[0140] It can be seen that in or When , the concentration product is reduced to the case of uniform concentration distribution or unimolecular reaction, that is,

[0141]

[0142] Considering the case of more adjacent grids in three dimensions, the average concentration product is simply replaced by the average of the integral results of three one-dimensional concentration products, and the expression of the chemical source term is obtained as follows:

[0143]

[0144] The summation term is the summation in the x, y, z directions, etc., and m is the number of summation terms. It can be seen that in a bimolecular reaction, when the concentration of any component is uniformly distributed in the grid, the reaction rate can be calculated using the concentration product of the average concentration.

[0145] (4) Trimolecular reaction and overall mechanism

[0146] There are a few trimolecular reactions in the multi-step reaction mechanism. This trimolecular reaction generally contains an M species with a concentration of c M is the trisomic concentration, since c M is the sum of the concentrations of all species and can be approximately considered unchanged. Therefore, the rate equation of the trimolecular reaction degenerates into a bimolecular reaction, and the rate constant can be treated as a bimolecular reaction (Equation 11).

[0147] Complex reaction mechanisms also have general package mechanisms. Due to the complexity of chemical reactions, the reaction order of this type of mechanism is determined by experiments or other methods. It is generally not an integer, and the reaction rate may be in the form of:

[0148]

[0149] At this time, it becomes more complicated to use the concentration product to integrate in the grid space to determine the chemical source term. At this time, a simpler method can be considered. For example, similar to the linear concentration assumption of bimolecular reaction, the average of several average concentrations can be taken to calculate the reaction rate instead of integrating in the grid. For example, the average concentration of x = 0.25 and x = 0.75 is taken on the L2 line, that is:

[0150]

[0151] Similar treatment of bimolecular reactions can be used to obtain the reaction rate equation when considering an m-dimensional grid. i The chemical source term can be obtained by summing the reaction rates of all steps r of generation and consumption, that is:

[0152]

[0153] Reaction rate calculation using hydrogen and oxygen combustion as an example

[0154] Consider an induction step of hydrogen-oxygen combustion,

[0155] like Figure 2 , the distance between the two nozzles is a, and a one-dimensional grid is taken, with a grid length of a. Assume that the initial value of H2 and O2 in the one-dimensional grid is 1 mol, and the concentration of H2 (mol per unit length) decreases linearly from 2 mol / a at x = 0 to 0 at x = a. On the contrary, the concentration of O2 increases linearly from 0 at x = 0 to 2 mol / a at x = a. Then the average concentration

[0156] By definition, If the average concentration of a one-dimensional grid is used to represent the reaction rate of this step (only the forward direction is considered), when x = 0-a, only one grid point is used for the distance, and the average concentration is used to calculate the concentration product c H2 c O2 get

[0157]

[0158] The subscript “1” of B1 indicates that the length 0-a is divided into only one grid point.

[0159] If the length 0-a is divided into 10 grids, each grid length is a / 10. The average concentration of each grid point is Figure 3 The average concentration product of 10 grid points is shown in Figure 2.

[0160]

[0161] Similarly, when 0-a is divided into 100 grid points,

[0162] According to the concentration integration method of the present invention, Figure 2 The linear concentration assumption is

[0163]

[0164] Substituting the integral into

[0165]

[0166] B ∞ This is the exact value of the reaction rate obtained by integrating the concentration spatial distribution function.

[0167] B1>B 10 >B 100 >B ∞ (22) When the reaction rate is calculated using the grid average concentration with different grid points, the relative error with the exact value of the reaction rate (assuming a linear concentration distribution) is as follows:

[0168]

[0169] Similarly, δ 10 =0.5%,δ 100 =0.005%. It can be seen that as the mesh becomes denser, the error caused by simply calculating the concentration product using the average concentration decreases rapidly.

[0170] Step 3, temperature distribution correction of reaction rate within the grid;

[0171] The value of the chemical source term depends only on the temperature and concentration in the grid. In large eddy simulation or RANS numerical calculation, solving the energy equation can obtain the average temperature of the grid, but at locations such as the flame surface, the temperature distribution in the grid may be very uneven. At this time, the impact on the calculation of the chemical source term needs to be considered. In particular, when the thickness of the flame surface is comparable to the grid scale, using the average temperature to calculate the chemical source term may bring large errors, and it is necessary to introduce non-uniform temperature distribution correction. A feasible solution is to obtain the average rate constant and will Used for chemical source term calculations. For example, for bimolecular reactions, the reaction rate can be calculated as:

[0172]

[0173] However, in general, especially when the activation energy E a When the average temperature is larger Substituting into the Arrhenius equation does not give Right now:

[0174]

[0175] For the convenience of mathematical processing, according to Arrhenius's differential form:

[0176]

[0177] It can be seen that the uniformity of the spatial distribution of temperature in the grid (assuming one dimension) has an effect on the calculation of the chemical source term and the activation energy E a When the activation energy is low (E a ≈0), the effect of T on the rate constant k value weakens. At this time, the average temperature Substituting the rate constant into the reaction rate equation will not lead to large errors. a From the differential form of the Arrhenius equation, it can be seen that different values ​​of T will have a great impact on the calculation of k, even an order of magnitude difference. a When it is larger (such as E a >20kcal / mol), the influence of uneven temperature distribution on the calculation of chemical source terms needs to be considered.

[0178]

[0179] It can be obtained that the inflection point temperature of the k~T curve is If the activation energy E a =20kcal / mol,T inf is 5000K. The average grid temperature It is easy to find that the average temperature is used below the inflection point temperature. Chemical reaction rates will be underestimated. In order to obtain reasonable The present invention considers the weighted average method, assuming The corresponding equivalent temperature is T b ,Due to the exponential dependence of k on T, when the gradient of temperature distribution in the grid is large, considering the exponential nature of the effect of temperature on the rate constant, T is obtained by weighted average in the following way: b ,Right now:

[0180]

[0181] Where T x- and T x+ are the left and right boundary temperatures of grid G, respectively, and can be taken as the average temperature of grids G-1 and G+1. x+ The weighting factor is taken as Logarithmic term Middle molecular To ensure E a =0 o'clock No divergence and weighting factor Formula (30) is an empirical formula, and its reliability is easy to verify. Other methods can also be used to obtain T b To ensure:

[0182]

[0183] Table 1 shows the average rate constant calculated according to formula (30) when the temperature in the grid changes linearly from 300K to 2400K under different activation energy conditions.

[0184] Table 1 Calculated by weighted average method

[0185]

[0186] Based on the weighted average of the above, the equivalent temperature T is obtained b In the actual calculation, T can be simply taken as b and the average temperature of the grid G The chemical reaction rate is calculated by taking the arithmetic mean of

[0187]

[0188] Where c A c B Corrected by the above concentration product calculation method. a ≈0, we can simply take

[0189] When considering a two-dimensional or three-dimensional grid, a similar treatment to the concentration product is possible: In order to consider the dimension of the actual temperature distribution in space, the chemical source term for the correction of temperature inhomogeneity can be obtained by substituting it into equation (32).

Claims

1. A turbulent flow-chemical reaction decoupling modeling method, characterized in that: The following steps are involved: Step 1: Basic premises and assumptions; Step 2: Processing of different types of reactions; According to the law of mass action, only unimolecular and bimolecular reactions are considered; Step 3, temperature distribution correction of reaction rate within the grid; Introducing non-uniform temperature distribution correction to obtain the average rate constant and will Used for chemical source term calculations.

2. According to the turbulent flow-chemical reaction decoupling modeling method of claim 1, step 1 comprises: (1) The chemical reaction rate depends on the spatial distribution of species concentration and temperature within the grid; (2) The turbulence within the grid cannot be resolved, and the spatial distribution of species concentration within the grid cannot be obtained using a processing method similar to small eddy modeling; (3) Sub-grid vortex modeling provides the reaction temperature and pressure of grid G; (4) Obtain the grid average concentration by solving ODEs based on CKL (a minimization reaction network mechanism developed by the Combustion Dynamics Center of Sichuan University) or other chemical reaction mechanisms; (5) Accurate calculation of single-molecule reaction rate equations; (6) The bimolecular reaction concentration product of the reaction rate equation adopts the spatial linear distribution and the concentration spatial distribution approximation, and the concentration distribution curve of the grid G ​​is constructed according to the average concentration of the adjacent grids; the concentration product in three directions is spatially integrated to obtain the concentration factor of the reaction rate equation; (7) The trimolecular reaction and Lindemann unimolecular reaction are degenerated into bimolecular and unimolecular reactions respectively by assuming the three-body concentration is constant; (8) The concentration factor of the reaction rate equation for the total package reaction of non-integer order is calculated using a double midpoint concentration approximation; (9) Directly obtain the chemical source term of the grid G ​​species by summing up the multi-step reactions of the reaction rate equation; The rate of change of species concentration is as follows: Among them, y i For species X i The mass fraction, c i is the molar concentration, M i is the molecular weight, ρ is the density; For the reaction of N species L, the reaction mechanism of the elementary step list is: where v′ i，r and v″ i，r is the stoichiometric coefficient; according to the law of mass action, species X i The production rate in the rth reaction is Then species X i The generation rate, i.e. the source term, is expressed as 3. According to the turbulent flow-chemical reaction decoupling modeling method of claim 2, step 2 specifically comprises: (1) Single molecule reaction: Cracking and isomerization reactions fall into this category: Where P is the product; the rate of change of the number of molecules of A in the grid is only related to the total number of A and has nothing to do with the distribution of A in space; solving the ordinary differential equations (ODE) of the reaction mechanism gives the average concentration use The reaction rate is calculated exactly and is independent of the flow, namely: The reaction rate is expressed by the change in product concentration, the same below; (2) Lindemann unimolecular reaction mechanism, The true reverse reaction rate constants of the first reversible reaction are k1 and k -1 The rate constant of the second one-way reaction is k2, and the quasi-steady-state species A is obtained according to the first step of fast equilibrium. * The concentration of P is expressed as: Where c M is the three-body concentration, that is, the sum of the concentrations of all species, which is a constant. The above formula degenerates into a single-molecule reaction situation, that is: (3) Bimolecular reaction: Such as: The reaction rate of a bimolecular reaction is written as: Such reactions include oxidation and free radical recombination reactions; the contribution of species concentration to the chemical source term depends on the concentration product factor c A c B ; When A and B are evenly distributed in the grid volume V, that is, at any position: but: Of the two species included in the concentration product, if one has a uniform spatial distribution of concentration within the grid G, equation (11) holds true; The reaction rate is calculated by setting the spatial distribution of the concentration of each species in the grid; when the concentration distribution of A or B in the grid is uneven, the spatial distribution of the concentration needs to be set; in the actual flow field calculation, the average concentration of the grid G ​​and the adjacent grids G-1 and G+1 is obtained. First, a one-dimensional grid in the x direction is discussed. The concentration is assumed to change linearly along the x direction. The average concentration of the left and right adjacent grids G-1 and G+1 is taken as the actual concentration of the corresponding boundary of the grid G; in order to ensure the conservation of the number of species A particles in the grid G, that is, the number of mol of species A in the grid is unchanged, where V is the volume of the grid G; two methods are used to represent the concentration distribution in the grid; first consider the linear concentration distribution case, let the one-dimensional concentration distribution, straight line L1 is the straight line connecting the left and right boundaries of the grid G, and the concentrations at the left and right boundaries are set to the average concentrations of the left neighboring grid G-1 and the right neighboring grid G+1 respectively and In order to ensure the conservation of the number of moles of species A in the grid, the straight line is translated, and the concentration at x = 0.5 is In this way, the linear distribution of the spatial concentration of species A in grid G ​​is represented by straight line L2: Consider two chemical species A and B, then: Make the same assumptions in the y and z directions, then: Calculate the reaction rate when considering the concentrations of two adjacent grids under one-dimensional conditions; integrate along the x-axis within the grid, and integrate the concentration product of the bimolecular reaction rate along the x-direction within the grid G, that is: Substituting equations (12) and (13) into equation (18), the reaction rate under one-dimensional conditions is obtained as follows: exist or When , the concentration product is reduced to the case of uniform concentration distribution or unimolecular reaction, that is: The reaction rate is calculated by averaging the integration results of the three one-dimensional concentration products, and the expression of the chemical source term is: The summation term is the summation in the x, y, and z directions, and m is the number of summation terms; (4) Trimolecular reaction and overall mechanism The trimolecular reaction contains one species M, whose concentration is c M is the trisomic concentration, c M is the sum of the concentrations of all species. Assuming that is constant, the rate equation of the trimolecular reaction degenerates into the case of a bimolecular reaction, and the rate constant is treated as a bimolecular reaction, that is, equation (21); There is also a general package mechanism for complex reaction mechanisms, and the reaction rate is as follows: The total package reaction is a form of reaction rate equation with relatively low accuracy. The reaction rate is calculated by taking the average of multiple average concentrations to simplify the calculation. Considering the linear concentration distribution assumption, the average concentration of x = 0 to 0.5 is taken in the grid G. and the average value of x = 0.5 to 1.0 The simple average of c A，x ,Right now: The reaction rate equation considering the m-dimensional concentration spatial distribution is obtained by processing similar to equation (21); The chemical source term is obtained by summing the reaction rates of all steps r in the generation and consumption of species A in the reaction mechanism, that is:

4. According to the turbulent flow-chemical reaction decoupling modeling method of claim 3, step 3 specifically comprises: Introducing non-uniform temperature distribution correction to obtain the average rate constant and will Used for chemical source term calculation; for bimolecular reactions, calculate reaction rates such as: The rate constant k is expressed by the Arrhenius equation; the activation energy E a When the average temperature is Substituting the Arrhenius equation will increase the deviation Right now: For the convenience of mathematical processing, according to Arrhenius's differential form: When the activation energy E a ≈0, the effect of T on the rate constant k value weakens. At this time, the average temperature Substituting the reaction rate equation into the calculation of the reaction rate will not lead to a large error; In E a >20kcal / mol, the effect of uneven temperature distribution on the reaction rate calculation needs to be considered, from the second-order derivative: The inflection point temperature of the k~T curve is If the activation energy E a =20 kcal·mol -1 ,T inf is 5000K, grid G The value is much lower than this temperature; it is easy to find that the average temperature is used below the inflection point temperature. Chemical reaction rates will be underestimated; in order to obtain reasonable Considering the weighted average method, let The corresponding equivalent temperature is T b ,Due to the exponential dependence of k on T, T is obtained by weighted average in the following way b ,Right now: Where T x- and T x+ are the left and right boundary temperatures of grid G, respectively, and the average temperature of grids G-1 and G+1 is taken; T x+ The weighting factor is taken as Logarithmic term Middle molecular To ensure E a =0 o'clock No divergence and weighting factor Get T b To ensure: Weighted average to obtain equivalent temperature T b , that is, the rate constant term is obtained; in addition to considering the average temperature influence of the two adjacent grids G-1 and G+1, the average temperature of the grid G ​​itself must also be considered The influence of T b and The chemical reaction rate is calculated by taking the arithmetic mean of Where c A c B It is obtained by the concentration product calculation method; a ≈0 reaction, simply take When considering a two-dimensional or three-dimensional grid, take: m = 1, 2, 3 is the dimension of the actual temperature distribution in space, which is substituted into equation (35) to obtain the chemical source term for the correction of temperature inhomogeneity.

Images