CN115952748A - Numerical calculation method for two-dimensional rarefied gas flow - Google Patents
Numerical calculation method for two-dimensional rarefied gas flow Download PDFInfo
- Publication number
- CN115952748A CN115952748A CN202211703473.9A CN202211703473A CN115952748A CN 115952748 A CN115952748 A CN 115952748A CN 202211703473 A CN202211703473 A CN 202211703473A CN 115952748 A CN115952748 A CN 115952748A
- Authority
- CN
- China
- Prior art keywords
- calculation
- equation
- numerical
- dimensional
- flow
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000004364 calculation method Methods 0.000 title claims abstract description 81
- 238000000034 method Methods 0.000 claims abstract description 45
- 238000001514 detection method Methods 0.000 claims abstract description 6
- 230000006870 function Effects 0.000 claims description 20
- 230000004907 flux Effects 0.000 claims description 13
- 238000004590 computer program Methods 0.000 claims description 4
- 230000000670 limiting effect Effects 0.000 claims description 3
- 125000004122 cyclic group Chemical group 0.000 claims 1
- 230000035939 shock Effects 0.000 abstract description 9
- 230000008901 benefit Effects 0.000 abstract description 5
- 238000000329 molecular dynamics simulation Methods 0.000 abstract description 5
- 230000007704 transition Effects 0.000 abstract description 5
- 239000007789 gas Substances 0.000 description 27
- 230000014509 gene expression Effects 0.000 description 12
- 230000010354 integration Effects 0.000 description 9
- 239000012530 fluid Substances 0.000 description 7
- 238000010586 diagram Methods 0.000 description 5
- 230000008878 coupling Effects 0.000 description 4
- 238000010168 coupling process Methods 0.000 description 4
- 238000005859 coupling reaction Methods 0.000 description 4
- 238000005206 flow analysis Methods 0.000 description 4
- 239000011159 matrix material Substances 0.000 description 4
- 230000008569 process Effects 0.000 description 4
- 230000009466 transformation Effects 0.000 description 4
- 238000005315 distribution function Methods 0.000 description 3
- 239000000853 adhesive Substances 0.000 description 2
- 230000001070 adhesive effect Effects 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 238000011438 discrete method Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- FKLFBQCQQYDUAM-UHFFFAOYSA-N fenpiclonil Chemical compound ClC1=CC=CC(C=2C(=CNC=2)C#N)=C1Cl FKLFBQCQQYDUAM-UHFFFAOYSA-N 0.000 description 2
- 230000003993 interaction Effects 0.000 description 2
- 230000010355 oscillation Effects 0.000 description 2
- 238000001179 sorption measurement Methods 0.000 description 2
- 235000015842 Hesperis Nutrition 0.000 description 1
- 235000012633 Iberis amara Nutrition 0.000 description 1
- 238000000342 Monte Carlo simulation Methods 0.000 description 1
- 238000010521 absorption reaction Methods 0.000 description 1
- 230000004075 alteration Effects 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 125000004432 carbon atom Chemical group C* 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 230000006835 compression Effects 0.000 description 1
- 238000007906 compression Methods 0.000 description 1
- 238000009792 diffusion process Methods 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 230000005484 gravity Effects 0.000 description 1
- 230000002401 inhibitory effect Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000002360 preparation method Methods 0.000 description 1
- 238000007789 sealing Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000003313 weakening effect Effects 0.000 description 1
Images
Classifications
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
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
Technical Field
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.
Background
Conventional fluid mechanics equations have been studied focused on continuous regions, i.e., the mean free path of fluid molecules is much smaller than the flow characteristic dimension, and the navier-stokes-fourier equation (NSF equation) is a typical representation of this type of equation. The method is characterized in that three flow conservation equations of continuity, momentum and energy are closed by coupling Stokes assumption of the constitutive relation of Newton viscous stress and Fourier heat conduction. In nearly two centuries, the NSF equation has enjoyed great success in the continuous flow field, driving the development of hydrodynamics. However, the continuity assumption fails when the mean free path of the molecule gradually increases, i.e., the gas becomes thinner and less dense, at which time it is necessarily inappropriate to continue to employ NSF equations based on a continuous medium or equilibrium state (when the system's macroscopic thermal observations no longer change over time, at which point the system is in equilibrium).
In recent years, with the continuous progress of science and human exploration of unknown fields, the field of thin gases and its engineering applications, such as the prediction of the aerodynamic force of an aircraft flying in the atmosphere where the air is thin (several tens of kilometers of high altitude) have been of increasing interest. The gas density at this condition is very low, roughly 10 -7 To 10 -10 kg/m 3 At this time, the distance that a molecule has to travel twice due to the larger molecular distance is larger, that is, the mean free path of the molecule is larger, and the collision effect of the molecule gradually appears, so the traditional macroscopic equation, NSF equation, is no longer applicable, and the Boltzmann equation with collision terms plays an important role in solving the rarefied flow.
According to Knudsen K n The size of the number (important parameters for measuring the gas rareness degree and the macroscopic model correctness degree) can be divided into four types:
continuous flow domain (Kn < 0.001), conventional euler's equation or NS equation is applicable;
in the slip flow field (0.001-Kn-0.1), NS equation is established in the main flow outside the boundary, but the boundary conditions of speed slip and temperature jump are considered at the boundary;
in the transition flow field (0.1-Kn-10), the NS equation is no longer true;
free molecular flow (10 > Kn) directly using molecular kinematics methods such as DSMC.
Among the various methods for solving the flow problems described above, the NS equation and the molecular kinematics theory have succeeded in the continuous flow and free molecular flow fields, respectively, but for the slip flow domain and the transition flow domain, there is no unified and effective control equation for solving the flow domain of the region, for example, the gas flow in the critical space layer (altitude 20-100 km) which has been vigorously developed in recent years is exactly in the two flow domains, so it is extremely important to develop a gas dynamics theory applicable to both continuous flow and thin flow.
To study lean gas flow, a number of mathematical physical equations have been proposed, which can be basically divided into three categories: the molecular dynamic model solved from the microscopic angle comprises a direct simulation Monte Carlo method (DSMC), a DSMC information storage method (DSMC-IP) and a Lattice Boltzmann Method (LBM); there are also methods for directly performing discrete modeling solution on the distribution function in the Boltzmann equation, such as the BGK model equation; and substituting the macroscopic statistical expression into a microscopic Boltzmann equation to further obtain a macroscopic quantity conservation equation and an evolution equation, wherein the representative methods are a Grad method, a Chapman-Enskog method and a Eu method.
For the molecular dynamics model, since it uses simulated molecules to simulate the movement of real molecules and the collision of gas flow, it costs a lot of computing resources at low Kn number (large concentration of molecules). For the BGK model equation, the collision term in the Boltzmann equation is simplified, so that the BGK model equation can obtain an accurate result only in an equilibrium state or a near-equilibrium state, and the BGK model equation also has inaccuracy in the calculation of the transport coefficient. The prandtl (Pr) number for the monatomic gas, as calculated by the BGK model, is 1, rather than the correct value of 2/3. In the non-NS fluid control equation method evolved by Boltzmann equation, the R-13 equation based on the moment method proposed by Grad and the Burnet equation based on Chapman-Enskog expansion are not considered to satisfy the entropy condition of the thermodynamic law, i.e. the definite solution of Boltzmann equation is difficult to satisfy, so the development of the two methods is limited to a certain extent. Different from the two methods, the Eu method strictly meets the entropy condition, and starts from Boltzmann theorem, grasps the characteristic of entropy increase from the unbalanced state to the balanced state, constructs an exponential-form distribution function, obtains an entropy increase dissipation model from the unbalanced state to the balanced state, converges to a Rayleigh-Onsager dissipation function when the model approaches to the balanced state, and processes the collision item by constructing a unified nonlinear entropy increase model from the unbalanced state to the balanced state on the basis of the Rayleigh-Onsager dissipation function. Therefore, three conservation equations and constitutive equations of viscous stress, heat flux and the like are derived through the Boltzmann equation, and the conservation equations and the non-conservation equations are coupled together to complete the sealing of the gas dynamic equation to form the Eu unified fluid equation. Thereafter, myong focuses on the higher order terms of the constitutive equation based on Eu and constructs the Nonlinear Coupling Constitutive Relation (NCCR) based thereon.
Therefore, if a numerical calculation method for two-dimensional lean gas flow can be provided, the method has a good application prospect.
Disclosure of Invention
In order to solve the technical problems, the invention provides a numerical calculation method suitable for the flowing of the rarefied gas based on a two-dimensional NCCR equation derived by an Eu method, in particular to the non-equilibrium heat flow analysis of hypersonic flight with shock wave phenomenon.
In order to achieve the purpose of the invention, the technical scheme adopted by the invention is as follows: a numerical calculation method for two-dimensional rarefied gas flow introduces a DG numerical method into an NCCR equation to carry out numerical discrete calculation, and introduces a limiter to carry out numerical discontinuity detection and limitation.
Preferably: the NCCR equation is a dimensionless formal equation obtained by simplifying dimensionless parameters and similarity criterion numbers.
Preferably: by approximating the solution S h 、U h To express the global solutions S and U within the local unit omega,
in the formula, N is the number of basis functions;
is a base function->
The DG-NCCR discrete equation is as follows:
preferably: the limiter is a TVB slope limiter.
Preferably: the method comprises the following steps:
s01, carrying out grid division on the two-dimensional solution flow field through a grid program, and generating a grid file;
s02, reading a grid file, and recording grid node coordinates;
s03, obtaining standard square grid units according to the node coordinates;
s04, giving a pressure far-field boundary condition and a Langmuir wall surface sliding boundary condition;
s05, initializing calculation conditions, and determining time step length according to CFL conditions; assigning an initial value x to the flow field parameter according to the boundary condition 0 ，x n-1 ＝x 0 ；
S06, mixing x n-1 Substituting into a numerical integral term in a DG-NCCR discrete equation for calculation;
s07, calculating the flux in the discrete equation by adopting different flux calculation formats;
s08, according to the calculation results of the steps S06 and S07, adopting third-order TVD-RK to perform dispersionSolving an equation to obtain a flow field parameter x n 。
Preferably: in the step S01, the generated two-dimensional mesh is a triangular mesh unit; in step S03, the triangular mesh unit in step S01 is converted into a standard triangular mesh unit, and then converted into a standard square mesh unit.
Preferably: further comprising:
s09, traversing all the solving results, searching a problem unit, limiting the solving parameters in the problem unit by adopting a limiter, and updating x n ；
S10, calculating an iteration error x n -x n-1 And judging the magnitude of the iteration error and the target calculation precision epsilon:
when x is n -x n-1 <ε, proceed to step S11;
when x is n -x n-1 Not less than epsilon, let x n-1 ＝x n Step S06 is carried out for circular calculation until the calculation result meets the target calculation precision requirement epsilon;
and S11, finishing the loop calculation.
Correspondingly: the application of the numerical calculation method for the two-dimensional rarefied gas flow in the flight of the adjacent space aircraft.
Correspondingly: an electronic device, characterized in that: the method comprises the following steps:
one or more processors;
storage means 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 lean gas flow.
Correspondingly: a computer readable medium, said readable medium storing a computer program, characterized in that: the computer program, when executed by a processor, implements a numerical calculation method for two-dimensional lean gas flow.
Compared with the prior art, the invention has the following beneficial effects:
the DG method is introduced into an NCCR equation to carry out numerical discrete calculation, and meanwhile, the limiter is introduced into the numerical method to carry out numerical discontinuous detection and limitation, so that the shock wave problem of the hypersonic speed in a rarefied state is well solved. 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.
Drawings
FIG. 1 is a block diagram of a flow chart of a two-dimensional DG-NCCR discrete equation calculation method of the present invention;
FIG. 2 is a two-dimensional cylindrical flow-around unstructured grid diagram according to an embodiment of the present invention;
FIG. 3 is a schematic diagram of coordinate transformation between an arbitrary triangle unit Ω and a standard triangle unit T according to the present invention;
FIG. 4 is a schematic diagram of coordinate transformation between a standard triangle cell T and a standard square cell R according to the present invention;
FIG. 5 is a schematic diagram of a two-dimensional restriction unit K and its neighboring units according to the present invention.
Detailed Description
The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. Unless otherwise specified, the technical means used in the examples are conventional means well known to those skilled in the art.
The invention provides a numerical calculation method suitable for two-dimensional rarefied gas flow analysis, which has the core idea that: and introducing a DG numerical method into an NCCR equation to carry out numerical discrete calculation, and introducing a limiter to carry out numerical discontinuity detection and limitation. The fluid transport equation used in the method is derived from the Boltzmann equation, which is generally known as the full-flow-domain gas dynamics control equation. Therefore, the invention aims to break through the disadvantages that the traditional NS equation cannot solve discontinuous flow and the molecular dynamic model has large calculated amount in the transition flow and the slip flow, and provides a thermal aerodynamic analysis method suitable for the flight process of the near space vehicle.
The numerical calculation method of the invention adopts a DG (interrupted Galerkin) numerical method in the process of dispersing a two-dimensional NCCR equation, combines the advantages of the traditional Finite Element (FEM) and the Finite Volume (FVM), and specifically comprises the following steps: (1) The method has high adaptability to unstructured grids and is suitable for complex configurations; (2) When the grid is encrypted or thickened, continuity is not considered, and the self-adaptive technology of grid encryption can be easily realized; (3) Each computing unit is independent, which has the advantages of easy coding and parallelism; (4) Due to the operation of discontinuity approximation in the cell, the DG method is well suited for solving numerical discontinuity problems. The DG-NCCR discrete equation has good solving advantage for the problem of the hypersonic shock wave in the near space, the shock wave is a compression wave with suddenly changed airflow parameters in the hypersonic streaming process, and the shock wave widely exists in hypersonic flight of shells, rockets, airplanes and the like, and the hypersonic problem in the near space is more difficult to solve and analyze because the air in the near space is thin and the atmospheric temperature and pressure are lower and the air environment is lower.
To better explain the numerical calculation method of the present invention, the source of the two-dimensional NCCR equation is briefly given first. For a two-dimensional Boltzmann equation, by introducing an exponential distribution function and a corresponding Rayleigh-Onsager dissipation function, a nonlinear coupling constitutive equation and NCCR constitutive relation meeting the principle of entropy increase is obtained:
wherein,
in the formula, Λ k Is a molecular collision dissipation term, Z k The kinetic energy term that causes the streamline effect of the fluid for molecular diffusion. U is a conservative variable, F inv (U) is no adhesive, F vis (U) is a viscosity term, Π is viscosityStress, Q is thermal conduction, and I is the unit tensor.
q(k)＝sinhk/k，k 2 Is a Rayleigh-Onsager dissipation function, gamma' = (5-3 gamma)/2 (gamma is specific heat ratio), eta b For adding stress viscosity coefficient, C p The specific heat capacity at constant pressure is shown.
p is pressure, ρ is density, T is temperature, Δ is volume added normal stress, and E is energy per unit mass.
Here, both the velocity and divergence are two-dimensional cases. The first equation in the above equations is the conservation equation of mass, momentum and energy, and the second equation is the constitutive equation, which gives the expressions of viscous stress pi and heat flux Q for enclosing the conservation equation set.
To simplify the calculation, the above equation is dimensionless:
wherein the subscript r denotes the reference state, f b Expressing the ratio of bulk viscosity to shear viscosity, which can be measured using a sonic absorption experiment, f for a monoatomic molecule b =0, the additional viscous stress is zero. t is time, x is length, L is flow field characteristic length, η is hydrodynamic viscosity coefficient, λ is thermal conductivity coefficient, and u is velocity.
The following similarity criteria numbers are defined:
and N delta represents the magnitude of viscous stress relative to static pressure, the magnitude of the viscous stress can measure the degree of the system far away from an equilibrium state, ma is Mach number, re is Reynolds number, ec is a function of the Mach number Ma, and Pr is a Prandtl number. Then the dimensionless vector form of the system of equations can be obtained by the above dimensionless parameters and the similarity criterion numbers:
wherein,
the solution of equation set (2) is a numerical calculation solution for the present solution involving two-dimensional lean gas flow. The temperature field, the speed field, the pressure field and the density field of the two-dimensional rarefaction can be solved by solving the flow parameters of the equation, such as temperature, speed, pressure, density and the like, so that the thermal flow analysis of rarefied hypersonic flight is realized.
A numerical calculation method for two-dimensional rarefied gas flow is used for solving a DG-NCCR discrete equation, and comprises the following steps:
and S01, carrying out grid division on the two-dimensional solution flow field through a grid program, and generating a grid file. Commercial grid generation software can be selected to generate the grid file, and the grid generation software can be Gambit or ICEM. The mesh generated in this step is typically a two-dimensional triangular mesh.
And S02, reading the grid file and recording the grid node coordinate parameters. The grid file is read by the two-dimensional DG-NCCR program, and the numerical calculation method for the two-dimensional rarefied gas flow disclosed by the invention is written in the two-dimensional DG-NCCR program.
And S03, obtaining a standard square grid unit according to the node coordinates. Specifically, the triangular mesh unit in step S01 is converted into a standard triangular mesh unit, and then converted into a standard square mesh unit.
S04, setting a pressure far-field boundary condition, calculating the speed through a given flow-around Mach number, and enabling the rest parameters to be consistent with the surrounding atmospheric parameters. The wall is set to Langmuir slip boundary conditions that determine velocity and temperature on the wall based on Langmuir adsorption isotherms and taking into account the interaction of gas molecules with the wall.
S05, carrying out iterative calculation: initializing calculation conditions, such as the maximum iteration step number M and the calculation precision requirement epsilon, and determining the time step length according to the CFL condition; assigning an initial value x to the flow field parameter according to the boundary condition 0 ，x n-1 ＝x 0 . The CFL condition is an important concept in stability and convergence analysis in finite difference and finite volume methods. n denotes the nth calculation, n =1,2,3, … ….
S06, mixing x n-1 And substituting into a numerical integral term in a DG-NCCR discrete equation for calculation.
And S07, calculating the flux in the DG-NCCR discrete equation by adopting different flux calculation formats.
S08, solving the DG-NCCR discrete equation by adopting a third-order TVD-RK according to the calculation results of the steps S06 and S07 to obtain a flow field parameter result x of the current solving step n 。
S09, traversing all the solving results, searching a 'problem unit', limiting solving parameters in the problem unit by adopting a limiter, and updating x n . The limiter is a TVB slope limiter.
S10, calculating iteration error x n -x n-1 And judging the magnitude of the iteration error and the target calculation precision epsilon: when x is n -x n-1 <ε, proceed to step S11; when x is n -x n-1 Not less than epsilon, let x n-1 ＝x n And step S06 is carried out for circular calculation until the calculation result meets the target calculation precision requirement epsilon.
And S11, finishing loop calculation. According to x of each unit n And performing inverse dimensionless solving on the original value of the flow field parameter, and substituting the result data into post-processing software to obtain a flow field cloud chart and a flow chart, wherein the flow field cloud chart and the flow chart are used for performing flow analysis by using related fluid post-processing software such as techniclot, origin and the like.
Because NCCR is highly nonlinear and coupled, the NCCR is difficult to obtain a good solution in a common discrete method, the numerical calculation method of the invention introduces a DG method into an NCCR equation to carry out numerical discrete calculation, and simultaneously introduces a limiter into the numerical method to carry out numerical discontinuity detection and limitation, so that the invention has a good solution to the shock wave problem of hypersonic velocity in a rarefied state.
Examples
The expression form of the DG-NCCR discrete equation is determined first.
By approximating the solution S h 、U h To approximate the global solutions S and U within the local unit omega,
in the formula,
is a basis function which is used for numerical integration in the DG discrete method due to the following properties in integration:
it is to be noted that each sheetThe base function form in the element is the same, u i (t) and s i (t) are respectively the in-cell approximate solutions S h And U h Knowing their values, the approximate solution of the cell is obtained.
N is the number of the basis functions, the size of the basis functions is restricted by the calculation order I, and the specific relationship is as follows:
in order to ensure the calculation accuracy and save the calculation resources, the order I is preferably 2, and the order N is preferably 6. By approximating the solution S h And U h Instead of coupling S and U in equation set (2), multiplying the basis functions at the left and right ends of the equation, then using fractional integration over the entire unit Ω and applying the Gaussian theorem, one can get a discrete integral form of DG-NCCR:
the flow in FIG. 1 is basically all around the above formula (3) u i (t) and s i And (t) solving and expanding.
S01: commercial grid generation software such as Gambit or ICEM is adopted to read a flow field geometric file to be solved, and then a structural grid or an unstructured grid is generated in a flow area. As shown in fig. 2, the mesh generated for the flow field in the flight of a cylindrical aircraft, it can be seen from fig. 2 that the two-dimensional meshes generated are all generally triangular.
S02: the geometric information of each grid is read according to the grid graph shown in fig. 2, mainly the node coordinates of each grid cell.
S03: in order to simplify the calculation, the basis function is the Dubiner basis function, and the calculation thereof needs to be performed under the standard mesh, so that it is necessary to perform step S03 for each of fig. 2The grid cells are converted into standard triangular cells and standard quadrilateral cells, thereby simplifying integral calculation. Specifically, for a common unstructured triangle unit Ω, its three vertices 1,2,3 have coordinates of (x) 1 ,y 1 ),(x 2 ,y 2 ),(x 3 ,y 3 ) The area of the triangle can be calculated from the vertex coordinates:
the horizontal and vertical coordinates of the gravity center of the triangular unit are respectively as follows:
with the above preparation, the original triangle unit Ω is converted into the standard triangle unit T, and as shown in fig. 3, the relationship between the coordinate system of T and the cartesian coordinate system is:
Ω→T:
since the integral in the whole unit is performed in the standard square unit R = { (a, b), -1 ≦ a, b ≦ 1}, the coordinate transformation is detailed as shown in FIG. 4, and the transformation relationship between the coordinate system of T and the coordinate system of R is as follows:
T→R:
boundary integration for equation (3)
And &>
Dubiner basis letterThe solution of the number on this type of integral is performed in a standard trigonometric unit T, the calculation of which requires the use of the jacobian determinant between two coordinate systems: />
And for the rest unit area integration, the integration is carried out in a standard square unit, and the corresponding Jacobian is as follows:
s04: the constraint equation set (3) sets the outer flow field boundary to the pressure far field boundary by giving the boundary condition when the actual lean gas flows, calculates the velocity by giving the flow around mach number, and the rest of the parameters are consistent with the ambient atmospheric parameters. The wall is set to a Langmuir slip boundary, which boundary condition is based on Langmuir adsorption isotherms and takes into account the interaction of gas molecules with the wall, thereby determining the velocity and temperature on the wall. For the initial conditions, the calculation is performed according to the boundary conditions, and for the parameters which cannot be calculated, the values are assigned to 0.
S05: the expression of the surface integral can be obtained through the grid coordinate conversion formula given above and the Jacobian:
a in the above formula i For the multiplication coefficient, the multiplication formula of the integration in other grids can be obtained by the same method.
S06: for boundary integration, the scheme of the invention adopts a numerical flux function to solve the solution on the boundary, a plurality of flux formats exist, and in order to make the density and the pressure in the calculation constant, the invention adopts a constant flux format to calculate the flux without an adhesive term in a conservation equation:
here, ,
the superscripts-and + represent the values inside and outside the cell interface, respectively, i.e. the calculation cell Ω k Solution value for a point on the boundary and the point in the adjacent grid Ω k+1 The solution values in (1) are respectively marked as U - And U + . The flux expression can be obtained
The product formula (other boundary integrals are solved in a similar fashion):
s07: the semi-discrete form can be obtained by solving the integral term of the discrete equation (3) in the previous step:
in the formula,
for the quality matrix, since the basis functions are orthogonal, the quality matrix is diagonal and its inverse is very easy to obtain. R on the right of equal sign U (U,S),R S The (S) matrix is the numerical integration value of the discrete equation (3) integrated in addition to the time term. The scheme of the invention adopts 3-order TVD-RK to carry out time dispersion on the formula:
wherein
In the same way, for s i The time differential equation of (a) is also time-discrete in the above-described form. The time step Δ t can be calculated by:
here, CFL is the number of CFLs, which satisfies the condition CFL ≦ 1. X in step S07 n I.e. solved in equation (5)
x n-1 I.e., based on: (5)>
By iterating equation (5) until x n -x n-1 <Epsilon completes the calculation of the two-dimensional lean flow field. Where x is n Not all solution variables are necessary, and the density and the velocities in the x and y directions may be selected as the convergence determination conditions.
S08: in order to solve the problem that the DG method can generate non-physical numerical value oscillation near an intermittent solution (such as shock waves) so that the calculated pressure and density have negative values and the result is diffused, the method adopts a limiter technology to capture a problem unit (namely a unit near the interruption) and reconstructs the solution of the problem unit, thereby inhibiting the non-physical oscillation. Specifically, the limiter adopted by the invention is a TVB slope limiter, the core of the TVB slope limiter is to introduce rmin mod, so that the precision at a critical point is prevented from weakening to first-order precision, and the specific expression of the rmin mod function is as follows:
wherein, the expression of the function sign (x) is as follows:
wherein, Δ L is a parameter for judging whether limitation is needed, and is related to the solution of the geometric parameters of the grid and the surrounding grids thereof, and the expression thereof in the invention is as follows:
wherein e is the total number of the triangle unit sides respectively,
the average values of the solution values within the cell and the cell adjacent to the i boundary are respectively. Δ x i ,Δy i The expression of the row element in the unit vertex coordinate difference matrix Δ x, Δ y is:
and Δ x c,i ,Δy c,i The expression of (a) is:
Δx c,i ＝x c,i -x c ,Δy c,i ＝y c,i -y c
(x c,i ,y c,i ) Is the barycentric coordinates of the i-boundary neighboring grid.
For the average value of solution values in any unit under the second order condition (the number of basis functions is 6), the expression is as follows:
similarly, U can be found h Mean of the first and second derivatives within the cell:
the solution coefficient for the constraint can be established by the rmin mod function and the mean expression of the solution values in the cell K (i, j) and its surrounding cells
The system of solved equations of (c):
in the above formula, the first and second carbon atoms are,
through the above equation, a constrained approximate local solution of the "problem unit" can be obtained:
thus, can be aligned with x n And performing limit updating of the discontinuity points.
The above-described embodiments are merely illustrative of the preferred embodiments of the present invention, and do not limit the scope of the present invention, and various changes, modifications, alterations, and substitutions which may be made by those skilled in the art without departing from the spirit of the present invention shall fall within the protection scope defined by the claims of the present invention.
Claims (10)
1. A numerical calculation method for two-dimensional lean gas flow, characterized by: and introducing a DG numerical method into an NCCR equation to carry out numerical discrete calculation, and introducing a limiter to carry out numerical discontinuity detection and limitation.
2. A numerical calculation method according to claim 1, characterized in that: the NCCR equation is a dimensionless formal equation obtained by simplifying dimensionless parameters and similarity criterion numbers.
3. A numerical calculation method according to claim 2, characterized in that: by approximating the solution S h 、U h To express the global solutions S and U within the local unit omega,
in the formula, N is the number of basis functions;
is a base function->
The DG-NCCR discrete equation is as follows:
4. a numerical calculation method according to claim 1, characterized in that: the limiter is a TVB slope limiter.
5. A numerical calculation method according to claim 3, characterized in that: the method comprises the following steps:
s01, carrying out grid division on the two-dimensional solution flow field through a grid program, and generating a grid file;
s02, reading a grid file, and recording grid node coordinates;
s03, obtaining standard square grid units according to the node coordinates;
s04, giving a pressure far-field boundary condition and a Langmuir wall surface sliding boundary condition;
s05, initializing calculation conditions, and determining time step length according to CFL conditions; assigning an initial value x to the flow field parameter according to the boundary condition 0 ，x n-1 ＝x 0 ；
S06, mixing x n-1 Substituting into a numerical integral term in a DG-NCCR discrete equation for calculation;
s07, calculating the flux in the discrete equation by adopting different flux calculation formats;
s08, solving the discrete equation by adopting a third-order TVD-RK according to the calculation results of the steps S06 and S07 to obtain a flow field parameter x n 。
6. The numerical calculation method according to claim 5, characterized in that: in the step S01, the generated two-dimensional mesh is a triangular mesh unit; in step S03, the triangular mesh unit in step S01 is converted into a standard triangular mesh unit, and then converted into a standard square mesh unit.
7. A numerical calculation method according to claim 5, characterized in that: the method comprises the following steps: further comprising:
s09, traversing all the solving results, searching a problem unit, limiting the solving parameters in the problem unit by adopting a limiter, and updating x n ；
S10, calculating iteration error x n -x n-1 And judging the magnitude of the iteration error and the target calculation precision epsilon:
when x is n -x n-1 <ε, proceed to step S11;
when x is n -x n-1 Not less than epsilon, let x n-1 ＝x n Step (c) of enteringS06, performing cyclic calculation until the calculation result meets the target calculation precision requirement epsilon;
and S11, finishing the loop calculation.
8. Use of a method according to any one of claims 1 to 7 for the numerical calculation of a two-dimensional rarefied gas flow in the flight of a near space vehicle.
9. An electronic device, characterized in that: the method comprises the following steps:
one or more processors;
storage means for storing one or more programs;
when executed by the one or more processors, cause the one or more processors to implement the numerical calculation method for two-dimensional lean gas flow of any one of claims 1 to 7.
10. A computer readable medium, said readable medium storing a computer program, characterized in that: the computer program, when executed by a processor, implements a numerical calculation method for two-dimensional lean gas flow as claimed in any one of claims 1 to 7.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211703473.9A CN115952748A (en) | 2022-12-29 | 2022-12-29 | Numerical calculation method for two-dimensional rarefied gas flow |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211703473.9A CN115952748A (en) | 2022-12-29 | 2022-12-29 | Numerical calculation method for two-dimensional rarefied gas flow |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN115952748A true CN115952748A (en) | 2023-04-11 |
Family
ID=87287223
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202211703473.9A Pending CN115952748A (en) | 2022-12-29 | 2022-12-29 | Numerical calculation method for two-dimensional rarefied gas flow |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115952748A (en) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116757126A (en) * | 2023-08-21 | 2023-09-15 | 中国空气动力研究与发展中心计算空气动力研究所 | Method for determining flow stability of lean gas based on gas dynamic theory |
| CN117933144A (en) * | 2024-03-22 | 2024-04-26 | 中国空气动力研究与发展中心超高速空气动力研究所 | Multiple grid method for solving grid flow field of complex topological structure |
-
2022
- 2022-12-29 CN CN202211703473.9A patent/CN115952748A/en active Pending
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116757126A (en) * | 2023-08-21 | 2023-09-15 | 中国空气动力研究与发展中心计算空气动力研究所 | Method for determining flow stability of lean gas based on gas dynamic theory |
| CN116757126B (en) * | 2023-08-21 | 2023-12-12 | 中国空气动力研究与发展中心计算空气动力研究所 | Method for determining flow stability of lean gas based on gas dynamic theory |
| CN117933144A (en) * | 2024-03-22 | 2024-04-26 | 中国空气动力研究与发展中心超高速空气动力研究所 | Multiple grid method for solving grid flow field of complex topological structure |
| CN117933144B (en) * | 2024-03-22 | 2024-05-28 | 中国空气动力研究与发展中心超高速空气动力研究所 | Multiple grid method for solving grid flow field of complex topological structure |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN115952748A (en) | Numerical calculation method for two-dimensional rarefied gas flow | |
| Singh et al. | General drag coefficient for flow over spherical particles | |
| Fei et al. | A benchmark study of kinetic models for shock waves | |
| Liu et al. | Three-dimensional high-order least square-based finite difference-finite volume method on unstructured grids | |
| Xiaoquan et al. | Robust implicit direct discontinuous Galerkin method for simulating the compressible turbulent flows | |
| Yuan et al. | Modified nonlinear coupled constitutive relations model for hypersonic nonequilibrium flows | |
| Chinnappan et al. | Transport dynamics of an ellipsoidal particle in free molecular gas flow regime | |
| Zakeri et al. | New chemical-DSMC method in numerical simulation of axisymmetric rarefied reactive flow | |
| Ching et al. | Sensitivity of hypersonic dusty flows to physical modeling of the particle phase | |
| Lobbia | Rapid supersonic/hypersonic aerodynamics analysis model for arbitrary geometries | |
| Zhan et al. | Discrete gas-kinetic scheme-based arbitrary Lagrangian–Eulerian method for moving boundary problems | |
| Zangeneh | Development of a new algorithm for modeling viscous transonic flow on unstructured grids at high Reynolds numbers | |
| Zhao et al. | Numerical simulation of lateral jet interaction with rarefied hypersonic flow over a two-dimensional blunt body | |
| Dobrov et al. | Simulation of high-temperature flowfield around hypersonic waverider using graphics processor units | |
| Hinkle et al. | Efficient Solution of Surface Erosion in Particle-Laden Hypersonic Flows | |
| Pekardan et al. | Rarefaction effects for transonic airfoil flows at low Reynolds numbers | |
| Tan et al. | Improved mode multigrid method for accelerating turbulence flows | |
| Li et al. | Novel hybrid hard sphere model for direct simulation Monte Carlo computations | |
| Zhao et al. | Computation of rarefied hypersonic flows using modified form of conventional Burnett equations | |
| Hepp et al. | A kinetic Fokker–Planck approach for modeling variable hard-sphere gas mixtures | |
| Jain et al. | Second moment closure modeling and direct numerical simulation of stratified shear layers | |
| Açıkgöz et al. | Dynamic mesh analyses of helicopter rotor–fuselage flow interaction in forward flight | |
| Cooper et al. | Implementation and verification of a surface recession module in a finite volume ablation solver | |
| Tensuda et al. | Application of a maximum-entropy-based 14-moment closure for multi-dimensional non-equilibrium flows | |
| Raeisi et al. | Numerical investigation of interaction of counter flow jet and hypersonic capsule flow via modified DsmcFoam |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |