Numerical simulation method suitable for three-dimensional spherical flame

Abstract

The invention provides a numerical simulation method suitable for a three-dimensional spherical flame, and belongs to the field of electrical data processing, and the method comprises the steps: building a computational domain model according to the requirements of a computational scene, dividing a grid unit, and obtaining the node information of the required grid unit; establishing a Navier-Stokes equation set of a coupling chemical reaction and diffusion process model; estimating convection and diffusion term variable gradients in the conservation equation; calculating a time step length according to a set Kulang number, and carrying out time advancing; estimating an expected value of a velocity field, and solving a predicted velocity; the divergence of the velocity field and the correction value of the pressure field are calculated; obtaining a corrected velocity field; a self-adaptive grid means is introduced to dynamically adjust the grid resolution; performing iterative correction for multiple times until the velocity field and the pressure field converge to convergence conditions; performing spherical flame numerical simulation, and monitoring the convergence of the solution. According to the invention, the flame dynamics and cellular instability of the premixed spherical expansion flame can be reproduced.

Description

Technical Field

[0001] The invention belongs to the field of electrical data processing, and in particular relates to a numerical simulation method suitable for three-dimensional spherical flames. Background Art

[0002] Exploring the development and instability process of spherical flames, revealing the complex flow, heat transfer and chemical reaction mechanisms inside the flame, will help optimize the design of combustion systems, improve energy efficiency, promote the innovative development of new combustion technologies, deepen the theoretical understanding of flame dynamics, and provide important theoretical support and practical guidance for solving related engineering problems. At present, some inherent limitations in the experiment have not been completely resolved: first, it is difficult for the schlieren system to detect flame pre-combustion caused by pressure oscillation, temperature superposition or excessive wall temperature at the corner of the wall surface inside the constant volume combustion bomb; second, the uniformity of the temperature field in the constant volume combustion bomb cavity during heating is poor, the quartz windows at both ends of the cylindrical cavity cannot be wrapped with insulation materials, and the temperature field inside the cavity is uneven during heating. Due to the thermal convection in the closed space, the temperature field of the mixed gas in the bomb can never be absolutely uniform regardless of the heating method used, and the laminar combustion velocity of the mixed gas is a function of the initial parameters. The inhomogeneity of the temperature field will cause a large uncertainty in the measurement value of the laminar combustion velocity. Then, due to the limitation of the number and layout area of ​​the temperature probes in the projectile, it is difficult to fully reflect the temperature field distribution of the fluid medium inside the projectile, which leads to a certain deviation in the measurement of the initial temperature of the fluid medium during the experiment. Finally, the data obtained in the experiment is limited by the limited frame rate of the camera and the spatial resolution of the optical device and the influence of noise. In addition, the accurate detection of the flame front in the experiment is very complicated, which will also have a certain impact on the measurement results.

[0003] During transient flame propagation, fluid flow is complex and time-dependent. The flame propagation process is affected by multiple instabilities, and the flame dynamics behavior is violent. Therefore, it is of great significance to conduct CFD (computational fluid dynamics) numerical simulation research on the spherical flame development process of premixed gas in a constant volume combustion bomb, which helps to deeply understand the combustion process, optimize the design of the bomb body, evaluate the experimental risks, and promote the development of combustion science and technology.

[0004] In the past, computational studies on the premixed combustion process of gas fuels mainly used reaction kinetics software to perform numerical simulations of one-dimensional plane flames, or used RANS (Reynolds average method) and two-dimensional DNS (direct numerical simulation) to perform computational studies on the premixed combustion process of gas fuels, analyze the temperature field distribution, the generated free radicals and their dominant chemical reaction mechanisms, but there are the following disadvantages:

[0005] 1. Even for the one-dimensional flame of the PREMIX model in the CHEMKIN software based on the detailed chemical reaction mechanism of GRI-Mesh 3.0, the error of the methane / air laminar flame with a chemically appropriate ratio at room temperature and pressure still reaches 7.6%. =0.6) and rich fuel ( =1.3), the maximum difference is significantly larger: 40.6% and 26.4%, respectively.

[0006] 2. The research object of the two-dimensional DNS method is a plane flame rather than a spherical flame. The research using asymptotic models and spectral numerical algorithms shows that the acceleration rate of a three-dimensional spherical flame is much higher than that of a two-dimensional flame. At the same time, the wrinkles of a three-dimensional flame are much stronger and the acceleration occurs earlier. Due to the large amount of calculation, the DNS method needs to set multiple time steps to accurately reflect the changing characteristics of the flame front with time and space structure, which means that direct solution requires a lot of computing resources, so it is less used.

