A CFD Method for Hypersonic Near-Continuous Flow Simulation

Abstract

The present invention relates to a CFD method for hypersonic near - continuum flow simulation, comprising the following steps: Step 1, based on the macroscopic conservation equations, introducing a linear constitutive correction model to close the shear stress and heat flux, and introducing wall slip conditions to describe the wall velocity slip and temperature jump phenomena, so as to obtain the macroscopic control equations and boundary conditions applicable to hypersonic near - continuum flow; Step 2, discretizing the macroscopic control equations obtained in Step 1 using a finite - difference scheme; Step 3, based on the discretization format and boundary treatment in Step 2, performing CFD solution for the computational configuration; Step 4, post - processing the solution results in Step 3 to obtain the flow field around the computational configuration, the wall aerodynamic force and heat, and the calculation of the wall heat flux takes into account the sliding friction work term.

Description

Technical Field

[0001] The present invention relates to the field of rarefied gas dynamics, and in particular to a computational fluid dynamics (CFD) method for hypersonic near-continuous flow simulation. Background Art

[0002] The flight Mach number of a near-space hypersonic vehicle is usually between 5 and 20, and the flight altitude is between 20 and 100 km. During flight, due to the rarefied oncoming gas and strong shear in the boundary layer, local rarefied gas effects often exist in its flow field around the vehicle, which will seriously affect the aerodynamic characteristics and thermal protection design of the vehicle. Accurately and quickly predicting the wall aerodynamic force and heat of the vehicle is crucial for vehicle design.

[0003] For continuous medium flow, the Navier-Stokes-Fourier (NSF) equations can usually be used for accurate description and solved by CFD methods. However, as the gas becomes rarefied, the continuous medium hypothesis gradually fails, making the linear constitutive relation and no-slip boundary hypothesis of the NSF equations no longer valid. At this time, the errors of the vehicle drag and heat flux predicted by traditional CFD methods gradually increase. Generally speaking, the kinetic theory Boltzmann equation can be directly solved or the DSMC particle method can be used to directly simulate rarefied gas flow. However, in the near-continuous flow regime, the computational amount and storage required by kinetic theory methods are often very large, exceeding the limit that can be borne by engineering applications. Other rarefied gas simulation methods also include solving higher-order hydrodynamic equations, such as Grad moment equations and Burnett equations. However, such higher-order equations may have instability problems or it is difficult to propose higher-order boundary conditions by themselves, which also limits their application in engineering. Generally speaking, for the flow field simulation of hypersonic near-continuous flow and the prediction of wall aerodynamic force and heat, traditional CFD methods are no longer accurate and kinetic theory methods require huge computational amounts. Currently, there is a lack of a method with accurate models and computational amounts that meet engineering practice. Summary of the Invention

[0004] The purpose of the present invention is to provide a CFD method that can be used for hypersonic near-continuous flow simulation and whose computational amount meets engineering practice.

[0005] The technical solution of the present invention is as follows:

[0006] A CFD method for hypersonic near-continuous flow simulation includes the following steps:

[0007] Step 1, based on the macroscopic conservation equations, introduce a linear constitutive correction model to close the shear stress and heat flux, and introduce a wall slip condition to describe the wall velocity slip and temperature jump phenomena, so as to obtain the macroscopic control equations and boundary conditions applicable to hypersonic near-continuous flow;

[0008] Step 2: Discretize the macroscopic control equation obtained in Step 1 using the finite difference scheme;

[0009] Step 3: Perform CFD solution for the computational configuration based on the discretization format and boundary treatment in Step 2;

[0010] Step 4: Post-process the solution results in Step 3 to obtain the flow field around the computational configuration, the wall aerodynamic force and heat, and the wall heat flux calculation takes into account the sliding friction work term.

[0011] Step 1 is specifically as follows:

[0012] 1) The conservation equations describing the motion of compressible gas are

[0013]

[0014] where t is time, ρ is density, T is temperature, p is pressure, E is total energy, the subscripts i, j take 1, 2, 3, (x1, x2, x3) = (x, y, z) are Cartesian coordinates, (u1, u2, u3) = (u, v, w) is velocity, τ ij and q j are the shear stress tensor and heat flux vector respectively, and δ is the Kronecker operator;

