CN115952748A - Numerical calculation method for two-dimensional rarefied gas flow - Google Patents

Abstract

The invention belongs to the technical field of rarefied gas dynamics, and particularly relates to a numerical calculation method for two-dimensional rarefied gas flow. The specific technical scheme is as follows: and introducing a DG numerical method into an NCCR equation to perform numerical discrete calculation, and introducing a limiter to perform numerical discontinuity detection and limitation. The method solves the technical problems that the traditional NS equation can not solve the technical problems of discontinuous flow and large calculation amount of a molecular dynamic model in transition flow and slip flow, and has great advantages for solving the problems of ultrahigh-speed flow and thermodynamic calculation of various aircrafts in the near space, particularly shock waves.

Description

A numerical method for two-dimensional rarefied gas flows

Technical Field

The invention belongs to the technical field of rarefied gas dynamics, and in particular relates to a numerical calculation method for two-dimensional rarefied gas flow.

Background Art

The study of traditional fluid mechanics equations focuses on the continuous region, that is, the mean free path of fluid molecules is much smaller than the characteristic size of the flow. The Navier-Stokes-Fourier equations (NSF equations) are a typical representative of this type of equations. It closes the three major flow conservation equations of continuity, momentum, and energy by coupling the constitutive relationship of Newtonian viscous stress and Fourier heat conduction with the Stokes hypothesis. In the past two centuries, the NSF equation has achieved great success in the field of continuous flow and promoted the development of fluid mechanics. However, when the mean free path of molecules gradually increases, that is, the gas becomes thinner and no longer dense, the continuity assumption fails. At this time, if we continue to use the NSF equation based on continuous media or equilibrium state (when the macroscopic thermal observation quantity of the system no longer changes with time, the system is in equilibrium state), it must be inappropriate.

In recent years, with the continuous progress of science and human exploration of unknown areas, people are paying more and more attention to the field of rarefied gases and their engineering applications, such as the prediction of aerodynamic forces of aircraft flying in the thin atmosphere (tens of kilometers above the ground). Under this condition, the gas density is very low, about ^10-7 to ^10-10 kg/ ^m3 . At this time, the larger distance between molecules will lead to a larger distance traveled by a molecule after two collisions, that is, the average free path of the molecule is larger. As the collision effect of molecules gradually emerges, the traditional macroscopic equation - NSF equation will no longer be applicable. At this time, the Boltzmann equation with collision terms takes on the important task of solving rarefied flow.

According to the size of Knudsen _Kn number (an important parameter to measure the rarefaction of gas and the correctness of macroscopic model), the existing gas flow and its corresponding solution methods can be divided into four types:

In the field of continuous flow (Kn<0.001), the traditional Euler equation or NS equation is applicable;

In the slip flow domain (0.001<Kn<0.1), the NS equations hold in the mainstream outside the boundary, but the velocity slip and temperature jump boundary conditions must be considered at the boundary;

In the transitional flow domain (0.1<Kn<10), the NS equations no longer hold;

For free molecular flow (10>Kn), directly use the molecular kinetic theory method, such as DSMC.

Among the above-mentioned methods for solving fluid problems, the NS equations and the molecular kinetic theory have achieved success in the fields of continuous flow and free molecular flow, respectively. However, for the slip flow domain and the transition flow domain, there is no unified and effective control equation to solve the flow domain in this area. For example, the gas flow in the critical space layer (20-100km above sea level) that has been vigorously developed in recent years happens to be in these two domains. Therefore, it is particularly important to develop a set of gas dynamics theories that are applicable to both continuous flow and rarefied flow.

In order to study the flow of rarefied gases, a large number of mathematical physics equations have been proposed, which can basically be divided into three categories: molecular dynamics models that are solved from a microscopic perspective, including the direct simulation Monte Carlo method (DSMC), the information preservation method of DSMC (DSMC-IP) and the lattice Boltzmann method (LBM); there are also methods that directly perform discrete modeling and solve the distribution function in the Boltzmann equation, such as the BGK model equation; and substituting the macroscopic statistical expression into the microscopic Boltzmann equation to obtain the conservation equations and evolution equations of macroscopic quantities. Representative methods are the Grad method, the Chapman-Enskog method and the Eu method.