[0007] 3. Although the RANS method studies three-dimensional spherical flames, the obtained two-dimensional temperature field and density field distribution cloud map of the center section of the bomb cannot well reflect the cracks and cellular structures on the flame surface during the flame instability stage. Therefore, previous studies have only extracted the flame contour of the two-dimensional cross section during the flame stable propagation process. However, in the actual process, the streamlines of the flow field gradually become complex when the three-dimensional spherical flame is unstable. The dimensionality reduction process will lose a lot of information and cause inaccurate results. At the same time, the RANS method solves the time-averaged flow field, and its accuracy is limited when predicting flame dynamics.

[0008] 4. Although the numerical calculation method using the LES method and the premixed combustion model in the past restored the three-dimensional spherical flame, the calculation model error was still as high as 14.65%, and when the flame radius developed to 24 mm, a large number of cracks and slight wrinkles began to appear on the simulated flame surface, which means that the simulated flame reaction zone is narrower and the flame stability is poor. The simulated flame surface will begin to split cracks earlier, and the crack length will gradually increase with the growth of the flame. The calculation error will further increase with the growth of the flame radius. Due to the different adiabatic temperatures of the flame under different initial conditions and components, the threshold of the adaptive grid method based on the temperature gradient is different, which will lead to certain calculation errors and inaccuracies. In summary, the numerical calculation method of the prior art still needs to be further improved. Summary of the invention

[0009] In order to overcome the shortcomings of large computational complexity and uncertainty in the existing numerical simulation results of spherical flame propagation processes, the present invention provides a numerical simulation method suitable for three-dimensional spherical flames, which is suitable for the calculation of most premixed combustion scenarios, and can be used for computational fluid dynamics simulation and theoretical analysis of phenomena such as the stable propagation process of laminar spherical flames, the instability process and the turbulent spherical flame propagation process. It can reproduce the flame dynamics and cellular instability of premixed spherical expansion flames with higher accuracy and faster calculation rate, while retaining the three-dimensional spherical expansion flame front information as completely as possible, and has the characteristics of wide application range, etc.

[0010] In order to achieve the above object, the present invention adopts the following technical solution:

[0011] A numerical simulation method suitable for three-dimensional spherical flames comprises the following steps:

[0012] Step 1: construct a computational domain model and divide the grid units according to the requirements of the computational scenario to obtain the required grid unit node information, wherein the grid unit node information includes the grid unit center coordinates, the grid unit interface position coordinates, and the grid unit size;

[0013] Step 2: Initialize the physical variables in the computational domain according to the computational requirements, establish the Navier-Stokes equations for the coupled chemical reaction and diffusion process model, and set the boundary conditions according to the actual physical problem;

[0014] Step 3: Discretize the control equation space through the second-order upwind scheme, use the multidimensional linear reconstruction method combined with Taylor series to improve the accuracy, and use the bounded second-order implicit method to enhance the stability of the numerical solution. By solving the minimization problem of the non-square coefficient matrix, estimate the variable gradients of the convection and diffusion terms in the conservation equations, introduce a multidimensional limiter to suppress non-physical oscillations, and ensure that the results are stable and reliable;

[0015] Step 4: Calculate the time step according to the set Courant number and perform time advancement;

[0016] Step 5: Estimate the expected value of the velocity field based on the pressure field at the current time step; solve the predicted velocity based on the velocity field at the previous time step and the current pressure field;

[0017] Step 6: Calculate the divergence of the velocity field and the correction value of the pressure field to ensure the non-divergence condition of the velocity field when updating the velocity field;

[0018] Step 7: Based on the pressure correction, the pressure correction field obtained is applied to the predicted velocity field to obtain a corrected velocity field;

[0019] Step 8: Introduce adaptive grid means to dynamically adjust the grid resolution;

[0020] Step 9, repeat steps 5 to 8 for multiple iterations until the velocity field and the pressure field converge to the convergence condition;

[0021] Step 10: Perform numerical simulation of spherical flame and monitor the convergence of the solution.

[0022] Beneficial effects:

[0023] The present invention can highly restore the spherical flame expansion process, and study the flame dynamics and cellular instability of the spherical expansion flame in the premixed methane / air mixture, which not only helps to reveal the basic mechanism of flame formation and maintenance, but also provides a theoretical basis and technical support for the structural design and performance optimization of the low-emission gas turbine combustion chamber, and at the same time provides a theoretical and computational basis for improving energy efficiency and reducing environmental pollution. The present invention is suitable for the calculation of most premixed combustion scenarios, and can be used for computational fluid dynamics simulation and theoretical analysis of phenomena such as the stable propagation process of laminar spherical flames, the instability process, and the propagation process of turbulent spherical flames. It can reproduce the flame dynamics and cellular instability of the premixed spherical expansion flame with higher accuracy and faster calculation rate, while retaining the three-dimensional spherical expansion flame front information as completely as possible, and has the characteristics of a wide range of applications. BRIEF DESCRIPTION OF THE DRAWINGS