[0015] 2) Introduce a linear constitutive correction model to close the shear stress τ and heat flux q in the above conservation equations, which is

[0016]

[0017] where μ e and κ e are the equivalent viscosity coefficient and equivalent thermal conductivity coefficient respectively, and their calculation formulas are

[0018] μ e = A μ (Zh)·μ, κ e = A κ (Zh)·κ (3)

[0019] where μ and κ are the gas viscosity coefficient and thermal conductivity coefficient respectively, Zh is the characteristic parameter describing the degree of local shear non-equilibrium, A μ and A κ are the viscosity correction function and thermal conductivity correction function respectively, and their calculation formulas are as follows

[0020]

[0021] where coth is the hyperbolic cotangent function, and the coefficients a1, a2, a3 in the above formula are obtained by simulating rarefied shear non-equilibrium flow using the DSMC method. The expression of Zh is

[0022]

[0023] Among them is the most probable velocity of local molecular thermal motion, R is the gas constant, is the gas velocity is the derivative along the normal direction n of the streamline, λ is the mean free path of molecules, and the calculation formula is

[0024]

[0025] 3) The velocity slip and temperature jump phenomena at the wall are described by applying the Maxwell slip boundary condition on the solid wall surface, which is

[0026]

[0027] Let u be the tangential coordinate along the wall, and y be the normal coordinate along the wall. Among them, u slip and T jump are the gas slip velocity and jump temperature respectively, u w and T w are the wall velocity and temperature respectively, σ u and σ T are the wall tangential momentum accommodation coefficient and thermal accommodation coefficient respectively, and λ T is the mean free path of molecules in an equivalent form:

[0028]

[0029] Among them, γ is the specific heat ratio, and c v is the specific heat at constant volume;

[0030] 4) The above conservation equation (1) combined with the linear constitutive correction models (2)-(4) and the wall slip condition (7) constitutes the macroscopic control equation and boundary conditions for describing hypersonic near-continuous flow.

[0031] In step two, the wall slip boundary condition obtained in step one is processed by unilateral implicit differencing. The method is as follows:

[0032]

[0033] Among them, u1 and T1 are the velocity and temperature of the first grid layer close to the wall, Δy is the grid spacing of the first layer, and C u =λ(2 - σ u ) / σ u , and C T =λ T (2 - σ T ) / σ T .

[0034] In step two, the time marching adopts the LU - SGS format.

[0035] The process of Step 3 is as follows:

[0036] 1) Mesh generation, dividing the computational flow field into structured body-fitted meshes;

[0037] 2) Set initial parameters, including the incoming flow Mach number, Reynolds number, angle of attack, boundary type, and set the initial flow field;

[0038] 3) Discretize the macroscopic control equations in space and apply boundary conditions, and iterate along time;

[0039] 4) During the time marching process, judge whether the flow field converges or reaches the maximum calculation duration every certain number of steps. If 'yes', stop the calculation and output the calculation results.

[0040] In Step 4, the calculation formulas for wall shear stress and wall heat flux are as follows

[0041]

[0042] Advantages of the present invention compared with the prior art:

[0043] (1) The present invention establishes a CFD method applicable to hypersonic near-continuous flow simulation, which is more suitable for complex engineering practical problems in near-space flight compared with existing methods;

[0044] (2) The present invention takes into account the influence of local rarefaction effect of the boundary layer, and significantly improves the prediction ability for near-continuous flow compared with traditional CFD methods;

[0045] (3) The present invention has high computational efficiency, which is more than one order of magnitude faster than the kinetic theory DSMC particle simulation method;

[0046] (4) The present invention is easy to implement and can be conveniently implanted into any existing CFD solver. Description of the Drawings

[0047] Figure 1 is the block diagram of the steps of the present invention

[0048] Figure 2 is the velocity profile at x = 25 mm of the hypersonic rarefied flat plate flow calculated by the present invention and the comparison with other methods

[0049] Figure 3 is the wall cumulative drag of the hypersonic rarefied flat plate flow calculated by the present invention and the comparison with other methods

[0050] Figure 4 is the wall heat flux coefficient of the cylinder flow at Mach 10 and Knudsen number 0.05 calculated by the present invention and the comparison Detailed Implementation Modes