For the molecular dynamics model, since it uses simulated molecules to simulate the movement of real molecules and the collision of gas flows, it will consume a lot of computing resources at low Kn numbers (high molecular concentration). For the BGK model equation, it simplifies the collision terms in the Boltzmann equation, so it can only obtain accurate results in equilibrium or near-equilibrium states, and its calculation of transport coefficients is also inaccurate. For example, the Prandtl (Pr) number of monatomic gas calculated by the BGK model is 1, instead of the correct value of 2/3. In the non-NS fluid control equation method evolved from the Boltzmann equation, the R-13 equation based on the moment method proposed by Grad and the Burnett equation based on the Chapman-Enskog expansion are considered to not meet the entropy conditions of the laws of thermodynamics, that is, it is difficult to meet the definite solution of the Boltzmann equation, so the development of these two methods has been restricted to a certain extent. Different from the above two methods, the Eu method strictly satisfies the entropy condition. It starts from the Boltzmann solution condition H theorem, grasps the characteristics of entropy increase from non-equilibrium to equilibrium, constructs an exponential distribution function, and obtains the entropy increase dissipation model from non-equilibrium to equilibrium. When approaching the equilibrium state, the model converges to the Rayleigh-Onsager dissipation function. Based on this, the collision term is handled by constructing a unified nonlinear entropy increase model from non-equilibrium to equilibrium. Therefore, the three major conservation equations and constitutive equations such as viscous stress and heat flux are derived through the Boltzmann equation, and the equations of these conservation quantities and non-conservation quantities are coupled together to complete the closure of the gas dynamics equation to form the Eu unified fluid equation. Afterwards, Myong focused on dealing with the high-order terms of the constitutive equation based on Eu and constructed a nonlinear coupled constitutive relation (NCCR) based on this.

Therefore, if a numerical calculation method for two-dimensional rarefied gas flow can be provided, it will have excellent application prospects.

Summary of the invention

To solve the above technical problems, the present invention provides a numerical calculation method suitable for rarefied gas flow, especially non-equilibrium thermal flow analysis of hypersonic flight with shock wave phenomenon, based on the two-dimensional NCCR equation derived from the Eu method.

To achieve the above-mentioned purpose of the invention, the technical solution adopted by the present invention is: a numerical calculation method for two-dimensional rarefied gas flow, which introduces the DG numerical method into the NCCR equation to perform numerical discrete calculation, and introduces a limiter to perform numerical discontinuity detection and limitation.

Preferably, the NCCR equation is a dimensionless form equation obtained by simplifying dimensionless parameters and similarity criteria numbers.

Preferably: the global solutions _S and _U in the local unit Ω are expressed by approximate solutions Sh and Uh,

为基函数，

Where N is the number of basis functions;

is the basis function,

The DG-NCCR discrete equation is:

Preferably, the limiter is a TVB slope limiter.

Preferably, the method comprises the following steps:

S01. Meshing the two-dimensional solution flow field through a mesh program and generating a mesh file;

S02, read the grid file and record the grid node coordinates;

S03. Obtain standard square grid units according to the node coordinates;

S04, given pressure far-field boundary conditions, Langmuir wall slip boundary conditions;

S05, initializing calculation conditions, determining the time step according to CFL conditions; assigning initial values x ₀ , x _n-1 = x ₀ to flow field parameters according to boundary conditions;

S06, substituting x _n-1 into the numerical integration term in the DG-NCCR discrete equation for calculation;

S07, using different flux calculation formats to calculate the flux in the discrete equation;

S08. According to the calculation results of steps S06 and S07, the discrete equation is solved by using a third-order TVD-RK to obtain flow field parameters x _n .

Preferably: in the step S01, the generated two-dimensional grid is a triangular grid unit; in the step S03, the triangular grid unit in the step S01 is first converted into a standard triangular unit, and then converted into a standard square grid unit.

Preferably, it also includes:

S09, traverse all solution results, find the "problem unit" and use the limiter to limit the solution parameters in the unit, and update x _n ;

S10, calculating the iteration error x _n -x _n-1 , and determining the magnitude of the iteration error and the target calculation accuracy ε:

When x _n -x _n-1 <ε, go to step S11;

