CFD (Computational Fluid Dynamics) method for near-continuous hypersonic thermochemical non-equilibrium flow

Abstract

The invention relates to a CFD (Computational Fluid Dynamics) method for nearly continuous hypersonic high-temperature thermochemical unbalanced flow numerical simulation, which comprises the following steps of: 1, introducing a Park double-temperature model into an NS equation to describe a thermochemical unbalanced effect at high temperature, and closing a shear stress tensor, a heat flow vector and a mass diffusion flux by adopting an equivalent constitutive model to obtain a closed shear stress model; capturing speed slippage and temperature jump phenomena on the wall surface by adopting corresponding slippage boundary conditions to obtain a macroscopic control equation and boundary conditions suitable for near-continuous hypersonic high-temperature thermochemical non-equilibrium flow; 2, discretely solving the macroscopic control equation and the boundary conditions given in the step 1 by adopting a finite volume method; and step 3, performing post-processing on a solution result obtained in the step 2 to obtain a streaming flow field, a wall surface aerodynamic force and a wall surface heat flow of the calculated appearance.

Description

Technical Field

[0001] This invention belongs to the interdisciplinary field of rarefied gas dynamics and high-temperature gas dynamics, and particularly relates to a computational fluid dynamics (CFD) method suitable for simulating near-continuous hypersonic high-temperature thermochemical nonequilibrium flows. Background Technology

[0002] When near-space hypersonic vehicles traverse the atmosphere, the flow field exhibits significant rarefied gas effects due to the thin incoming gas and strong boundary layer shear. Simultaneously, intense aerodynamic heating causes vibrational excitation, dissociation, and even ionization of nitrogen and oxygen molecules in the air, exposing them to high-temperature real gas effects. These two effects are coupled, rendering traditional CFD methods inaccurate in predicting the aerodynamic and thermal environments of the vehicle surface. To accurately assess the aerodynamic and thermal environment of near-space hypersonic vehicles, computational models need to simultaneously consider both rarefied gas effects and high-temperature thermochemical nonequilibrium effects.

[0003] Generally, rarefied thermochemical nonequilibrium flows can be described by the generalized Boltzmann equation (GBE). However, due to the extremely high degrees of freedom and complex integral collision terms of the GBE, direct solutions are very difficult. Direct simulation Monte Carlo (DSMC) methods, as a probabilistic particle simulation method, can accurately handle physicochemical nonequilibrium processes at the molecular mean free path scale and have been successfully applied to the simulation of hypersonic rarefied thermochemical nonequilibrium flows. However, in the near-continuous flow region, DSMC is computationally very expensive. In continuous flow, thermochemical nonequilibrium flows are usually described by the extended Navier–Stokes (NS) equations based on the Park two-temperature model and solved using CFD methods. Combining the Park two-temperature model with slip boundary conditions remains the main approach for simulating near-continuous thermochemical nonequilibrium flows. However, as local rarefaction effects intensify, the linear constitutive relations upon which the NS equations based on the continuous medium assumption gradually fail, and the prediction errors of traditional CFD for aircraft drag and heat flux also increase. Currently, for the simulation of hypersonic near-continuous flow at high temperature thermochemical non-equilibrium, traditional CFD methods are no longer accurate, and theoretical methods are computationally extremely expensive. There is an urgent need to develop general simulation methods that balance accuracy and computational efficiency to achieve complex engineering applications. Summary of the Invention

[0004] The purpose of this invention is to provide a CFD method that can be used for numerical simulation of near-continuous hypersonic high-temperature thermochemical nonequilibrium flows, while meeting the computational requirements for engineering calculations. The technical solution of this invention is as follows:

[0005] A CFD method for numerical simulation of near-continuous hypersonic high-temperature thermochemical nonequilibrium flows includes the following steps:

[0006] Step 1: Introduce the Park dual-temperature model into the NS equation to describe the thermochemical non-equilibrium effect at high temperature. Use an equivalent constitutive model to close the shear stress tensor, heat flux vector and mass diffusion flux. Apply corresponding slip boundary conditions on the wall to capture velocity slip and temperature jump phenomena. Obtain macroscopic governing equations and boundary conditions applicable to near-continuous hypersonic high-temperature thermochemical non-equilibrium flow.

[0007] Step 2: Discretize and solve the macroscopic governing equations and boundary conditions given in Step 1 using the finite volume method;

[0008] Step 3: Post-process the solution results from Step 2 to obtain the flow field around the calculated shape, the aerodynamic forces on the wall, and the heat flux on the wall.

[0009] Furthermore, the method for step one is as follows:

[0010] 1) The Navier-Stokes equations of the three-dimensional compressible multi-component Park two-temperature model describing thermochemical nonequilibrium flow in Cartesian coordinates can be written in the following conservation form.

[0011]

[0012] Where U is the conserved variable; t is time; (x,y,z) are Cartesian coordinates; F x F y and F z The inviscid fluxes in the x, y, and z directions, respectively; F v,x F v,y and F v,z denoted as viscous fluxes in the x, y, and z directions, respectively; S is the thermochemical nonequilibrium source term, where,

[0013] U = [ρ s ,ρu,ρv,ρw,ρE,ρe ve ] T ,

[0014] F x =[ρ s u,ρu 2 +p,ρuv,ρuw,(ρE+p)u,ρe ve u] T ,

[0015] F y =[ρ s v,ρvu,ρv 2 +p,ρvw,(ρE+p)v,ρe ve v] T ,

[0016] F z =[ρ sw,ρwu,ρwv,ρw 2 +p,(ρE+p)w,ρe ve w] T ,

[0017]

[0018] In the formula, ρ s It is the density of component s, where s ranges from 1 to N. s Values, N s ρ is the total number of components in the system; p is the total density of the gas mixture; (u, v, w) are the velocity components in the x, y, and z directions, respectively; and E is the total energy of the gas mixture. e is the internal energy per unit mass of the gas mixture, which is composed of the translational-rotational energy of the gas mixture. tr Vibration-electron energy e ve and the zero-point formation enthalpy Δh of the mixed gas f,0 Composition; h s e is the absolute enthalpy per unit mass of component s; ve,s Let the vibrational-electronic energy of component s be denoted as , and the electronic energy e be neglected. el,s The influence of, and on the vibrational energy e of monatomic gases. vib,s It is zero; J s τ is the mass diffusion flux vector of component s, τ is the shear stress tensor, and q is the mass diffusion flux vector of component s. tr and q ve These are the translational-rotational and vibrational-electron heat flow vectors, respectively. It is the mass formation rate of component s. The source term in the vibrational energy equation for a mixed gas; the mass generation rate of component s. It can be written as:

[0019]

[0020] Among them, M s υ is the molar mass of component s; sr and υ s ′ r These are the stoichiometric coefficients of the reactants and products in the r-th reaction, respectively; [X] s ] r k is the molar concentration of the r-th reactant component s; fr and k br These are the forward and reverse reaction rate coefficients for the r-th reaction, respectively; the subscript "r" indicates the r-th chemical reaction, "f" indicates the forward reaction, and "b" indicates the reverse reaction; the forward reaction rate coefficient k of the r-th reaction is... f The calculation is performed using the modified Arrhenius equation:

[0021]

[0022] In the formula, C fr η is the exponential factor of the positive reaction of the r-th reaction. fr θ is the forward reaction temperature exponent of the r-th reaction. dr T is the characteristic temperature of the dissociation of reactant particles in the r-th reaction. tr Translational-rotational temperature; Z r It is the vibration-dissociation coupling nonequilibrium factor, used to characterize the effect of vibration nonequilibrium on the dissociation rate, expressed as:

[0023]

[0024] In the formula, T a The geometric mean temperature is defined as follows: Among them, T ve For vibrational-electron temperature, k br The reverse reaction rate coefficient;

[0025] 2) Introduce a linear constitutive correction model to close the above conservation equations (1) and (2), and give the shear stress tensor τ and the translational-rotational heat flux vector q. tr Vibrational-electron heat flow vector q ve And the mass diffusion flux vector J of component s s Represented in component form:

[0026]

[0027] Among them, u i These are the components of the velocity vector (u,v,w); x i These are the components of the spatial coordinates (x, y, z); the subscripts i, j, k can all take the values ​​1, 2, 3, (x1, x2, x3) = (x, y, z) are Cartesian coordinates, (u1, u2, u3) = (u, v, w) are velocities, and δ is the velocity. ij For Kronecker notation; y s μ is the mass fraction of component s; e , κ e,tran , κ e,rot and κ e,vib These are the equivalent viscosity coefficient, equivalent translational, rotational, and vibrational thermal conductivity coefficients of the gas mixture, and D. e,s The effective diffusion coefficient of component s is given by the DSMC method for high-speed thermochemical non-equilibrium shear flow analysis, and their calculation formulas are as follows:

[0028] μ e =A μ ·μ (8)

[0029] κ e,tran =Aκ ·κ tran ,κ e,rot =A μ ·κ rot ,κ e,vib =A μ ·κ vib (9)

[0030] D e,s =(κ) e,tran +κ e,rot )Le / (ρc p,tr (10)

[0031] Where μ is the viscosity coefficient of the gas mixture; κ tran κ rot and κ vib , respectively, represent the translational, rotational, and vibrational thermal conductivity coefficients of the gas mixture; Le is the Lewis number; c p,tr A represents the translational-rotational isobaric specific heat capacity of the gas mixture; μ and A κ These are the correction functions for viscous and translational thermal conductivity coefficients, expressed as functions of the local shear nonequilibrium parameter Zh;

[0032] 3) The generalized Maxwell-Smoluchowski model is used to describe wall velocity slip and temperature jump phenomena on solid walls. It is compatible with the above-mentioned linear constitutive correction model and Park's two-temperature thermochemical nonequilibrium model. For three-dimensional flow, its specific form is as follows:

[0033]