[0051] The present invention will be described in detail below in conjunction with specific embodiments and the accompanying drawings. In the following description, for purposes of explanation rather than limitation, specific details are set forth in order to provide a thorough understanding of the present invention. However, it will be apparent to those skilled in the art that the present invention may be practiced in other embodiments without these specific details.

[0052] It should be noted here that in order to avoid obscuring the present invention with unnecessary details, only the device structures and / or processing steps closely related to the solution according to the present invention are shown in the drawings, while other details less relevant to the present invention are omitted.

[0053] The embodiment of the present invention provides a CFD method for hypersonic near - continuous flow simulation. This method overcomes the deficiencies of the prior art and mainly solves the technical problem of: overcoming the failure problem of existing CFD methods for high - speed rarefied gas flow, giving a linear constitutive correction model that can consider the local rarefied gas effect of the boundary layer and the corresponding wall - slip boundary conditions, and achieving stable and high - efficiency solution. The implementation of the present invention is as Figure 1 shown, including the following steps:

[0054] Step 1: Based on the macroscopic conservation equations, introduce a linear constitutive correction model to close the shear stress and heat flux, and introduce wall - slip conditions to describe the wall velocity slip and temperature jump phenomena. Thus, the macroscopic control equations and boundary conditions applicable to hypersonic near - continuous flow are obtained, including the following:

[0055] 1) The conservation equations describing the motion of compressible gas are

[0056]

[0057] where t is time, ρ is density, T is temperature, p is pressure, E is total energy, the subscripts i, j take 1, 2, 3, (x1, x2, x3) = (x, y, z) are Cartesian coordinates, (u1, u2, u3) = (u, v, w) is velocity, τ ij and q j are the shear - stress tensor and heat - flux vector respectively, and δ is the Kronecker operator;

[0058] 2) Introduce a linear constitutive correction model to close the shear stress τ and heat flux q in the above conservation equations, which is

[0059]

[0060] where μ e and κ e are the equivalent viscosity coefficient and equivalent heat - conduction coefficient respectively, and their calculation formulas are

[0061] μ e = A μ(Zh)·μ, κ e =A κ (Zh)·κ (3)

[0062] where μ and κ are the gas viscosity coefficient and thermal conductivity coefficient, respectively, Zh is the characteristic parameter describing the degree of local shear non-equilibrium, and A μ and A κ They are the viscosity correction function and the heat conduction correction function, and the calculation formula is as follows

[0063]

[0064] Where coth is the hyperbolic cotangent function. The coefficients a1, a2, and a3 in the above formula are obtained by simulating a large number of rarefied shear non-equilibrium flows using the DSMC method. The expression of Zh is:

[0065]

[0066] in is the most probable speed of local molecular thermal motion, R is the gas constant, is the gas velocity The derivative along the streamline normal n, λ is the molecular mean free path, and the calculation formula is

[0067]

[0068] In this step, the viscosity coefficient μ and the thermal conductivity coefficient κ in formula (3) are only related to temperature and can be calculated by the power law formula:

[0069]

[0070] Where T ref is the reference temperature, μ ref is the viscosity coefficient at the reference temperature, ω is the temperature index, Pr is the Prandtl number, c p is the specific heat at constant pressure.

[0071] In this step, the equivalent viscosity coefficient μ in formula (3) is e and the equivalent thermal conductivity κ e The calculation requires simultaneous equations (3)(5)(6) to form the following implicit equation:

[0072]

[0073] Among them Zh NS It is directly calculated from the flow field variables and their gradients. According to equation (8), the Newton iteration method is first used to calculate Zh, and then the obtained Zh is substituted into equation (4) to obtain the viscosity correction function A and the heat conduction correction function A, and the equivalent viscosity coefficient and the equivalent heat conduction coefficient are further calculated by equation (3).

[0074] Furthermore, if the Newton iteration method is used to solve the equivalent viscosity coefficient and the equivalent heat conduction coefficient in each time advancement process, the computational efficiency will be reduced. The viscosity correction function and the heat conduction correction function can be pre-solved according to Equation (8) and made into an interpolation table, and then directly interpolated and used during the time advancement process, thus avoiding iteration in each step.