[0024] Figure 1 The present invention is a flow chart of a numerical simulation method suitable for three-dimensional spherical flames.

[0025] Figure 2 Schematic diagram of initial mesh division.

[0026] Figure 3 Comparison of the simulated cloud image and experimental schlieren image of methane / air premixed flame during flame propagation (T=300 K, P=0.1 MPa, =0.8).

[0027] Figure 4 The comparison of the flame radius of the methane / air premixed flame simulation results and experimental results over time (T=300 K, P=0.1 MPa, =0.8).

[0028] Figure 5 The density of methane gas / air spherical flame propagation is shown in the numerical calculation results (T=300 K, P=0.1 MPa, =0.8).

[0029] Figure 6 The temperature variation with radial distance of the numerical calculation results of methane / air spherical flame propagation (T=300K, P=0.1 MPa, =0.8).

[0030] Figure 7 For a methane / air spherical flame at T=300 K, P=1.0 MPa, =0.8 simulated cloud diagram of flame propagation process. DETAILED DESCRIPTION

[0031] In order to make the purpose, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not intended to limit the present invention. In addition, the technical features involved in the various embodiments of the present invention described below can be combined with each other as long as they do not conflict with each other.

[0032] like Figure 1 As shown, a numerical simulation method suitable for three-dimensional spherical flames of the present invention comprises the following steps:

[0033] Step 1: According to the requirements of the computing scenario, a computing domain model is constructed and grid units are divided to obtain the required grid unit node information; the grid unit node information includes the grid unit center coordinates, the grid unit interface position coordinates, and the grid unit size;

[0034] Step 2: Initialize the physical variables in the calculation domain according to the calculation requirements, establish the Navier-Stokes control equations for coupling chemical reactions and diffusion processes, and set boundary conditions according to the actual physical problems; the physical variables include temperature, pressure, velocity, and reactant components;

[0035] The Navier-Stokes governing equations are as follows:

[0036] (1)

[0037] (2)

[0038] (3)

[0039] (4)

[0040] (5)

[0041] (6)

[0042] (7)

[0043] in, is the density, t is the time, u, v, w represent the velocity components in the x, y, z directions respectively, , , represents the normal stress in the x, y, and z directions, , , , , , represents the shear stress in different directions; f x , f y , f z represents the unit mass force component in the x, y, and z directions, e is the specific internal energy, k is the thermal conductivity, q is the specific heat change rate, C S is the volume fraction of component S, D S is the diffusion coefficient, T is the fluid temperature, S S is the mass produced by the chemical reaction of component S per unit volume per unit time, is the divergence, is the velocity vector of the fluid, is the external force vector, is the stress tensor vector.

[0044] By directly calculating the motion of large-scale vortices, a mathematical model is used to describe the influence of small-scale vortices:

[0045] (8)

[0046] (9)

[0047] (10)

[0048] (11)

[0049] (12)

[0050] (13)

[0051] (14)

[0052] (15)

[0053] (16)

[0054] in, is the average value of the strain rate tensor at the decomposition scale, It is the stress tensor calculated by large eddy simulation. The superscript Te indicates a specific discretization scheme. The subscripts xy, yx, and zz mainly indicate the components of the physical quantity in different directions. It is the temporal or spatial average of the vortex tensor components, which is used to remove the influence of small-scale eddies. It represents the density of the fluid after filtering (filtering in large eddy simulation). Represents the product of fluid velocity and heat field after filtering. represents the filtering amount of the heat field, and the subscript S represents the position of the scalar field. is the diagonal part of the stress, is the trace of the velocity field, is the sub-mesh stress tensor. is the filtering amount of fluid velocity, is the large-scale velocity, Q i is the subgrid flux, and the subscripts i=1,2,3 correspond to the x, y, and z directions in the spatial coordinates. is the momentum viscosity coefficient of the subgrid turbulence, Pr SGS is the subgrid Prandtl number, the subscript SGS represents the part not captured after filtering, is the eddy viscosity coefficient, L S is the mixing length, C p is the specific heat capacity at constant pressure, is the von Karman constant, V is the grid volume, C w is the WALE constant, d is the distance to the nearest wall, is the value of the scalar field, is the filtered value of the scalar field, is the filtered value of the product of fluid density and scalar field, is the eddy diffusion coefficient.

[0055] The radius of the sphere volume of the initial fire core is r t The laminar combustion process variable c of the spark during the ignition process is calculated by the following formula:

[0056] (17)

[0057] (18)

[0058] Where r represents the spark radius, r0 represents the initial spark radius, L is the grid length scale, l t is the turbulence length scale.

[0059] After ignition, the position of the flame front moves as the combustion reaction proceeds. As the flame front moves, the reactants in the unburned zone gradually burn into combustion products. The flame front propagation is modeled by solving the transport equation of the density-weighted average reaction process variable:

[0060] (19)

[0061] Among them, c, and Sc t represent the laminar combustion process variable, the average reaction process variable and the laminar combustion gradient flow Schmidt number, respectively. is the average flow velocity vector, C p is the specific heat capacity at constant pressure, k is the diffusion coefficient, is the momentum coefficient, S c It is the source term of the reaction process.

[0062] The reaction process variable c represents the total mass fraction of the product, and the formula is:

[0063] (20)

[0064] Among them, the superscript mei is the unburned reactant, Y k is the volume fraction of reactant k, Y c is the mass fraction of the product, is the mass fraction of the product after complete adiabatic combustion, a k is a constant, the reaction process variable c of the unreacted product is 0, and the reaction process variable c of the product is 1.

[0065] Average Response Process Variable Calculate using the following formula:

[0066] (twenty one)

[0067] Among them, ρ u is the density of the unburned mixture, U t is the flame propagation speed, which is calculated by formula (13):

[0068] (twenty two)

[0069] (twenty three)

[0070] (twenty four)

[0071] (25)

[0072] (26)

[0073] Among them, U t is the flame propagation speed, A is the constant of the premixed combustion model, u' is the root mean square speed of premixed combustion, u L is the laminar flame propagation speed, is the thermal diffusion coefficient of the unburned gas, is the turbulence time scale, is the chemical reaction time scale, l t is the large eddy length scale, C W is the WALE constant, L is the grid size, is the sub-grid stress tensor, is the strain rate tensor, u x , u y represents the velocity components x and x in the x and y directions of the fluid j Represents the position of the fluid in space.

[0074] Step 3: Use the second-order upwind scheme to spatially discretize the control equations of transient simulation, use the multidimensional linear reconstruction method to reconstruct the conserved variable values ​​on both sides of the grid unit interface according to the mean value of the conserved variables in the grid unit in the computational domain, and perform Taylor series expansion on the grid-centered solution about the grid centroid to obtain higher accuracy on the grid surface. Use the bounded second-order implicit method to evolve the time variation of the flow field and combustion field, and improve the stability of the numerical solution by implicitly processing the time step.

[0075] The gradient calculation uses the Least Squares Cell-Based method to accurately estimate the gradient of the variables in the calculation of convection and diffusion terms in the flow conservation equation by solving the minimization problem of the non-square coefficient matrix in the sense of least squares. A multidimensional limiter is introduced as a gradient limiter to suppress the non-physical oscillation of the numerical solution by checking the gradient on each grid and trimming only the normal component of the gradient to the grid surface, thereby ensuring the stability and reliability of the simulation results.

[0076] The unsteady conservation equation of the scalar transport equation based on the control volume is discretized as follows:

[0077] (27)

[0078] Where V is the volume of any control volume, is the density, is the velocity vector, is the surface area vector, is a scalar, is the scalar diffusion coefficient, is the control volume source term per unit volume.

[0079] The linearized form is:

[0080] (28)

[0081] The subscript nb represents the adjacent cell, a P and a nb for and The linearization coefficient of .

[0082] The multidimensional linear reconstruction method is used to calculate the physical quantities on the unit surface. By performing a Taylor series expansion on the solution about the grid centroid centered on the grid, the second-order upwind scheme is selected to calculate the value on the grid surface:

[0083] (29)

[0084] in, The value of the physical quantity calculated on the mesh surface is: and is the value at the center of the grid and the gradient value in the upstream grid, is the displacement vector from the upstream grid center to the face center.

[0085] Any independent variable can be discretized in time as:

[0086] (30)

[0087] (31)

[0088] (32)

[0089] in, is the value at the center of the grid, t is the time, n, n-1, n-2, n +1 / 2, n-1 / 2 are different time points, and is a bounded factor of the variable at time points n+1 / 2 and n-1 / 2.

[0090] The calculation of the convection and diffusion terms in the flow conservation equation requires the use of the gradient of the variable. The cell value between grid c0 and grid c1 changes along the vector from the center of grid c0 to grid c1, which can be expressed as:

[0091] (33)

[0092] (34)

[0093] in, represents the gradient operator, represents the change, [J] is the coefficient matrix of the pure geometric function, represents the scalar field on the grid surface c0 The gradient of From mesh surface c0 to mesh surface c i The distance and Indicates that on the grid surface c i and the physical quantity at c0, is a scalar field Gradient operator, is a scalar field The amount of change.

[0094] The grid center gradient is:

[0095] (35)

[0096] (36)

[0097] (37)

[0098] in, , , is the value of the physical quantity of the grid center in the x, y, and z directions, i is the index variable, representing each adjacent node in the grid, and n is the total number of grid points adjacent to the grid center. , , are the weight coefficients in the x, y, and z directions, is the value of the physical quantity at the i-th grid point.

[0099] Numerical solution on 3D mesh elements , its gradient can be expressed as ,in, is the gradient of function u, and and Represent the rate of change of function u' in the x, y, and z directions respectively.

[0100] The gradient along each coordinate direction is calculated as:

[0101] (38)

[0102] (39)

[0103] (40)

[0104] in, Represents the amount of change, i, j, k are the coordinate indexes of the grid, indicating the position on the three-dimensional discrete grid. , , It is the differential component of the numerical solution in the x, y, and z directions. It represents the difference between the physical quantities of two adjacent grid points (i+1 and i) along the x direction when the i-th grid point, j-th and k-th coordinates are fixed. , , and are the physical quantity u at the grid points (i+1,j,k), (i,j+1,k), (i,j,k+1) and (i,j,k) respectively.

[0105] Apply the Minmod (minimum modulus) function to constrain the gradient in each direction:

[0106] (41)

[0107] (42)

[0108] (43)

[0109] (44)

[0110] in, , , It is the difference value after restriction corresponding to the x, y, and z directions, indicating the local gradient change in each direction. The two variables a and b usually represent the values ​​in adjacent grid points or discretization schemes. sign(a) is a sign function, which means returning the sign of variable a, that is, if a is positive, it returns 1; if a is negative, it returns -1; if a is 0, it returns 0. min(|a|, |b|) represents the minimum of the absolute values ​​of a and b. and The gradients or differences at the grid points (i, j, k) and (i-1, j, k) represent the rate of change of the point relative to the adjacent points. i, j, k represent the coordinates of the discrete grid points in three directions.

[0111] Reconstruct the corrected numerical solution for:

[0112] (45)

[0113] in, represents the updated value of the variable u at the grid point (i, j, k), is the original numerical solution at the grid point (i, j, k), , , is the finite difference calculation result of the gradient along the x, y, z directions. i, j, k are the coordinate indices of the grid, indicating the position on the three-dimensional discrete grid.

[0114] Step 4: Calculate the time step according to the set Courant number :

[0115] (46)

[0116] Where Co is the Courant number, S n is the fire propagation speed, is the time step, is the grid size.

[0117] Step 5: Estimate the expected value of the velocity field based on the pressure field at the current time step; solve the predicted velocity based on the velocity field at the previous time step and the current pressure field.

[0118] The predicted velocity field is:

[0119] (47)

[0120] Among them, a P is the coefficient matrix of the pressure equation, a N is the influence coefficient related to the neighboring point N when calculating the velocity field, is the grid point P at the The vector value of the velocity field predicted by the step, is the velocity prediction value vector at grid point N in the PISO algorithm, is the source term vector, P 0 is the initial pressure field, is the gradient operator.

[0121] Step 6: Calculate the divergence of the velocity field and the correction value of the pressure field to ensure the non-divergence condition of the velocity field when updating the velocity field.

[0122] Step 7: Based on the pressure correction, apply the solved pressure correction field to the predicted velocity field to obtain the corrected velocity field:

[0123] The pressure field to be sought:

[0124] (48)

[0125] The velocity field to be determined is:

[0126] (49)

[0127] in, is the pressure field coefficient, P n is the pressure field to be solved, is the gradient operator, Hb represents the influence of boundary conditions, Indicates the movement or pressure change of a fluid in a specific direction. is the time correction term in the discretization process, is the velocity field, is the velocity field to be determined, P is the pressure field, and n is the current time step.

[0128] Step 8. During the spherical flame combustion process, a strong chemical reaction takes place in the fast reaction zone, with large temperature gradients and component concentration gradients providing the driving force for flame self-sustaining. Therefore, it is necessary to ensure that there are enough grids in the flame surface. Adaptive grid means are introduced to dynamically adjust the grid resolution. First, the numerical error on the current grid is estimated. After successfully evaluating the grid metric tensor, the individual grids used for refinement and coarsening are marked according to the adaptive function created from the solution data. The grid with poor grid quality is found and modified to obtain the required tolerance. According to the adaptability standard, the grid resolution is dynamically adjusted to adapt to new computing needs using operations such as node insertion, node movement, node deletion, and edge exchange. The grid resolution is increased in areas with large errors or drastic physical changes (at the flame front), and the grid resolution is reduced in areas with small errors or slow changes (burned and unburned areas). Finally, the domain is interpolated from the old grid to the new grid.

[0129] The grid metric is a tensor that provides an estimate of the a posteriori error measure of the field. The expected error metric of the field is first estimated using a grid adaptation method:

[0130] (50)

[0131] (51)

[0132] in, is the error index, V is the grid unit volume, r is the gradient volume weight, is the Euclidean norm of the process variable gradient, c represents the laminar combustion process variable, and x, y, z represent different directions in three-dimensional space.

[0133] Step 9, repeat steps 5 to 8 for multiple iterations until the velocity field and pressure field converge to sufficient accuracy to meet the set convergence conditions. In each iteration cycle, the velocity and pressure are re-predicted, and then pressure correction and velocity correction are performed. The quality of the grid is measured through function iteration, the worst element in the grid is located and modified to obtain the required tolerance, and the numerical solution is re-performed after adjusting the grid according to the flow field characteristics.

[0134] Step 10: Perform numerical simulation of spherical flame and monitor the convergence of the solution. In each time step, the pressure-velocity coupling solver is used to iteratively solve the fluid dynamics equations until the convergence criterion is reached.

[0135] Example:

[0136] Taking the numerical calculation of the methane / air premixed three-dimensional spherical laminar flame propagation process under the conditions of initial temperature of 300 K, initial pressure of 1 atm and 10 atm, and equivalence ratio of 0.8 as an example, the present invention is described in full and detail. The present invention is not limited to laminar flow scenarios, but is applicable to premixed spherical flames of all flow states. The embodiment of the present invention is mainly divided into the following steps: step (1) constructing a computational domain model; step (2) computational domain grid division, establishment and description of initial conditions and boundary conditions; step (3) control equation discretization; step (4) setting the time step; step (5) predicting the velocity field according to the pressure field; step (6) calculating the velocity field divergence and correcting the pressure field; step (7) applying the pressure correction field to update the velocity field; step (8) adaptive grid adjustment; step (9) iterative correction and grid adjustment; step (10) flame numerical simulation and convergence monitoring.

[0137] Step 1: Build the computational domain model:

[0138] According to the requirements of the calculation scenario, the laminar flame propagation test process is relatively simple. After all fuels and air are introduced into the cavity and fully mixed, the ignition operation can be carried out and the high-speed camera and oscilloscope can be triggered to record the flame propagation image. In order to save computing resources and considering that spherical flame propagation has a certain symmetry, the inlet and exhaust systems and other devices that are not directly related to the premixed combustion of gas fuel are omitted, and the flame propagation area is simplified to a spherical calculation domain with a diameter of 160 mm, see Figure 2 In order to better reflect the boundary characteristics of the flame shape and the controllability of the grid node unit distribution, the fluid domain grid is divided into tetrahedral grid units. After the grid independence test, the final Spacing Ralative of the grid is set to 1.2, and the initial minimum grid volume is 1.0 e -11 m 3 , the initial number of grids is 200,000.

[0139] Step 2: Initialize the physical variables including temperature, pressure, reactant components, etc. in the calculation domain according to the calculation requirements, and set the boundary conditions according to the actual physical problem. The initial temperature is set to 300 K, the initial pressure is set to 1 atm, the laminar combustion velocity is obtained according to experimental measurements, and the physical parameters are calculated by Gaseq software. Since the outer wall of the constant volume combustion bomb is usually wrapped with insulation materials with good insulation performance and the combustion process time is extremely short during the experiment, the outer wall of the model is set to an adiabatic boundary condition of 298 K at ambient temperature during the numerical simulation. Since the generation of sparks is an instantaneous and tiny phenomenon, it is difficult to fully describe it in a conventional calculation grid. Therefore, the electric spark model is used to establish a fine grid inside the calculation grid to simulate and track the generation, diffusion and combustion process of sparks. At the same time, it can ensure that the initial fire core is a standard sphere, set the ignition energy in the model to 0 mJ and ensure that the electric spark diffusion is a laminar process, the initial radius of the spark is 1 mm, and the minimum grid volume is set to 10 -16 m 3 The iteration frequency is set to 1 time step.

[0140] Step 3. Discretization of control equations: Use the second-order upwind scheme to spatially discretize the control equations of transient simulation. Use the multidimensional linear reconstruction method to reconstruct the values ​​of the conserved variables on both sides of the grid unit interface based on the mean values ​​of the conserved variables in the grid units in the computational domain. Perform a Taylor series expansion on the grid-centered solution about the grid centroid to obtain higher accuracy on the grid surface. Use a bounded second-order implicit method to evolve the time variation of the flow field and combustion field, and improve the stability of the numerical solution by implicitly processing the time step. Use the finite volume method to numerically discretize the density equation, velocity equation, component equation, temperature equation, and pressure equation on all grids;

[0141] Step 4: Set the time step to 10 -6 s.

[0142] Step 5: Estimate the flow state by predicting the velocity field.

[0143] Step 6: Based on step 5, pressure correction is performed to ensure that the flow field meets the mass conservation condition, thereby providing accurate flow field data for subsequent calculations.

[0144] Step 7. After completing the pressure correction calculation, the pressure increment obtained by the pressure field correction will be used to correct the predicted velocity field to ensure that the fluid flow satisfies the mass conservation (i.e., the divergence is zero).

[0145] Step 8: By introducing adaptive mesh technology, the mesh resolution can be dynamically adjusted during the numerical simulation process, so that the details of the flame front area can be accurately captured while reducing the amount of calculation in other areas. The process variable gradient adaptive method is used to refine the mesh of the flame front surface, and the maximum refinement level is set to 3.

[0146] Step 9: Iterative correction and grid adjustment. Repeat steps 5 to 8 for multiple iterative corrections until the velocity field and pressure field converge to a sufficient accuracy to meet the set convergence conditions.

[0147] Step 10: Flame numerical simulation and convergence monitoring. Perform spherical flame numerical simulation and use the pressure-velocity coupling solver to iteratively solve the fluid dynamics equations in each time step until the solution meets the predetermined convergence criteria. At the same time, monitor the convergence of the solution to ensure that the numerical solution is stable and accurately captures the flame process.

[0148] Figure 3 Comparison between the simulated cloud map and the experimental schlieren map of the methane / air spherical flame evolution over time under the conditions of initial temperature of 300 K, initial pressure of 1 atm, and equivalence ratio of 0.8, where T is the initial temperature, P is the initial pressure, is the equivalence ratio, R a is the area equivalent radius of the spherical flame. Figure 3 It can be seen that the flame fronts of the experiment and simulation remain smooth and have no significant difference during the flame propagation process, and the flame morphology has been greatly improved compared with the existing technology. The combustible gas mixture starts to burn at the center of the bomb and quickly diffuses to the surrounding area. As the combustion process proceeds, the density at the center of the bomb gradually decreases and the corresponding area continues to expand, forming an approximately circular area with an increasing radius and wrinkles. That is, the space inside the bomb is divided into two parts by the flame surface: a high-density area and a low-density area. The density of the unburned area is larger, which is the density of the initial gas mixture. The density of the burned area is smaller. This is because the combustion releases a large amount of heat, which causes the temperature of the burned area to rise rapidly. Due to thermal expansion, the density of this area decreases and gradually diffuses to the surrounding unburned areas. Therefore, the burned area has a lower density than the unburned area, and the low-density area gradually expands with the development of the flame. Figure 4 As shown. Since the ignition position is set at the center of the constant volume combustion bomb, that is, at the initial moment of combustion, the center of the constant volume combustion bomb is the burned area, and the gas temperature is the highest. These high-temperature gases heat the surrounding low-temperature unburned gases through convection diffusion. When the temperature of the unburned gas reaches the ignition temperature, it begins to burn. The flame spreads to the surroundings in an approximately circular shape, and the temperature gradually decreases from the high-temperature area to the outside. Figure 5 The variation of flame radius with time under the corresponding working conditions is shown in Figure 6 As shown, it can be seen that during the stable flame propagation period, the simulated flame propagation speed is slightly lower than the experimental value. The results calculated by the method of the present invention are in good agreement with the experimental results, and the error is controlled within 10%. Figure 7 The simulated cloud diagram of the time evolution of a methane / air spherical flame with an initial pressure of 10 atm is compared with Figure 3 and Figure 7It can be seen that the time when cracks and protrusions appear during flame propagation becomes earlier with the increase of initial pressure. In the initial stage of flame development, the cracks and cellular structures of the working condition with higher initial pressure will increase significantly, and the flame surface will change from smooth to wrinkled. The development of the flame will aggravate the degree of wrinkling. The cracks will spread more rapidly to the entire surface of the fire core with the development of the flame. The flame is greatly affected by fluid mechanics instability. The method of the present invention can capture the detailed characteristics of the flow field during the combustion process while ensuring the reliability of the calculation results.

Claims

1. A numerical simulation method suitable for three-dimensional spherical flames, characterized in that: The steps include: Step 1: construct a computational domain model and divide the grid units according to the requirements of the computational scenario to obtain the required grid unit node information, wherein the grid unit node information includes the grid unit center coordinates, the grid unit interface position coordinates, and the grid unit size; Step 2: Initialize the physical variables in the computational domain according to the computational requirements, establish the Navier-Stokes equations for the coupled chemical reaction and diffusion process model, and set the boundary conditions according to the actual physical problem; Step 3: Discretize the control equation space through the second-order upwind scheme, use the multidimensional linear reconstruction method combined with Taylor series to improve the accuracy, use the bounded second-order implicit method to enhance the stability of the numerical solution, estimate the variable gradients of the convection and diffusion terms in the conservation equations by solving the minimization problem of the non-square coefficient matrix, introduce a multidimensional limiter to suppress non-physical oscillations, and ensure that the results are stable and reliable; Step 4: Calculate the time step according to the set Courant number and perform time advancement; Step 5: Estimate the expected value of the velocity field based on the pressure field at the current time step; solve the predicted velocity based on the velocity field at the previous time step and the current pressure field; Step 6: Calculate the divergence of the velocity field and the correction value of the pressure field to ensure the non-divergence condition of the velocity field when updating the velocity field; Step 7: Based on the pressure correction, the pressure correction field obtained is applied to the predicted velocity field to obtain a corrected velocity field; Step 8: Introduce adaptive grid means to dynamically adjust the grid resolution; Step 9, repeat steps 5 to 8 for multiple iterations until the velocity field and the pressure field converge to the convergence condition; Step 10: Perform numerical simulation of spherical flame and monitor the convergence of the solution.

2. A numerical simulation method for three-dimensional spherical flames according to claim 1, characterized in that: The step 2 comprises: By directly calculating the motion of large-scale vortices, a mathematical model is used to describe the influence of small-scale vortices: (8) (9) (11) (12) (13) (14) (15) (16) in, is the average value of the strain rate tensor at the decomposition scale, It is the stress tensor calculated by large eddy simulation. The superscript Te indicates a specific discretization scheme. The subscripts xy, yx, and zz mainly indicate the components of the physical quantity in different directions. It is the temporal or spatial average of the vortex tensor components, which is used to remove the influence of small-scale eddies. It represents the density of the fluid after filtering (filtering in large eddy simulation). Represents the product of fluid velocity and heat field after filtering. represents the filtering amount of the heat field, and the subscript S represents the position of the scalar field. is the diagonal part of the stress, is the trace of the velocity field, is the sub-mesh stress tensor. is the filtering amount of fluid velocity, is the large-scale velocity, Q i is the subgrid flux, and the subscripts i=1,2,3 correspond to the x, y, and z directions in the spatial coordinates. is the momentum viscosity coefficient of the subgrid turbulence, Pr SGS is the subgrid Prandtl number, the subscript SGS represents the part not captured after filtering, is the eddy viscosity coefficient, L S is the mixing length, C p is the specific heat capacity at constant pressure, is the von Karman constant, V is the grid volume, C w is the WALE constant, d is the distance to the nearest wall, is the value of the scalar field, is the eddy diffusion coefficient.

3. A numerical simulation method suitable for three-dimensional spherical flames according to claim 2, characterized in that: The step 3 comprises: The multidimensional linear reconstruction method is used to calculate the physical quantities on the unit surface. The second-order upwind calculation grid surface is selected by performing a Taylor series expansion of the solution about the grid centroid centered on the grid: Numerical solution on 3D mesh elements The gradient of ; represents the velocity component in the x direction; Compute the gradient along each coordinate direction and reconstruct the corrected numerical solution: (45) in, represents the updated value of the variable u at the grid point (i, j, k), is the original numerical solution at the grid point (i, j, k), , , is the finite difference calculation result of the gradient along the x, y, and z directions; i, j, and k are the coordinate indices of the grid, indicating the position on the three-dimensional discrete grid.

4. A numerical simulation method for three-dimensional spherical flames according to claim 1, characterized in that: The step 4 includes: determining the global time step and performing time advancement.

5. A numerical simulation method for three-dimensional spherical flames according to claim 1, characterized in that: The step 8 comprises: Expected error measure of the estimated field: (50) (51) in, is the error index, V is the grid unit volume, r is the gradient volume weight, is the Euclidean norm of the process variable gradient, c represents the laminar combustion process variable, and x, y, z represent different directions in three-dimensional space.

Images