When _xn - _xn- 1≥ε, let _xn-1 ＝ _xn , and enter step S06 to perform loop calculation until the calculation result meets the target calculation accuracy requirement ε;

S11, the loop calculation ends.

Corresponding: Numerical methods for two-dimensional rarefied gas flows with applications in near-space vehicle flight.

Corresponding: an electronic device, characterized in that: comprising:

one or more processors;

A storage device for storing one or more programs;

When the one or more programs are executed by the one or more processors, the one or more processors implement a numerical calculation method for two-dimensional rarefied gas flow.

Correspondingly: a computer-readable medium storing a computer program, characterized in that: when the computer program is executed by a processor, a numerical calculation method for two-dimensional rarefied gas flow is implemented.

Compared with the prior art, the present invention has the following beneficial effects:

The present invention introduces the DG method into the NCCR equation for numerical discrete calculation, and introduces a limiter into the numerical method for numerical discontinuous detection and limitation, so that a good solution is obtained for the shock wave problem of ultra-high sound speed in a rarefied state. This method solves the technical problems that the traditional NS equation cannot solve discontinuous flow and the molecular dynamics model has a large amount of calculation in transition flow and slip flow, and has great advantages for ultra-high-speed flow and thermal calculation of various aircraft in near space, especially for solving shock wave problems.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG1 is a flowchart of a two-dimensional DG-NCCR discrete equation calculation method according to the present invention;

FIG2 is an unstructured grid diagram of a two-dimensional flow around a cylinder according to an embodiment of the present invention;

FIG3 is a schematic diagram of coordinate transformation between an arbitrary triangular unit Ω and a standard triangular unit T according to the present invention;

FIG4 is a schematic diagram of coordinate conversion between the standard triangle unit T and the standard square unit R of the present invention;

FIG5 is a schematic diagram of a two-dimensional restriction unit K and its adjacent units according to the present invention.

DETAILED DESCRIPTION

The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Unless otherwise specified, the technical means used in the embodiments are conventional means well known to those skilled in the art.

The present invention proposes a numerical calculation method suitable for two-dimensional rarefied gas flow analysis. The core idea is to introduce the DG numerical method into the NCCR equation to perform numerical discrete calculations, and introduce a limiter to perform numerical discontinuity detection and limitation. The fluid transport equation used in this method is derived from the Boltzmann equation, and the Boltzmann equation is generally considered to be the gas dynamics control equation for the entire flow domain. Therefore, the present invention aims to break through the disadvantages of the traditional NS equation that it cannot solve discontinuous flows, and the large amount of calculation of the molecular dynamics model in transition flows and slip flows, and propose a set of thermoaerodynamic analysis methods suitable for the flight process of near-space vehicles.

The numerical calculation method of the present invention adopts the DG (discontinuous Galerkin) numerical method in the process of discretizing the two-dimensional NCCR equation. The method combines the advantages of traditional finite element (FEM) and finite volume (FVM), specifically including: (1) high adaptability to unstructured grids, suitable for complex configurations; (2) no need to consider continuity when encrypting or thickening the grid, and the adaptive technology of grid encryption can be easily realized; (3) each calculation unit is independent, which has the advantages of easy coding and parallelization; (4) since the discontinuous approximation operation is adopted in the unit, the DG method is very suitable for solving numerical discontinuous problems. The DG-NCCR discrete equation has a good solution advantage for the hypersonic shock wave problem in near space. The shock wave is a compression wave in which the airflow parameters undergo a sudden change in the hypersonic flow process. It is widely present in the hypersonic flight of artillery shells, rockets, aircraft, etc. Since the air in near space is thin and the atmospheric temperature and pressure are lower than those in the low-altitude environment, the hypersonic problem in near space is more difficult to solve and analyze.

In order to better explain the numerical calculation method of the present invention, the source of the two-dimensional NCCR equation is briefly given. For the two-dimensional Boltzmann equation, by introducing the exponential distribution function and its corresponding Rayleigh-Onsager dissipation function, the nonlinear coupling constitutive equation and NCCR constitutive relationship that satisfy the entropy increase principle are obtained:

in,