[0034] Among them, u slip It is the tangential sliding velocity vector of the wall, T w For wall temperature; σ u , σ T and σ T,ve These are the tangential momentum adaptation coefficient, translational-rotational, and vibrational-electron energy regulation coefficients, respectively; γ is the specific heat ratio of the gas mixture; c v,tr and c v,vib These are the translational-rotational and vibrational isochoric specific heat capacities of the gas mixture, respectively; τ w This represents the tangential shear stress vector of the wall.

[0035] 4) The above conservation equation (1) and the linear constitutive correction model (7) combined with the wall slip boundary condition (14) constitute the macroscopic control equation and boundary conditions describing near-continuous hypersonic high-temperature thermochemical non-equilibrium flow.

[0036] Furthermore, in sub-step 1) of step one,

[0037] The source term in the vibrational energy equation of the mixed gas in formula (2) Represented as:

[0038]

[0039] In the formula, the first term represents the energy exchange between the translational-rotational modes and vibrational modes of the molecular component calculated according to the Landau–Teller equation; the second term represents the vibrational energy change caused by the gas chemical reaction; τ vib,s It is the molar average relaxation time between the translational-rotational mode and the vibrational mode of component s;

[0040] Furthermore, in sub-step 1) of step one,

[0041] In formula (5), the constant a is taken as 0.5 or 0.7; the reverse reaction rate coefficient k br Or through the equilibrium constant K eq,r According to k br =K eq,r / k fr Formula calculation or based on k br The fitting equation was obtained.

[0042] Furthermore, in sub-step 2) of step one,

[0043] A μ and A κ It is expressed as a function of the local shear nonequilibrium parameter Zh as follows:

[0044]

[0045] In the formula, the expression for Zh is:

[0046]

[0047] in It is the most probable velocity of local molecular thermal motion, and R is the gas mixture constant. For gas rate The derivative along the streamline normal n; λ e The equivalent molecular mean free path,

[0048] For the equivalent viscosity coefficient (μ) given by equations (8)-(10) in step one e ) and equivalent thermal conductivity (κ) e,tran , κ e,rot and κ e,vib ) and the equivalent diffusion coefficient (D) of component s e,s The calculations all rely on the local shear nonequilibrium parameter Zh in the flow field. The calculations of equations (8) and (13) together form the following implicit equation.

[0049]

[0050] Where λ is the mean free path of the molecule. Zh NS Calculated from the flow field variables and their gradients;

[0051] In step two, the calculation method of equation (15) is as follows: Zh is obtained by Newton's iteration method, and the obtained Zh is substituted into equations (11)-(12) to obtain the viscosity correction function A. μ Parallel dynamic heat conduction correction function A κ The equivalent viscosity coefficient, equivalent translational, rotational, and vibrational thermal conductivity coefficients and the equivalent diffusion coefficient of component s are calculated using equations (8)-(10).

[0052] Furthermore, in step three, the flow field around the exterior includes gas density, velocity, temperature, and component mass fraction, obtained from simulation results; the wall aerodynamic forces include wall pressure p. w and wall tangential shear stress, wall pressure p w The simulation results yield the following formulas for calculating the tangential shear stress and heat flow at the wall:

[0053] τ w,i =τ ij n w,j -(τ jk n w,k n w,j )n w,i (16)

[0054]

[0055] In the formula, n w,i The component of the unit normal vector of the wall in the i direction; the right-hand side of the calculation of wall heat flow in equation (17) contains three parts: the first term is the sliding friction heat caused by the work done by the surface force; the second term is the conduction heat of different energy modes of translation, rotation and vibration; and the third term is the component diffusion heat.

[0056] The beneficial effects of this invention compared to the prior art are as follows:

[0057] (1) This invention establishes a CFD method suitable for simulating near-continuous hypersonic thermochemical non-equilibrium flow, which is more suitable for complex engineering problems in near-space flight than existing methods.

[0058] (2) This invention takes into account the influence of the local high temperature and thinning coupling effect of the boundary layer, which significantly improves the prediction ability of near-continuous flow compared with the traditional high temperature thermochemical non-equilibrium flow CFD method.

[0059] (3) This invention is easy to implement and can be conveniently embedded into any existing CFD solver that considers high-temperature thermochemical non-equilibrium effects. Attached Figure Description

[0060] Figure 1 This is a flowchart illustrating the implementation steps of the present invention.

[0061] Figure 2 This invention provides a velocity profile of nitrogen gas flowing around a nitrogen plate at x = 0.5 m with a nitrogen Mach 10 and Knudsen number of 0.01, calculated by this invention, compared with other methods.

[0062] Figure 3 The wall friction coefficient calculated by this invention for air around a blunt cone with a Mach 25 and Knudsen number of 0.16 is compared with other methods.

[0063] Figure 4 This invention compares the wall heat flux coefficients of air around a 5-component HTV2 at Mach 25 and 85 km with those calculated by other methods. Detailed Implementation

[0064] The present invention will now be described in detail with reference to specific embodiments and accompanying drawings. In the following description, specific details are set forth for purposes of explanation and not limitation, in order to aid in a comprehensive understanding of the invention. However, it will be apparent to those skilled in the art that the invention may be practiced in other embodiments departing from these specific details.