[0075] In this step, the velocity gradient in Equation (5) is calculated in the following manner. For two-dimensional flow, the direction of the streamline is s = (u / V, v / V), and the gradient of the velocity along the normal direction of the streamline can be directly written as

[0076]

[0077] For three-dimensional flow, the direction of the streamline is s = (u / V, v / V, w / V), and there are usually infinitely many directions perpendicular to the streamline. At this time, the maximum gradient of the velocity along the normal direction of the streamline can be given according to the velocity gradient and its derivative along the streamline

[0078]

[0079] Furthermore, for the simulation of rarefied boundary layer flow, since the gradient of the velocity along the normal direction of the streamline in the boundary layer is generally much larger than its gradient along the streamline, the local velocity gradient can be used to approximate the gradient of the velocity along the normal direction of the streamline at this time, that is In this way, the calculation of Zh NS is more convenient and direct.

[0080] 3) Apply the Maxwell slip boundary condition on the solid wall surface to describe the velocity slip and temperature jump phenomena at the wall surface, which is

[0081]

[0082] Here, it is assumed that u is along the tangential direction of the wall surface, and y is the coordinate along the normal direction of the wall surface, where u slip and T jump are the gas slip velocity and the jump temperature respectively, u w and T w are the wall surface velocity and temperature respectively, σ u and σ T are the wall surface tangential momentum accommodation coefficient and the heat accommodation coefficient respectively, and λ T is the molecular mean free path in an equivalent form

[0083]

[0084] where γ is the specific heat ratio, and c v is the specific heat at constant volume.

[0085] 4) The above conservation equations (1), combined with the linear constitutive correction models (2)-(4) and the wall slip condition (11), constitute the macroscopic control equations and boundary conditions for describing hypersonic near-continuous flows.

[0086] Step 2: Discretize the control equations obtained in Step 1 using a high-precision finite-difference scheme. Among them, the convection term is discretized using the fifth-order WENO scheme, the viscous term is discretized using the sixth-order central-difference scheme, and the time marching is performed using the LU-SGS scheme. The specific calculation scheme can be found in (Yan Chao. Computational Fluid Dynamics Methods and Applications [M]. Beijing University of Aeronautics and Astronautics Press, 2006). In addition, the wall slip boundary condition obtained in Step 1 is processed using a one-sided implicit difference to ensure the stability of the calculation.

[0087] In this step, other schemes can also be used for spatial discretization and time marching, such as using the finite volume scheme for spatial discretization and the Runge-Kutta scheme for time marching.

[0088] In this step, if the wall slip condition is implemented directly by Equation (11), it may cause calculation instability. Here, a one-sided implicit difference scheme is adopted, that is, let

[0089]

[0090] Substituting into Equation (11), the implementation form of the slip boundary is obtained as

[0091]

[0092] where u1 and T1 are the velocity and temperature of the first grid layer near the wall, Δy is the grid spacing of the first layer, and C u = λ(2 - σ u ) / σ u , C T = λ T (2 - σ T ) / σ T .

[0093] Step 3: Based on the discretization scheme and boundary treatment in Step 2, perform CFD solutions for the given calculation configuration. The process is as follows:

[0094] 1) Mesh generation, dividing the computational flow field into structured body-fitted meshes;

[0095] 2) Set the initial parameters, such as the freestream Mach number, Reynolds number, angle of attack, boundary type, etc., and set the initial flow field;

[0096] 3) Perform spatial discretization of the control equations and apply the boundary conditions, and iterate along time;

[0097] 4) During the time advancement process, every certain number of steps is used to determine whether the flow field converges or reaches the maximum calculation duration. If 'yes', the calculation is stopped and the calculation results are output.

[0098] In this step, the calculated shape can be a simple shape, such as a flat plate, a cylinder, a blunt cone, or it can also be a three-dimensional complex shape or a flight shape.

[0099] Step 4: Post-process the solution results of Step 3 to obtain the flow field around the calculated shape and the aerodynamic force / heat on the wall surface. Among them, the flow field quantities (gas density, velocity, temperature, etc.) and the wall pressure p w are directly obtained from the simulation results. The calculation formulas for the wall shear stress and the wall heat flux are as follows

[0100]

[0101] Here, the calculation of the wall heat flux takes into account the sliding friction work term, that is, u slip τ w .