Where Λ _k is the molecular collision dissipation term, Z _k is the kinetic energy term of fluid streamline effect caused by molecular diffusion, U is the conserved variable, _Finv (U) is the inviscid term, _Fvis (U) is the viscous term, Π is the viscous stress, Q is the heat conduction, and I is the unit tensor.

q(k)=sinhk/k, ^k2 is the Rayleigh-Onsager dissipation function, γ`=(5-3γ)/2 (γ is the specific heat ratio), _ηb is the additional stress viscosity coefficient, and _Cp represents the constant pressure specific heat capacity.

p is pressure, ρ is density, T is temperature, Δ is volume additional normal stress, and E is energy per unit mass.

Here, the velocity and divergence are both two-dimensional. The first equation in the above equation is the three major conservation equations of mass, momentum and energy, and the second equation is the constitutive equation, which gives the expressions of viscous stress Π and heat flux Q, which is used to close the conservation equations.

In order to simplify the calculation, the above equation is dimensionless:

Wherein, the subscript r indicates the reference state, _fb indicates the ratio of bulk viscosity to shear viscosity, which can be measured by acoustic absorption experiment. For monatomic molecules, _fb = 0, and the additional viscous stress is zero. t is time, x is length, L is the characteristic length of the flow field, η is the fluid dynamic viscosity coefficient, λ is the thermal conductivity coefficient, and u is the velocity.

Define the following similarity criteria:

Nδ represents the magnitude of viscous stress relative to static pressure, and its magnitude can measure the degree to which the system is far from equilibrium. Ma is the Mach number, Re is the Reynolds number, Ec is a function of the Mach number Ma, and Pr is the Prandtl number. The dimensionless vector form of the equation group can be obtained through the above dimensionless parameters and similarity criteria:

in,

The solution of equation group (2) is the numerical calculation solution of the two-dimensional rarefied gas flow involved in the scheme of the present invention. By solving the flow parameters such as temperature, velocity, pressure and density of the above equations, the temperature field, velocity field, pressure field and density field of the two-dimensional rarefied flow can be solved, thereby realizing the thermal flow analysis of rarefied hypersonic flight.

A numerical calculation method for two-dimensional rarefied gas flow is used to solve the DG-NCCR discrete equation, which specifically includes the following steps:

S01. Mesh the two-dimensional solution flow field through a mesh program and generate a mesh file. Commercial mesh generation software can be used to generate the mesh file, and the mesh generation software can be Gambit or ICEM. The mesh generated in this step is generally a two-dimensional triangular mesh.

S02, read the grid file and record the grid node coordinate parameters. The two-dimensional DG-NCCR program reads the grid file, and the two-dimensional DG-NCCR program is written with the numerical calculation method for two-dimensional rarefied gas flow disclosed in the present invention.

S03, obtaining standard square grid units according to the node coordinates. Specifically, first convert the triangular grid units in step S01 into standard triangular grid units, and then convert them into standard square grid units.

S04. Given the pressure far-field boundary condition, the velocity is calculated by the given Mach number of the flow, and the other parameters are consistent with the surrounding atmospheric parameters. The wall is set to the Langmuir slip boundary condition, which is based on the Langmuir adsorption isotherm and takes into account the interaction between gas molecules and the wall to determine the velocity and temperature on the wall.

S05. Perform iterative calculation: Initialize the calculation conditions, such as the maximum number of iterations M, the calculation accuracy requirement ε, and determine the time step according to the CFL condition; assign initial values x ₀ , x _n-1 = x ₀ to the flow field parameters according to the boundary conditions. The CFL condition is an important concept in the stability and convergence analysis of the finite difference and finite volume methods. n represents the nth calculation, n = 1, 2, 3, ...

S06. Substitute x _n-1 into the numerical integration term in the DG-NCCR discrete equation for calculation.

S07. Use different flux calculation formats to calculate the flux in the DG-NCCR discrete equation.

S08. According to the calculation results of steps S06 and S07, the DG-NCCR discrete equation is solved by using the third-order TVD-RK to obtain the flow field parameter result x _n of the current solution step.

S09, traverse all solution results, find the "problem unit" and use a limiter to limit the solution parameters in the unit, and update x _n . The limiter is a TVB slope limiter.

S10, calculate the iteration error _xn - _xn-1 , and determine the size of the iteration error and the target calculation accuracy ε: when _xn - _xn-1 <ε, proceed to step S11; when _xn - _xn-1 ≥ε, let _xn-1 = _xn , proceed to step S06 for loop calculation until the calculation result meets the target calculation accuracy requirement ε.

S11, the cycle calculation ends. According to _xn of each unit, the original value of the flow field parameters is solved dimensionlessly, and the result data is substituted into the post-processing software to obtain the flow field cloud map and streamline map, such as tecplot, origin and other related fluid post-processing software, for flow analysis.

Due to the high nonlinearity and coupling of NCCR, it is difficult to obtain a good solution using general discrete methods. The numerical calculation method of the present invention introduces the DG method into the NCCR equation for numerical discrete calculation, and at the same time introduces a limiter into the numerical method for numerical discontinuity detection and limitation, so that the present invention has a good solution to the ultra-high sound shock wave problem in a rarefied state.

Example

First determine the expression of the DG-NCCR discrete equation.

The global solutions _S and _U in the local unit Ω are approximated by the approximate solutions Sh and Uh.

为基函数，其由于在积分上有如下特性而在DG离散方法中被用于数值积分：In the formula,

is a basis function, which is used for numerical integration in the DG discrete method due to its following properties in integration:

It should be noted that the form of the basis functions in each unit is the same. u _i (t) and s _i (t) are the coefficients to be solved of the approximate solutions _Sh and U _h in the unit, respectively. Knowing their values, the approximate solution of the unit is obtained.

N is the number of basis functions, and its size is restricted by the calculation order I. The specific relationship is:

In order to ensure the calculation accuracy and save computing resources, the preferred value of the order I is 2, and the preferred value of N is 6. The approximate solutions _Sh and _Uh are used to replace S and U in the coupled equations (2), and the left and right sides of the equations are multiplied by the basis functions. Then, the integral by parts is used in the entire unit Ω and Gauss's theorem is applied to obtain the discrete integral form of DG-NCCR:

The process in Figure 1 is basically centered around solving equation (3) u _i (t) and s _i (t).

S01: Use commercial mesh generation software, such as Gambit or ICEM, to read the geometry file of the flow field to be solved, and then generate a structured or unstructured mesh in the flow area. Since the generation of a two-dimensional unstructured mesh is simple, and the unstructured mesh will be converted into a structured mesh in the solution of the present invention, the mesh generated in this step is a general two-dimensional triangular mesh. As shown in Figure 2, it is a mesh generated for the flow domain of a cylindrical aircraft in flight. As can be seen from Figure 2, the generated two-dimensional meshes are all ordinary triangles.

S02: According to the grid diagram shown in FIG2 , the geometric information of each grid is read, mainly the node coordinates of each grid unit.

S03: In order to simplify the calculation, the basis function adopts the Dubiner basis function, and its calculation needs to be performed under the standard grid. Therefore, it is necessary to convert each grid unit in Figure 2 into a standard triangle unit and a standard quadrilateral unit through step S03 to simplify the integral calculation. Specifically, for an ordinary unstructured triangle unit Ω, the coordinates of its three vertices 1, 2, and 3 are (x ₁ , y ₁ ), (x ₂ , y ₂ ), and (x ₃ , y ₃ ), respectively. The area of the triangle can be calculated based on the vertex coordinates:

The horizontal and vertical coordinates of the centroid of the triangular unit are:

With the above preparation, the original triangular unit Ω is converted to the standard triangular unit T, as shown in Figure 3. The relationship between the coordinate system of T and the Cartesian coordinate system is:

Ω→T:

Since the integration in the entire unit is performed in the standard square unit R = {(a, b), -1 ≤ a, b ≤ 1}, the detailed coordinate transformation is shown in Figure 4. The transformation relationship between the coordinate system where T is located and the coordinate system where R is located is:

T→R:

和

Dubiner基函数在该类型积分上的求解是在标准三角单元T内进行的，其计算过程需要用到两个坐标系之间的雅克比行列式：For the boundary integral of equation (3)

and

The solution of the Dubiner basis function on this type of integral is performed in the standard triangular unit T, and the calculation process requires the Jacobian determinant between the two coordinate systems:

For the remaining unit surface integrals, they are performed in the standard square unit, and the corresponding Jacobian determinant is:

S04: Constrain the equation group (3) by giving the boundary conditions of the actual rarefied gas flow. Set the outer boundary of the flow field as the pressure far field boundary, calculate the velocity by the given Mach number of the flow, and the other parameters are consistent with the surrounding atmospheric parameters. Set the wall surface as the Langmuir slip boundary. This boundary condition is based on the Langmuir adsorption isotherm and takes into account the interaction between gas molecules and the wall surface, thereby determining the velocity and temperature on the wall surface. For the initial conditions, calculate according to the boundary conditions, and assign 0 to the parameters that cannot be calculated.

S05: The expression of the surface integral can be obtained through the grid coordinate conversion formula and Jacobian determinant given above:

_Ai in the above formula is the integration coefficient. Similarly, the integration formula for integrals in other grids can be obtained.

S06: For boundary integrals, the solution of the present invention uses a numerical flux function to solve the solution on the boundary. There are many flux formats. In order to make the density and pressure in the calculation always positive, the present invention uses a always positive flux format to calculate the flux without viscosity term in the conservation equation:

here,

的求积公式(其他边界积分采用类似的形式求解)：The superscripts - and + represent the values inside and outside the unit interface, respectively. That is, the solution value of the calculation unit Ω ^k at a point on the boundary and the solution value of the point in the adjacent grid Ω ^k+1 are recorded as U ^- and U ⁺ respectively. The above flux expression can be obtained

The integral formula of (other boundary integrals are solved in a similar way):

S07: After solving the integral term of the discrete equation (3) in the above steps, its semi-discrete form can be obtained:

为质量矩阵，由于基函数是正交的，故质量矩阵是对角的且非常容易得到其逆矩阵。等号右边的R_U(U,S),R_S(S)矩阵为离散方程(3)除时间项以外其他积分的数值积分值。本发明方案采用3阶TVD-RK对上式进行时间离散：In the formula,

is the mass matrix. Since the basis functions are orthogonal, the mass matrix is diagonal and its inverse matrix is very easy to obtain. The R _U (U, S), R _S (S) matrices on the right side of the equal sign are the numerical integral values of the discrete equation (3) except the time term. The solution of the present invention uses the third-order TVD-RK to perform time discretization on the above equation:

in

Similarly, the time differential equation of _si is also discretized in the above form. The time step Δt can be calculated by the following formula:

x_n-1即为公式(5)中求解的

通过不断进行公式(5)的迭代，直至x_n-x_n-1<ε便完成了二维稀薄流场的计算。这里x_n不一定为所有求解变量，也可以选取密度和x、y两个方向的速度作为收敛判定条件。Here, CFL is the CFL number, which satisfies the condition CFL≤1. The _xn in step S07 is the solution of formula (5).

x _n-1 is the solution to formula (5)

By continuously iterating formula (5), the calculation of the two-dimensional rarefied flow field is completed until _xn - _xn-1 <ε. Here, _xn is not necessarily all the variables to be solved, and the density and the speed in the x and y directions can also be selected as the convergence judgment conditions.

S08: Since the DG method will produce non-physical numerical oscillations near discontinuous solutions (such as shock waves), the calculated pressure and density will appear negative, thus causing the results to diverge. In order to solve this problem, the present invention uses a limiter technology to capture the "problem unit" (i.e., the unit near the discontinuity) and reconstruct its solution to suppress the occurrence of non-physical oscillations. Specifically, the limiter used in the present invention is the TVB oblique limiter, the core of which is the introduction of rmin mod, thereby avoiding the accuracy at the critical point from being weakened to the first-order accuracy. The specific expression of the rmin mod function is:

Among them, the expression of the function sign(x) is:

Among them, ΔL is a parameter for judging whether restriction is needed, which is related to the geometric parameters and solutions of the grid and its surrounding grids. In the present invention, its expression is:

分别为求解值在单元内及与i边界相邻单元的平均值。Δx_i,Δy_i为单元顶点坐标差值矩阵Δx,Δy中的行元素，其表达式为：Among them, e is the total number of triangle unit edges,

are the average values of the solution values in the cell and the cells adjacent to the i boundary, respectively. _Δxi , _Δyi are the row elements in the cell vertex coordinate difference matrix Δx, Δy, and their expressions are:

The expressions of Δx _c,i ,Δy _c,i are:

Δx _c,i =x _c,i -x _c ,Δy _c,i =y _c,i -y _c

(x _c,i ,y _c,i ) are the centroid coordinates of the adjacent grid on the i boundary.

For the second-order case (the number of basis functions is 6), the average value of the solution in any unit is expressed as:

Similarly, the average values of the first-order and second-order derivatives of U _h in the unit can be calculated:

的求解方程组：The restricted solution coefficients can be established by using the rmin mod function and the average expression of the solution values in the element K(i,j) and its surrounding elements.

The solution of the equation system:

In the above formula,

Through the above formula, we can get the restricted approximate local solution of the "problem unit":

In this way, _xn can be updated with limited discontinuity points.

The embodiments described above are only descriptions of the preferred modes of the present invention and are not intended to limit the scope of the present invention. Without departing from the design spirit of the present invention, various deformations, modifications, and substitutions made to the technical solutions of the present invention by ordinary technicians in this field should all fall within the protection scope determined by the claims of the present invention.

Claims (10)

1. A numerical calculation method for two-dimensional rarefied gas flow, characterized in that: the DG numerical method is introduced into the NCCR equation to perform numerical discrete calculations, and a limiter is introduced to perform numerical discontinuity detection and limitation. 2. The numerical calculation method according to claim 1 is characterized in that: the NCCR equation is a dimensionless form equation obtained by simplifying dimensionless parameters and similarity criterion numbers. 3. The numerical calculation method according to claim 2, characterized in that: the global solutions _S and _U in the local unit Ω are expressed by approximate solutions Sh and Uh,

Where N is the number of basis functions;

is the basis function,

The DG-NCCR discrete equation is:

4 . The numerical calculation method according to claim 1 , wherein the limiter is a TVB slope limiter. 5. The numerical calculation method according to claim 3, characterized in that it comprises the following steps: S01. Meshing the two-dimensional solution flow field through a mesh program and generating a mesh file; S02, read the grid file and record the grid node coordinates; S03. Obtain standard square grid units according to the node coordinates; S04, given pressure far-field boundary conditions, Langmuir wall slip boundary conditions; S05, initializing calculation conditions, determining the time step according to CFL conditions; assigning initial values x ₀ , x _n-1 = x ₀ to flow field parameters according to boundary conditions; S06, substituting x _n-1 into the numerical integration term in the DG-NCCR discrete equation for calculation; S07, using different flux calculation formats to calculate the flux in the discrete equation; S08. According to the calculation results of steps S06 and S07, the discrete equation is solved by using a third-order TVD-RK to obtain flow field parameters x _n . 6. The numerical calculation method according to claim 5 is characterized in that: in the step S01, the generated two-dimensional grid is a triangular grid unit; in the step S03, the triangular grid unit in step S01 is first converted into a standard triangular unit, and then converted into a standard square grid unit. 7. The numerical calculation method according to claim 5, characterized in that it comprises the following steps: S09, traverse all solution results, find the "problem unit" and use the limiter to limit the solution parameters in the unit, and update x _n ; S10, calculating the iteration error x _n -x _n-1 , and determining the magnitude of the iteration error and the target calculation accuracy ε: When x _n -x _n-1 <ε, go to step S11; When _xn - _xn- 1≥ε, let _xn-1 ＝ _xn , and enter step S06 to perform loop calculation until the calculation result meets the target calculation accuracy requirement ε; S11, the loop calculation ends. 8. Application of the numerical calculation method for two-dimensional rarefied gas flow as claimed in any one of claims 1 to 7 in the flight of near-space aircraft. 9. An electronic device, characterized in that it comprises: one or more processors; A storage device for storing one or more programs; When the one or more programs are executed by the one or more processors, the one or more processors implement the numerical calculation method for two-dimensional rarefied gas flow as described in any one of claims 1 to 7. 10. A computer-readable medium storing a computer program, wherein the computer program, when executed by a processor, implements the numerical calculation method for two-dimensional rarefied gas flow according to any one of claims 1 to 7.

Images