[0065] It should be noted that, in order to avoid obscuring the invention with unnecessary details, only the device structure and / or processing steps closely related to the solution according to the invention are shown in the accompanying drawings, while other details that are not closely related to the invention are omitted.

[0066] This invention provides a CFD method for simulating near-continuous hypersonic thermochemical nonequilibrium flows. This method overcomes the shortcomings of existing technologies, primarily addressing the failure of existing high-temperature thermochemical nonequilibrium CFD methods for flows with localized rarefied gas effects. It provides a linear constitutive correction model that considers localized high-temperature rarefied coupling effects within the boundary layer, along with appropriate wall-slip boundary conditions, achieving stable and efficient solutions. The invention is implemented as follows: Figure 1 As shown, it includes the following steps:

[0067] Step one involves introducing the Park two-temperature model into the traditional Navier-Stokes equations to describe the thermochemical nonequilibrium effects at high temperatures. An equivalent constitutive model is used to close the shear stress tensor, heat flux vector, and mass diffusion flux. Corresponding slip boundary conditions are applied at the wall surface to capture velocity slip and temperature jump phenomena. Through these steps, macroscopic governing equations and boundary conditions applicable to near-continuous hypersonic high-temperature thermochemical nonequilibrium flow are obtained, in the following form:

[0068] 1) The three-dimensional compressible multi-component two-temperature Navier-Stokes equations describing thermochemical nonequilibrium flow in rectangular coordinates can be written in the following conservation form.

[0069]

[0070] Where U is the conserved variable; t is time; (x,y,z) are Cartesian coordinates; F x F y and F z The inviscid fluxes in the x, y, and z directions, respectively; F v,x F v,y and F v,z Let be the viscous fluxes in the x, y, and z directions, respectively; and S be the thermochemical nonequilibrium source term. Its specific expression is as follows:

[0071] U = [ρ s ,ρu,ρv,ρw,ρE,ρe ve ] T ,

[0072] F x =[ρ s u,ρu 2 +p,ρuv,ρuw,(ρE+p)u,ρe ve u] T ,

[0073] F y =[ρ s v,ρvu,ρv 2 +p,ρvw,(ρE+p)v,ρe ve v] T ,

[0074] F z =[ρ s w,ρwu,ρwv,ρw 2 +p,(ρE+p)w,ρe ve w] T ,

[0075]

[0076]

[0077] In the formula, the superscript 'T' indicates vector transpose, ρ s It is the density of component s, where s ranges from 1 to N. s Values, N s ρ is the total number of components in the system; p is the total density of the gas mixture; (u, v, w) are the velocity components in the x, y, and z directions, respectively; and E is the total energy of the gas mixture. e is the internal energy per unit mass of the gas mixture, which is composed of the translational-rotational energy of the gas mixture. tr Vibration-electron energy e ve and the zero-point formation enthalpy Δh of the mixed gas f,0 Composition; h s e is the absolute enthalpy per unit mass of component s; ve,s为组分s Vibrational-electron energy, neglecting electron energy e el,s The influence of, and on the vibrational energy e of monatomic gases. vib,s It is zero. J s τ is the mass diffusion flux vector of component s, τ is the shear stress tensor, and q is the mass diffusion flux vector of component s. tr and q ve These are the translational-rotational and vibrational-electron heat flow vectors, respectively. It is the mass formation rate of component s. This is the source term in the vibrational energy equation for a mixed gas. The mass formation rate of component s. It can be written as:

[0078]

[0079] Among them, M s υ is the molar mass of component s; sr and υ′ sr , respectively, are the stoichiometric coefficients of the reactants and products in the r-th reaction. [X] s ] r k is the molar concentration of the r-th reactant component s; fr and k br These are the forward and reverse reaction rate coefficients for the r-th reaction, respectively. The subscript "r" indicates the r-th chemical reaction, "f" represents the forward reaction, and "b" represents the reverse reaction. The forward reaction rate coefficient k for the r-th reaction is... f The calculation is performed using the modified Arrhenius equation:

[0080]

[0081] In the formula, C fr η is the exponential factor of the positive reaction of the r-th reaction. fr θ is the forward reaction temperature exponent of the r-th reaction. dr T is the characteristic temperature of the dissociation of reactant particles in the r-th reaction.tr The temperature is the translational-rotational temperature. Z r It is the vibration-dissociation coupling nonequilibrium factor, used to characterize the effect of vibration nonequilibrium on the dissociation rate, and can be expressed as:

[0082]

[0083] In the formula, T a The geometric mean temperature is defined as follows: Among them, T ve The vibrational-electron temperature is given by the constant α, which is typically taken as 0.5 or 0.7. The reverse reaction rate coefficient k... br Then it can be determined by the equilibrium constant K. eq,r Obtain (k) br =K eq,r / k fr ) or based on k br The fitted equation is obtained. The source term in the vibrational energy equation of the mixed gas in equation (2) is... It can be represented as:

[0084]

[0085] In the formula, the first term represents the energy exchange between the translational-rotational modes and vibrational modes of the molecular component calculated according to the Landau–Teller equation; the second term represents the vibrational energy change caused by the gas chemical reaction. τ vib,s It is the molar average relaxation time between the translational-rotational mode and the vibrational mode of component s.

[0086] 2) Introduce a linear constitutive correction model to close the above conservation equations (1) and (2), and give the shear stress tensor τ and the translational-rotational heat flux vector q. tr Vibrational-electron heat flow vector q ve And the mass diffusion flux vector J of component s s Represent it in terms of components:

[0087]

[0088] Among them, u i These are the components of the velocity vector (u,v,w); x i These are the components of the spatial coordinates (x, y, z); the subscripts i, j, k can all take the values ​​1, 2, 3, (x1, x2, x3) = (x, y, z) are Cartesian coordinates, (u1, u2, u3) = (u, v, w) are velocities, and δ is the velocity. ij For Kronecker notation (if i = j, δ) ij =1; if i≠j, δ ij =0). y s μ represents the mass fraction of component s.e κ e,tran κ e,rot and κ e,vib These are the equivalent viscosity coefficient, equivalent translational, rotational, and vibrational thermal conductivity coefficients of the gas mixture, and D. e,s The effective diffusion coefficient of component s can be obtained by analyzing high-speed thermochemical non-equilibrium shear flow using the DSMC method, and the final calculation formula is as follows:

[0089] μ e =A μ ·μ (8)

[0090] κ e,tran =A κ ·κ tran , κ e,rot =A μ ·κ rot , κ e,vib =A μ ·κ vib (9)

[0091] D e,s =(κ) e,tran +κ e,rot )Le / (ρc p,tr (10)

[0092] Where μ is the viscosity coefficient of the gas mixture; κ tran κ rot and κ vib , respectively, are the translational, rotational, and vibrational thermal conductivity coefficients of the gas mixture; Le is the Lewis number, generally considered to be Le = 1.4; c p,tr A represents the translational-rotational isobaric specific heat capacity of the gas mixture; μ and A κ These are the correction functions for viscous and translational thermal conductivity coefficients, which can be expressed as functions of the local shear nonequilibrium parameter Zh:

[0093]

[0094] In the formula, the expression for Zh is:

[0095]

[0096] in It is the most probable velocity of local molecular thermal motion, and R is the gas mixture constant. For gas rate The derivative of λ along the streamline normal n. e The equivalent molecular mean free path,

[0097] 3) The generalized Maxwell-Smoluchowski model is used to describe wall velocity slip and temperature jump phenomena on solid walls. It is compatible with the above-mentioned linear constitutive correction model and two-temperature thermochemical nonequilibrium model. For three-dimensional flow, its specific form is as follows:

[0098]

[0099] Among them, u slip It is the tangential sliding velocity vector of the wall, T w For wall temperature; σ u , σ T and σ T,ve These are the tangential momentum adaptation coefficient, translational-rotational, and vibrational-electron energy regulation coefficients, respectively; γ is the specific heat ratio of the gas mixture; c v,tr and c v,vib These are the translational-rotational and vibrational isochoric specific heat capacities of the gas mixture, respectively; τ w This represents the tangential shear stress vector of the wall surface.

[0100] 4) The above conservation equation (1) and the linear constitutive correction model (7) combined with the wall slip boundary condition (14) constitute the macroscopic control equation and boundary conditions describing near-continuous hypersonic high-temperature thermochemical non-equilibrium flow, also known as the NS model (DiNS) based on DSMC data improvement.

[0101] Step 2: The macroscopic control equations (1)-(14) from Step 1 are discretized and solved using the conventional finite volume method. Specifically, the inviscid flux is reconstructed using the third-order MUSCL (Monotonic Upstream-centered Schemes for Conservation Laws, MUSCL) scheme, and the interface flux is calculated using the AUSM+-up scheme; the viscous flux is discretized using the central difference scheme; to overcome the rigidity problem in thermochemical non-equilibrium flow, the thermochemical non-equilibrium source term is processed using the point-implicit scheme; the time-progression scheme is the LU-SGS (Lower-Upper Symmetric Gauss-Seidel, LU-SGS) scheme, and the specific calculation scheme can be found in (Liu Wei, Zhang Theory, Wang Yongxian, et al. Fundamentals of Parallel Programming in Computational Aerodynamics [M]. Beijing: National Defense Industry Press, 2013). In particular, the equivalent viscosity coefficient (μ) given by equations (8)-(10) in Step 1... e ) and equivalent thermal conductivity (κ) e,tran κ e,rot and κ e,vib ) and the equivalent diffusion coefficient (D) of component s e,sThe values ​​of all these parameters depend on the local shear nonequilibrium parameter Zh in the flow field. To obtain Zh, equations (8) and (13) can be combined to form the following implicit equation:

[0102]

[0103] Where λ is the mean free path of the molecule. Zh NS The viscosity correction function A is obtained directly from the flow field variables and their gradients. Based on equation (15), Zh is first obtained by Newton's iteration method, and then the obtained Zh is substituted into equations (11)-(12) to obtain the viscosity correction function A. μ Parallel dynamic heat conduction correction function A κ The equivalent viscosity coefficient, equivalent translational, rotational, and vibrational thermal conductivity coefficients and the equivalent diffusion coefficient of component s are further calculated using equations (8)-(10).

[0104] Step 3: Based on the discrete scheme in Step 2, perform CFD solution for the calculated shape. The process is as follows:

[0105] 1) Mesh generation: Divide the computational flow field into a structured mesh;

[0106] 2) Set initial parameters, such as incoming flow velocity, temperature, angle of attack, component ratio, etc., and generate initial flow field information;

[0107] 3) Set the numerical format and boundary type for the equation discretization, spatially discretize the governing equations and apply the corresponding boundary conditions, read the chemical reaction database, and solve the flow by time progression.

[0108] 4) During the time progression, output the flow field residual at regular time steps to determine whether the flow field has converged or reached the maximum time iteration step. If it is 'yes', stop the calculation and output the calculation result.

[0109] Step four involves post-processing the results from step three to obtain the flow field around the calculated shape and the aerodynamic / thermal forces on the wall, including flow field quantities (gas density, velocity, temperature, component mass fraction, etc.) and wall pressure p. w The calculation formulas for wall shear stress and wall heat flow are obtained directly from the simulation results as follows:

[0110] τ w,i =τ ij n w,j -(τ jk n w,k n w,j )n w,i (16)

[0111]

[0112] In the formula, nw,i Let i be the component of the unit normal vector of the wall in the i-direction. The wall heat flow in equation (17) consists of three parts: the first term is the sliding friction heat caused by the work done by the surface force; the second term is the conduction heat of different energy modes of translation, rotation and vibration; and the third term is the component diffusion heat.

[0113] The following are three examples of calculating near-continuous hypersonic thermochemical nonequilibrium flow using the method of the present invention:

[0114] In the first embodiment, the object of calculation is a rarefied hypersonic nitrogen flat plate flow.

[0115] Specifically, in this embodiment, the operating conditions are zero angle of attack, incoming Mach number 10, and Knudsen number 0.10. Calculations are performed using the traditional Neumann method, the method of this invention (DiNS-slip), and the kinetic theory-based DSMC method, respectively. The calculation results are as follows: Figure 2 As shown, a velocity profile comparison at x = 0.5 m is presented. It can be seen that the boundary layer flow field profile calculated by this invention agrees very well with the DSMC results, representing a significant improvement over the traditional NS method.

[0116] In the second specific embodiment, the calculation object is a five-component rarefied hypersonic flow around an air cylinder.

[0117] Specifically, in this embodiment, the operating conditions are an incoming Mach number of 27.84 and a Knudsen number of 0.16. The wall friction coefficient calculated by this invention is compared with the results of the traditional NS thermochemical non-equilibrium model and the DSMC method. Figure 3 As shown in the figure, the wall friction of the cylinder flow calculated by this invention agrees very well with DSMC, which is a significant improvement over the traditional NS method.

[0118] In the third specific embodiment, the calculation object is the five-component rarefied hypersonic HTV2 flow around air.

[0119] Specifically, in this embodiment, the operating condition is a flow Mach number of 25 at a speed of 85 km / h. The wall heat flux coefficient calculated by this invention is compared with the results of the traditional NS thermochemical nonequilibrium model and the DSMC method. Figure 3 As shown in the figure, the predicted heat fluxes on both the windward and windward sides of the HTV2 flow are in excellent agreement with the DSMC, representing a significant improvement over the traditional NS method. In summary, the CFD method proposed in this invention can accurately simulate hypersonic near-continuous flow, including boundary layer profiles, wall resistance, and heat flux, and can be rapidly applied in engineering.

[0120] The features described and / or illustrated above with respect to one embodiment may be used in the same or similar manner in one or more other embodiments, and / or in combination with or in lieu of features in other embodiments.

[0121] It should be emphasized that the term "including / comprises" as used herein refers to the presence of a feature, whole, step, or component, but does not exclude the presence or addition of one or more other features, wholes, steps, components, or combinations thereof.

[0122] The apparatus and methods described above can be implemented in hardware or in combination with software. This invention relates to computer-readable programs that, when executed by a logic component, enable that logic component to implement the apparatus or constituent parts described above, or to implement the various methods or steps described above. This invention also relates to storage media for storing the above programs, such as hard disks, magnetic disks, optical disks, DVDs, flash memory, etc.

[0123] Many features and advantages of these embodiments are apparent from this detailed description, and therefore the appended claims are intended to cover all such features and advantages of these embodiments that fall within their true spirit and scope. Furthermore, since many modifications and alterations will readily occur to those skilled in the art, the embodiments of the invention are not intended to be limited to the precise structures and operations illustrated and described, but rather to encompass all suitable modifications and equivalents falling within their scope.

[0124] The parts of this invention not described in detail are techniques known to those skilled in the art.

Claims

1. A CFD method for numerical simulation of near-continuous hypersonic high-temperature thermochemical nonequilibrium flows, comprising the following steps: Step 1: Introduce the Park dual-temperature model into the NS equation to describe the thermochemical non-equilibrium effect at high temperature. Use an equivalent constitutive model to close the shear stress tensor, heat flux vector and mass diffusion flux. Apply corresponding slip boundary conditions on the wall to capture velocity slip and temperature jump phenomena. Obtain macroscopic governing equations and boundary conditions applicable to near-continuous hypersonic high-temperature thermochemical non-equilibrium flow. Step 2: Discretize and solve the macroscopic governing equations and boundary conditions given in Step 1 using the finite volume method; Step 3: Post-process the solution results from Step 2 to obtain the flow field around the calculated shape, the aerodynamic forces on the wall, and the heat flux on the wall.

2. The CFD method for numerical simulation of near-continuous hypersonic high-temperature thermochemical nonequilibrium flow according to claim 1, characterized in that, The method for step one is as follows: 1) The Navier-Stokes equations of the three-dimensional compressible multi-component Park two-temperature model describing thermochemical nonequilibrium flow in Cartesian coordinates can be written in the following conservation form. Where U is the conserved variable; t is time; (x,y,z) are Cartesian coordinates; F x F y and F z The inviscid fluxes in the x, y, and z directions, respectively; F v,x F v,y and F v,z denoted as viscous fluxes in the x, y, and z directions, respectively; S is the thermochemical nonequilibrium source term, where, U=[ρ s ,ρu,ρv,ρw,ρE,ρe ve ] T , F x =[ρ s u,ru 2 +p,ρuv,ρuw,(ρE+p)u,ρe ve [u] T , F y =[ρ s v,pvu,pv 2 +p,ρvw,(ρE+p)v,ρe ve [v] T , F z =[ρ s w,rwu,rwv,rw 2 +p,(ρE+p)w,ρe ve [w] T , In the formula, ρ s It is the density of component s, where s ranges from 1 to N. s Values, N s ρ is the total number of components in the system; p is the total density of the gas mixture; (u, v, w) are the velocity components in the x, y, and z directions, respectively; and E is the total energy of the gas mixture. e is the internal energy per unit mass of the gas mixture, which is composed of the translational-rotational energy of the gas mixture. tr Vibration-electron energy e ve and the zero-point formation enthalpy Δh of the mixed gas f,0 Composition; h s e is the absolute enthalpy per unit mass of component s; ve,s Let the vibrational-electronic energy of component s be denoted as , and the electronic energy e be neglected. el,s The influence of, and on the vibrational energy e of monatomic gases. vib,s It is zero; J s τ is the mass diffusion flux vector of component s, τ is the shear stress tensor, and q is the mass diffusion flux vector of component s. tr and q ve These are the translational-rotational and vibrational-electron heat flow vectors, respectively. It is the mass formation rate of component s. The source term in the vibrational energy equation for a mixed gas; the mass generation rate of component s. It can be written as: Among them, M s υ is the molar mass of component s; sr and υ′ sr These are the stoichiometric coefficients of the reactants and products in the r-th reaction, respectively; [X] s ] r k is the molar concentration of the r-th reactant component s; fr and k br These are the forward and reverse reaction rate coefficients for the r-th reaction, respectively; the subscript "r" indicates the r-th chemical reaction, "f" indicates the forward reaction, and "b" indicates the reverse reaction; the forward reaction rate coefficient k of the r-th reaction is... f The calculation is performed using the modified Arrhenius equation: In the formula, C fr η is the exponential factor of the positive reaction of the r-th reaction. fr θ is the forward reaction temperature exponent of the r-th reaction. dr T is the characteristic temperature of the dissociation of reactant particles in the r-th reaction. tr The translational-rotational temperature; Z r It is the vibration-dissociation coupling nonequilibrium factor, used to characterize the effect of vibration nonequilibrium on the dissociation rate, expressed as: In the formula, T a The geometric mean temperature is defined as follows: Among them, T ve For vibrational-electron temperature, k br The reverse reaction rate coefficient; 2) Introduce a linear constitutive correction model to close the above conservation equations (1) and (2), and give the shear stress tensor τ and the translational-rotational heat flux vector q. tr Vibrational-electron heat flow vector q ve And the mass diffusion flux vector J of component s s Represented in terms of components: Among them, u i These are the components of the velocity vector (u,v,w); x i These are the components of the spatial coordinates (x, y, z); the subscripts i, j, k can all take the values ​​1, 2, 3, (x1, x2, x3) = (x, y, z) are Cartesian coordinates, (u1, u2, u3) = (u, v, w) are velocities, and δ is the velocity. ij For Kronecker notation; y s μ is the mass fraction of component s; e , κ e,tran , κ e,rot and κ e,vib These are the equivalent viscosity coefficient, equivalent translational, rotational, and vibrational thermal conductivity coefficients of the gas mixture, and D. e,s The effective diffusion coefficient of component s is given by the DSMC method for high-speed thermochemical non-equilibrium shear flow analysis, and their calculation formulas are as follows: m e =A μ ·m (8) k e,tran =A κ ·k tran ,k e,rot =A μ ·k rot ,k e,vib =A μ ·k vib (9) D e,s =(k e,tran +k e,rot )Le / (ρc p,tr ) (10) Where μ is the viscosity coefficient of the gas mixture; κ tran , κ rot and κ vib , respectively, represent the translational, rotational, and vibrational thermal conductivity coefficients of the gas mixture; Le is the Lewis number; c p,tr A represents the translational-rotational isobaric specific heat capacity of the gas mixture; μ and A κ These are the correction functions for viscous and translational thermal conductivity coefficients, expressed as functions of the local shear nonequilibrium parameter Zh; 3) The generalized Maxwell-Smoluchowski model is used to describe wall velocity slip and temperature jump phenomena on solid walls. It is compatible with the above-mentioned linear constitutive correction model and Park's two-temperature thermochemical nonequilibrium model. For three-dimensional flow, its specific form is as follows: Among them, u slip It is the tangential sliding velocity vector of the wall, T w For wall temperature; σ u , σ T and σ T,ve These are the tangential momentum adaptation coefficient, translational-rotational, and vibrational-electron energy regulation coefficients, respectively; γ is the specific heat ratio of the gas mixture; c v,tr and c v,vib These are the translational-rotational and vibrational isochoric specific heat capacities of the gas mixture, respectively; τ w This represents the tangential shear stress vector of the wall. 4) The above conservation equation (1) and the linear constitutive correction model (7) combined with the wall slip boundary condition (14) constitute the macroscopic control equation and boundary conditions describing near-continuous hypersonic high-temperature thermochemical non-equilibrium flow.

3. The CFD method for numerical simulation of near-continuous hypersonic high-temperature thermochemical nonequilibrium flow according to claim 2, characterized in that, In sub-step 1) of step one, The source term in the vibrational energy equation of the mixed gas in formula (2) Represented as: In the formula, the first term represents the energy exchange between the translational-rotational modes and vibrational modes of the molecular component calculated according to the Landau–Teller equation; the second term represents the vibrational energy change caused by the gas chemical reaction; τ vib,s It is the molar average relaxation time between the translational-rotational mode and the vibrational mode of component s.

4. The CFD method for numerical simulation of near-continuous hypersonic high-temperature thermochemical nonequilibrium flow according to claim 2, characterized in that, In sub-step 1) of step one, the constant a in formula (5) is taken as 0.5 or 0.7; the reverse reaction rate coefficient k br Or through the equilibrium constant K eq,r According to k br =K eq,r / k fr Formula calculation or based on k br The fitting equation was obtained.

5. The CFD method for numerical simulation of near-continuous hypersonic high-temperature thermochemical nonequilibrium flow according to claim 3, characterized in that, In sub-step 2) of step one, A μ and A κ It is expressed as a function of the local shear nonequilibrium parameter Zh as follows: In the formula, the expression for Zh is: in It is the most probable velocity of local molecular thermal motion, and R is the gas mixture constant. For gas rate The derivative along the streamline normal n; λ e The equivalent molecular mean free path, 6. The CFD method for numerical simulation of near-continuous hypersonic high-temperature thermochemical nonequilibrium flow according to claim 5, characterized in that, In step two, For the equivalent viscosity coefficient (μ) given by equations (8)-(10) in step one e ) and equivalent thermal conductivity (κ) e,tran , κ e,rot and κ e,vib ) and the equivalent diffusion coefficient (D) of component s e,s The calculations all rely on the local shear nonequilibrium parameter Zh in the flow field. The calculations of equations (8) and (13) together form the following implicit equation. Where λ is the mean free path of the molecule. Zh NS It is calculated from the flow field variables and their gradients.

7. The CFD method for numerical simulation of near-continuous hypersonic high-temperature thermochemical nonequilibrium flow according to claim 6, characterized in that, In step two, the calculation method of equation (15) is as follows: Zh is obtained by Newton's iteration method, and the obtained Zh is substituted into equations (11)-(12) to obtain the viscosity correction function A. μ Parallel dynamic heat conduction correction function A κ The equivalent viscosity coefficient, equivalent translational, rotational, and vibrational thermal conductivity coefficients and the equivalent diffusion coefficient of component s are calculated using equations (8)-(10).

8. The CFD method for numerical simulation of near-continuous hypersonic high-temperature thermochemical nonequilibrium flow according to claim 2, characterized in that, In step three, the flow field around the exterior includes gas density, velocity, temperature, and component mass fraction, obtained from simulation results; the wall aerodynamic forces include wall pressure p. w and wall tangential shear stress, wall pressure p w The simulation results yield the following formulas for calculating the tangential shear stress and heat flow at the wall: t w,i =t ij n w,j -(t jk n w,k n w,j )n w,i (16) In the formula, n w,i The component of the unit normal vector of the wall in the i direction; the right-hand side of the calculation of wall heat flow in equation (17) contains three parts: the first term is the sliding friction heat caused by the work done by the surface force; the second term is the conduction heat of different energy modes of translation, rotation and vibration; and the third term is the component diffusion heat.

Images