[0102] The following gives two embodiments of calculating hypersonic near-continuous flow using the method of the present invention:

[0103] In the first embodiment, the calculation object is the rarefied hypersonic flat plate flow.

[0104] Specifically, the working conditions of this embodiment are zero angle of attack, and the incoming flow Mach numbers are 12.7 and 15 respectively. The traditional NSF method, the method of the present invention (DiNS), and the kinetic theory DSMC method are used for calculation respectively. The calculation results are as Figure 2 and Figure 3 shown. The velocity profile comparison at x = 25 mm for the Mach 12.7 case and the wall drag comparison for the Mach 15 case are given respectively. It can be seen that the boundary layer flow field profile and the wall drag calculated by the present invention are in good agreement with the DSMC results, and are greatly improved compared with the traditional NSF method. In addition, in terms of calculation time, both the traditional NSF method and the method of the present invention are about 2.5 hours, and the calculation time of the DSMC method is about 41 hours. The calculation efficiency of the present invention is more than one order of magnitude faster than the DSMC particle simulation method. For three-dimensional near-continuous flow, the improvement of this calculation efficiency will be more obvious, reaching more than two or three orders of magnitude.

[0105] In the second specific embodiment, the calculation object is the rarefied hypersonic flow around a cylinder.

[0106] Specifically, the working conditions of this embodiment are the incoming flow Mach number of 10 and the Knudsen number of 0.05. The wall heat flux coefficient calculated by the present invention and the comparison with the results of the traditional NSF equation and the DSMC method are as Figure 4As shown, it can be seen that the wall heat flux of the cylinder flow calculated by the present invention is in good agreement with DSMC, and is greatly improved compared with the traditional NSF method. Generally speaking, the CFD method proposed by the present invention can accurately simulate hypersonic near-continuous flow, including boundary layer profiles and wall drag and heat flux, and is much higher than the DSMC particle simulation method in terms of computational efficiency, and can quickly achieve engineering applications.

[0107] Features described and / or illustrated above for one embodiment can be used in the same or similar manner in one or more other embodiments, and / or combined with or replace features in other embodiments.

[0108] It should be emphasized that the term "comprising / including" as used herein refers to the presence of features, whole things, steps or components, but does not exclude the presence or addition of one or more other features, whole things, steps, components or combinations thereof.

[0109] The above devices and methods of the present invention can be implemented by hardware, or can be implemented by hardware in combination with software. The present invention relates to such a computer-readable program, which when executed by a logic component, can enable the logic component to implement the above-mentioned device or component, or enable the logic component to implement the above-mentioned various methods or steps. The present invention also relates to a storage medium for storing the above program, such as a hard disk, a magnetic disk, an optical disk, a DVD, a flash memory, etc.

[0110] Many features and advantages of these embodiments are clear from this detailed description, so the appended claims are intended to cover all such features and advantages of these embodiments that fall within their true spirit and scope. In addition, since many modifications and changes are easily conceivable by those skilled in the art, the embodiments of the present invention are not to be limited to the exact structures and operations illustrated and described, but may cover all suitable modifications and equivalents falling within their scope.

[0111] The parts not detailed in the present invention are well-known technologies to those skilled in the art.

Claims

1. A CFD method for hypersonic near - continuous flow simulation, comprising the following steps: Step 1: Based on the macroscopic conservation equations, introduce a linear constitutive correction model to close the shear stress and heat flux, and introduce the wall slip condition to describe the wall velocity slip and temperature jump phenomena, obtaining the macroscopic control equations and boundary conditions applicable to hypersonic near - continuous flow, as follows: 1) The conservation equations describing the motion of compressible gas are where \(t\) is time, \(\rho\) is density, \(T\) is temperature, \(p\) is pressure, \(E\) is total energy, the subscripts \(i,j\) take values \(1,2,3\), \((x_1,x_2,x_3)=(x,y,z)\) are Cartesian coordinates, \((u_1,u_2,u_3)=(u,v,w)\) is velocity, \(\tau\) ij and \(q\) j are the shear stress tensor and the heat flux vector respectively, \(\delta\) is the Kronecker operator; 2) Introduce a linear constitutive correction model to close the shear stress τ and heat flux q in the above conservation equations, which is where μ e and κ e are the equivalent viscosity coefficient and the equivalent heat conduction coefficient respectively, and their calculation formulas are μ e = A μ (Zh)·μ, κ e = A κ (Zh)·κ (3) where μ and κ are the gas viscosity coefficient and the thermal conductivity coefficient respectively, Zh is the characteristic parameter describing the degree of local shear non-equilibrium, A μ and A κ are the viscosity correction function and the thermal conductivity correction function respectively, and their calculation formulas are as follows where coth is the hyperbolic cotangent function, and the coefficients a1, a2, a3 in the above formula are obtained by simulating rarefied shear non - equilibrium flow using the DSMC method. The expression of Zh is wherein is the most probable velocity of local molecular thermal motion, R is the gas constant, is the gas velocity is the derivative along the normal n of the streamline, λ is the mean free path of molecules, and the calculation formula is For three - dimensional flow, the streamline direction is s=(u / V, v / V, w / V), and the maximum gradient along the normal direction of the streamline is given according to the velocity gradient and its derivative along the streamline: The equivalent viscosity coefficient μ in Equation (3) e and the equivalent heat conduction coefficient κ e are calculated by simultaneously solving Equations (3), (5), and (6) to form the following implicit equation: where Zh NS is directly calculated from the flow field variables and their gradients. According to Equation (11), first, the Newton iteration method is used to find Zh, and then the obtained Zh is substituted into Equation (4) to find the viscous correction function A μ and the heat conduction correction function A κ , and further, the equivalent viscous coefficient and the equivalent heat conduction coefficient are calculated by Equation (3); First, solve the viscous correction function and heat conduction correction function according to equation (8) and make an interpolation table, and then directly interpolate and use it during the time - marching process, thus avoiding iteration at each step; 3) Describe the velocity slip and temperature jump phenomena at the wall by applying the Maxwell slip boundary condition on the solid wall, which is Let \(u\) be the tangential coordinate along the wall surface and \(y\) be the normal coordinate along the wall surface, where \(u\) slip and \(T\) jump are the gas slip velocity and the jump temperature respectively, \(u\) w and \(T\) w are the wall velocity and the wall temperature respectively, \(\sigma\) u and \(\sigma\) T are the tangential momentum accommodation coefficient and the thermal accommodation coefficient of the wall surface respectively, \(\lambda\) T is the molecular mean free path in an equivalent form: where γ is the specific heat ratio and c v is the specific heat capacity at constant volume; 4) The above conservation equation (1) combined with the linear constitutive correction model (2) - (4) and the wall slip condition (7) constitutes the macroscopic control equations and boundary conditions describing hypersonic near - continuous flow; Step 2: Discretize the macroscopic control equations obtained in Step 1 using a finite - difference scheme, and use a one - sided implicit difference to process the wall slip boundary condition, as follows: where u1 and T1 are the velocity and temperature of the first grid layer near the wall, Δy is the first grid layer spacing, and C u = λ(2 - σ u ) / σ u , C T = λ T (2 - σ T ) / σ T ; Step 3: Based on the discretization scheme and boundary treatment in Step 2, perform CFD solution for the computational configuration; Step 4: Post - process the solution results in Step 3 to obtain the flow field around the computational configuration, the wall aerodynamic force and heat, and the calculation of the wall heat flux considers the sliding friction work term.

2. The CFD method for hypersonic near-continuous flow simulation according to claim 1, characterized in that, In Step 2, the LU - SGS format is used for time - marching.

3. The CFD method for hypersonic near-continuous flow simulation according to claim 1, wherein The process of Step 3 is as follows: 1) Mesh generation, divide the computational flow field into structured body - fitted meshes; 2) Set the initial parameters, including the incoming - flow Mach number, Reynolds number, angle of attack, boundary type, and set the initial flow field; 3) Discretize the macroscopic control equations in space and apply boundary conditions, and perform iteration along time; 4) During the time - marching process, judge whether the flow field converges or reaches the maximum calculation duration every certain number of steps. If 'yes', stop the calculation and output the calculation results.

4. The CFD method for hypersonic near - continuous flow simulation according to claim 1, characterized in that, In Step 4, the calculation formulas for the wall shear stress and wall heat flux are as follows

Images