A high-precision numerical simulation method suitable for hypersonic aircraft

Abstract

The application provides a high-precision numerical simulation method suitable for a hypersonic aircraft, first, a geometric model and a calculation domain of the hypersonic aircraft are established, and the calculation grid of the hypersonic aircraft is divided and read to obtain grid information and initial values according to inflow conditions; then, the control equation is discretized to obtain a semi-discrete format of a finite volume; then, the reconstruction values of both sides of the interface of the grid element are obtained according to the average of the grid element in the calculation domain of the hypersonic aircraft; the reconstruction values of all the Gauss points of the interface of the grid element in the hypersonic flow field are obtained by using the reconstruction values of both sides of the interface of all the grid elements, and are used for a multi-dimensional Riemann solver; the interface flux is obtained by using the obtained multi-dimensional Riemann solver; the residual is determined according to the interface flux, and time advancing is solved to obtain the final flow field of the hypersonic aircraft. The application can provide technical support for more accurate hypersonic numerical simulation tasks and hypersonic aircraft design work.

Description

Technical Field

[0001] This invention relates to the field of computational fluid dynamics, specifically to a high-precision numerical simulation method applicable to hypersonic vehicles. Background Technology

[0002] Hypersonic vehicles possess immense military and civilian value due to their high speed and strong penetration capabilities. For example, hypersonic cruise weapons offer greater lethality and accuracy compared to conventional weapons, while being less likely to be intercepted; hypersonic fighter jets boast high attack speeds and strong penetration capabilities, enabling them to efficiently complete global strike missions in a short time; and hypersonic air defense missiles can quickly reach interception airspace, shortening reaction time and improving interception efficiency. Therefore, the United States, Russia, the European Union, Japan, and other countries have invested heavily in research on hypersonic technologies, developing a series of research programs such as Hyper-X, HyFly, HIFiRE, HAWC, and HSSW.

[0003] Related studies have revealed that the aerodynamic environment of hypersonic vehicles is extremely complex. Their flow fields contain shock waves, vortices, viscous layers, shear layers, and their interactions and interferences. These not only increase drag and affect aerodynamic efficiency but also generate intense aerodynamic heating, impacting flight safety. For example, in October 1967, during a record-breaking flight at Mach 6.72, the X-15 experienced severe interference between shock waves generated by its engine nacelles and support shock waves, resulting in multiple areas of serious damage to the aircraft. The complex interactions and interferences in hypersonic flows pose significant challenges to the research and development of hypersonic vehicles.

[0004] Currently, research on the aerodynamic problems of hypersonic vehicles is mainly conducted through numerical simulation, wind tunnel experiments, and flight tests. However, wind tunnel experiments and flight tests are extremely costly, and experimental measurement techniques are limited, making them unsuitable as routine research methods for complex flows. Numerical simulation, due to its convenience, flexibility, and lower cost, has become an important tool for studying complex flow problems in recent years.

[0005] However, traditional numerical simulation methods suffer from low simulation accuracy and poor robustness in simulating hypersonic vehicles, leaving considerable room for improvement. On one hand, for flux schemes, traditional one-dimensional Riemann solvers exhibit insufficient computational accuracy and reduced stability when simulating complex multidimensional flows. Balsara proposed a true multidimensional Riemann solver based on a "corner frame" model. This solver achieves the multidimensional effects of the scheme simply and efficiently by solving the two-dimensional Riemann flux at the corners of the grid cells. However, this scheme only possesses first-order accuracy and cannot meet the needs of solving complex flows. On the other hand, scheme reconstruction is also an important method for improving computational accuracy. Traditional schemes such as MUSCL, TVD, and WENO have been widely applied in various flow solvers after years of development. However, multidimensional Riemann solvers have four discontinuous initial values, requiring the reconstruction of four state variables at the corners. Traditional high-order reconstruction methods suitable for structured meshes cannot achieve variable reconstruction at cell corners, thus limiting their application to higher-order accuracy. The WENO scheme, due to its high accuracy and essentially non-oscillating characteristics, is an ideal solution for reconstruction processes; however, apart from the ADER-WENO scheme developed by Balsara, it has been rarely studied. Therefore, researching multidimensional higher-order reconstruction schemes is beneficial for further enhancing the advantages of multidimensional Riemann solvers and improving the computational accuracy for complex flows. Summary of the Invention

[0006] To address the problems existing in the current technology and improve the simulation accuracy of shock waves and complex flow disturbances in hypersonic vehicles, this invention proposes a high-precision numerical simulation method for hypersonic vehicles, which includes a multi-dimensional Riemannian solver and a multi-dimensional fifth-order WENO scheme. This method applies a multi-dimensional high-order reconstruction scheme to a multi-dimensional Riemannian solver, providing technical support for more accurate hypersonic numerical simulation tasks and hypersonic vehicle design.

[0007] This invention first addresses the shortcomings of traditional CFD solvers, such as low computational accuracy and a small number of stable CFLs, by employing a multidimensional Riemannian solver to calculate interface flux. Secondly, to address the difficulty of applying traditional high-order reconstruction schemes to multidimensional Riemannian solvers, a fifth-order WENO scheme based on dimension splitting is developed to complete multidimensional reconstruction. Finally, to address issues such as insufficient computational efficiency and oscillations in high-order schemes, discontinuity probing techniques and accuracy-preserving limiters are employed to improve reconstruction efficiency and robustness.

[0008] The technical solution of this invention is as follows:

[0009] The high-precision numerical simulation method applicable to hypersonic vehicles includes the following steps:

[0010] Step 1: Establish the geometric model and computational domain of the hypersonic vehicle, divide and read the computational mesh of the hypersonic vehicle, obtain mesh information, and assign initial values ​​according to the incoming flow conditions;

[0011] Step 2: Discretize the governing equations to obtain a semi-discrete scheme in finite volume form;

[0012] Step 3: Obtain the reconstructed values ​​on both sides of the grid cell interface based on the mean value of the grid cells in the computational domain of the hypersonic vehicle;

[0013] Step 4: Using the reconstructed values ​​on both sides of all mesh element interfaces in Step 3, obtain the reconstructed values ​​at the Gaussian points of all mesh element interfaces in the hypersonic flow field for use in the multidimensional Riemann solver.

[0014] Step 5: Use the multidimensional Riemann solver obtained in Step 4 to obtain the interface flux;

[0015] Step 6: Determine the residual based on the interface flux obtained in Step 5, and perform time-progression solution to obtain the final hypersonic vehicle flow field.

[0016] Furthermore, step 2 includes the following:

[0017] The differential form governing equations for the flow around the aircraft are as follows:

[0018]

[0019] in,

[0020]

[0021] t represents time, x and y represent the abscissa and ordinate of the hypersonic vehicle in the computational domain, respectively, q is a conserved variable, f and g represent the flux in the x and y directions, and ρ, u, v, p, and E represent the hypersonic flow field density, velocity in the x direction, velocity in the y direction, pressure, and energy, respectively.

[0022] Spatial discretization of the flux term yields:

[0023]

[0024] Where i and j are the grid cell node numbers, and Δx and Δy represent the width of a single grid cell in the x and y directions, respectively; f i+1 / 2,j and g i,j+1 / 2 The interface fluxes in the x and y directions are respectively obtained by a higher-order reconstruction scheme and a multidimensional Riemann solver.

[0025] Furthermore, in step 3, for the mesh element V with node numbers i,j i,j Its unit mean is Then, by reconstructing along the x-direction using the unit mean, the mean values ​​on both sides of the interface i+1 / 2 are obtained. and The mean values ​​on both sides of the interface j+1 / 2 are obtained by reconstructing along the y-direction using the unit mean values. and

[0026] Furthermore, in step 3, for the mesh element V with node numbers i,j i,j ,calculate At that time, three sub-templates are created along the x-direction as follows:

[0027] S0={V i-2,j V i-1,j V i,j},S1={V i-1,j V i,j V i+1,j},S2={V i,j V i+1,j V i+2,j}

[0028] The global coordinates of the computational domain are converted into local coordinates within the cell, and the x-coordinate of the interpolation point is denoted as x. G =x i-1 / 2 +αΔx, where x i-1 / 2,j Representation unit V i,j The horizontal axis on the left side of the interface, Δx, represents the unit V. i,j Width in the x-direction; α in element V i,j The range of values ​​for α is [0,1]. After coordinate transformation, the interpolation polynomial of each sub-template is denoted as... for Taking α = 1, the interpolation polynomials are as follows:

[0029]

[0030]

[0031]

[0032] Calculate the smoothness factors IS0, IS1, IS2 for each polynomial:

[0033]

[0034]

[0035]

[0036] Calculate the ideal weights γ0, γ1, γ2, and the non-oscillating weights.

[0037]

[0038]

[0039]

[0040]

[0041]

[0042] τ5=|IS0-IS2|

[0043] in, All represent intermediate calculation variables with no particular physical meaning. ε represents a small constant used to avoid the denominator being zero. Finally, the reconstructed value is obtained.

[0044]

[0045] calculate At that time, establish V along the x-direction i+1,j Using the template centered at α, and setting α = 0; calculate At that time, establish V along the y direction i,j Using the template centered at α, and setting α = 1; calculate At that time, establish V along the y direction i,j+1 The template is centered, and α = 0.

[0046] Furthermore, step 4 includes the following steps:

[0047] Step 4.1: Obtain the preliminary reconstruction values ​​at the Gaussian points of the cell interface using the multidimensional fifth-order WENO method;

[0048] Step 4.2: Use a limiter to restrict the initial reconstruction value obtained in step 4.1, thereby obtaining the reconstruction value at the Gaussian point of the element interface, and use it for the multidimensional Riemann solver.

[0049] Furthermore, in step 4.1, for the grid interface i+1 / 2, the left side of the interface has a mean value of The unit is The right side of the interface has a mean value. The unit is For the fifth-order scheme, the interface has three Gaussian points G1, G2, and G3, with ordinates of G1, G2, and G3 respectively. y j2 =y j+1 / 2 ,

[0050] For the second Gaussian point G2, the flux is calculated using a two-dimensional Riemann solver, let in These represent the initial reconstructed values ​​of the lower left, upper left, lower right, and upper right corners of the Gaussian point G2, respectively.

[0051] For the first Gaussian point G1 and the third Gaussian point G3, calculate the preliminary reconstruction values ​​for the four regions respectively.

[0052] Reconstruct the variable at the lower left corner of Gaussian point G1 The reconstructed variable at the top left corner of Gaussian point G3 Create three sub-templates as follows:

[0053]

[0054]

[0055] The global coordinates of the computational domain are converted into local coordinates within the cell, and the coordinates of the interpolation point are denoted as y. G =y j-1 / 2 +αΔy, where y j-1 / 2 Representation unit The bottommost ordinate, Δy, represents the element. Width in the y-direction; α in the element The range of values ​​for α is [0,1]. After coordinate transformation, the interpolation polynomial of each sub-template is denoted as... for for The reconstruction polynomial built on the sub-template is as follows:

[0056]

[0057]

[0058]

[0059] Calculate the smoothing factors IS0, IS1, IS2 for each sub-template:

[0060]

[0061]

[0062]

[0063] Calculate the ideal weights γ0, γ1, γ2:

[0064]

[0065]

[0066]

[0067] Determine whether special handling is needed for negative weights based on the sign of the ideal weights. If min(γ0,γ1,γ2)>0, then directly calculate the non-oscillating weights.

[0068]

[0069]

[0070] τ5=|IS0-IS2|

[0071] in, All variables are intermediate calculation variables, with ε representing a small constant used to avoid zero denominators; the reconstructed values ​​are calculated based on non-oscillating weights and interpolation polynomials. and

[0072]

[0073] If min(γ0,γ1,γ2)<0, then the ideal weights are divided into positive and negative parts.

[0074]

[0075] θ represents the weighting parameter, from which the new positive and negative ideal weights are calculated.

[0076]

[0077] Where, σ ± The principle is to sum the positive and negative weights, and then calculate the non-oscillating weights of the positive and negative parts.

[0078]

[0079] in, Using intermediate calculation parameters, a reconstructed polynomial Q with positive and negative weights is formed. ± (α):

[0080]

[0081] Finally, the reassembled values ​​at G1 and G3 are calculated using interpolation polynomials:

[0082]

[0083] For the reconstructed variable at the top left corner of Gaussian point G1 The reconstructed variable at the lower left corner of Gaussian point G3 Establish separately and Reconstruct the template centered on α, where the values ​​of α are respectively

[0084] The required reconstruction values ​​on the right side of the interface are reconstructed using the average value on the right side of the cell.

[0085] Furthermore, in step 4.2, for This value is located in grid cell V. i,j Inside, firstly, with grid cell V i,j Define a 3×3 grid cell range centered on the shock detector and calculate the shock detector indicator factor (Detector):

[0086]

[0087]

[0088] in, σ represents the mean of the flow field variables within the 9 cells, and σ represents the variance of the flow field variables within the 9 cells. The density in the flow field is used to calculate the shock detector indicator factor. The size of the detector determines whether a limiter should be applied.

[0089]

[0090]

[0091]

[0092]

[0093] in Indicates the limiter indicator factor, when When indicates that the limiter is not used, The time indicates that the limiter is activated, M V and m V The maximum and minimum values ​​of the reconstructed value are represented by M and m, respectively, where M and m represent the maximum and minimum values ​​of the original variables within the 9 cells. The final reconstructed value after constraints is:

[0094]

[0095] Furthermore, in step 5, the flux at all Gaussian points is calculated and combined using one-dimensional and two-dimensional Riemann solvers to obtain the interface flux.

[0096] Furthermore, in step 5, for interface i+1 / 2, first calculate the flux at the first Gaussian point G1 of interface i+1 / 2; make the following variable substitutions:

[0097]

[0098] Where, q LD ,q LU ,q RD ,q RU To replace the reconstructed variables at the lower left, upper left, lower right, and upper right corners of the Gaussian point, calculate the wave velocity based on the reconstructed variables:

[0099]

[0100] Among them, S R and S L S represents the maximum and minimum wave velocities along the x-direction in the multidimensional wave propagation model. U and S D This represents the maximum and minimum wave velocities propagating along the y-direction in a multidimensional wave propagation model. Represents the state variable q RU Maximum wave velocity in the x-direction, Where u RU c represents the velocity in the x-direction. RU Indicates the speed of sound; Represents state q RU The minimum wave velocity in the x-direction at point X. It means (q) LU ,q RU The maximum wave velocity along the x-direction in the average state of Roe between ) This represents the average velocity of Roe in the x-direction. Let Roe represent the average speed of sound. It means (q) LU ,q RU The minimum wave velocity along the x-direction of the Roe average state between )

[0101] In the multidimensional wave propagation model, the two-dimensional Riemann flux is expressed as:

[0102]

[0103] Among them, f LD ,f LU ,f RD ,f RU The multidimensional wave propagation model represents the flux in the x-direction at the lower left, upper left, lower right, and upper right corners, respectively, where g is the flux. LD ,g LU ,g RD ,g RU These represent the flux in the y-direction at the lower left, upper left, lower right, and upper right corners of the multidimensional wave propagation model, respectively.

[0104] The flux at the third Gaussian point G3 It was obtained using the same method as Gaussian point G1;

[0105] For the second Gaussian point G2, the flux is calculated using a one-dimensional Riemann solver:

[0106]

[0107] Among them, f L ,f R These represent the flux on the left and right sides, respectively. These represent the reconstructed values ​​on the left and right sides of the Gaussian point after constraints, respectively. The superscript "m" represents the physical quantity related to the midpoint of the interface. and The maximum wave speeds representing leftward and rightward propagation in the respective wave propagation models are defined as follows:

[0108]

[0109] In the formula, 'a' represents the local speed of sound, and the superscript '~' indicates the Roe average.

[0110] Finally, the interface flux f is obtained by applying Gaussian integral. i+1 / 2,j :

[0111]

[0112] Furthermore, step 6 includes the following steps:

[0113] The semi-discrete finite volume scheme is transformed into a fully discrete finite volume scheme in spacetime using the third-order Runge-Kutta discretization formula for the time variable:

[0114] q (1) =q n +ΔtL(q n )

[0115]

[0116]

[0117] L represents the residual, and the superscript "n" represents the time step. The flow field variable values ​​at the next time step are solved using a spatiotemporally discrete finite volume scheme. This process is repeated to obtain the numerical simulation results of the full flow field stability of the hypersonic vehicle.

[0118] Beneficial effects

[0119] Compared with the prior art, the present invention has the following advantages:

[0120] 1. This invention utilizes a two-dimensional Riemann solver to calculate interface flux. Compared to a one-dimensional Riemann solver, it provides higher resolution for shock discontinuities and contact discontinuities in the flow field of hypersonic vehicles with the same mesh size.

[0121] 2. This invention utilizes a multi-dimensional reconstruction strategy to apply the fifth-order WENO reconstruction method to a two-dimensional Riemann solver. Compared with current low-order reconstruction methods, it has higher accuracy in capturing complex flow structures such as shock waves and flow disturbances, and can still maintain high accuracy when the mesh size is small, enabling more precise simulation of the flow field of hypersonic vehicles.

[0122] Additional aspects and advantages of the invention will be set forth in part in the description which follows, and in part will be obvious from the description, or may be learned by practice of the invention. Attached Figure Description

[0123] The above and / or additional aspects and advantages of the present invention will become apparent and readily understood from the description of the embodiments taken in conjunction with the following drawings, in which:

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

[0125] Figure 2 This is a schematic diagram of the reconstructed left and right average values ​​of the grid interface.

[0126] Figure 3 The shock wave detection factor is calculated within the selected grid range.

[0127] Figure 4 This is a schematic diagram of the multidimensional reconstruction process.

[0128] Figure 5 This is a schematic diagram of a multidimensional wave propagation model using a multidimensional Riemann solver.

[0129] Figure 6 This is a schematic diagram of a one-dimensional Riemann solver wave propagation model.

[0130] Figure 7 These are two-dimensional Riemann problem density contour maps calculated based on a low-order scheme and the present invention in Example 1; (a) second-order reconstruction scheme, (b) fifth-order reconstruction scheme.

[0131] Figure 8 The following are the radial Riemann problem density contour maps obtained by using a one-dimensional Riemann solver and the scheme of the present invention in Example 2: (a) one-dimensional Riemann solver and second-order reconstruction scheme, (b) two-dimensional Riemann solver and fifth-order reconstruction scheme.

[0132] Figure 9These are the supersonic flow pressure cloud diagrams of the NACA0012 airfoil obtained using the low-order reconstruction scheme and the scheme of the present invention in Example 3; (a) second-order reconstruction scheme, (b) fifth-order reconstruction scheme.

[0133] Figure 10 The images are Mach number cloud diagrams of hypersonic double elliptical flow obtained by using the low-order reconstruction scheme and the scheme of the present invention in Example 4; (a) second-order reconstruction scheme, (b) fifth-order reconstruction scheme.

[0134] Figure 11 This is the hypersonic double elliptical flow shock wave detector cloud map obtained by using the low-order reconstruction scheme and the scheme of the present invention in Example 4. Detailed Implementation

[0135] The present invention will be further described in detail below using hypersonic vehicle flow simulation as an example in conjunction with the accompanying drawings. The embodiments described are exemplary and intended to explain the present invention, and should not be construed as limiting the present invention.

[0136] Step 1: Establish the geometric model and computational domain of the hypersonic vehicle, divide and read the computational grid of the hypersonic vehicle, obtain grid information, and assign initial values ​​according to the incoming flow conditions;

[0137] A hypersonic computational mesh is generated based on the hypersonic vehicle flow problem to be analyzed. Mesh information is read to obtain the vehicle mesh scale, wall node coordinates, etc. Then, initial values ​​are assigned to the computational domain based on the incoming flow conditions, such as the mesh element V. i,j The unit mean is

[0138] Step 2: Construct a semi-discrete scheme in finite volume form

[0139] The governing equations in differential form are as follows:

[0140]

[0141] in,

[0142]

[0143] Where t represents time, x and y represent the horizontal and vertical axes respectively, q is a conserved variable, f and g represent the flux in the x and y directions respectively, and ρ, u, v, p, and E represent the density, velocity in the x direction, velocity in the y direction, pressure, and energy of the hypersonic vehicle flow field respectively.

[0144] Spatial discretization of the flux term yields:

[0145]

[0146] Where i and j are the cell node numbers; f i+1 / 2,j and g i,j+1 / 2 The interface numerical fluxes in the x and y directions are respectively obtained by a higher-order reconstruction scheme and a multidimensional Riemann solver. The specific solution process is detailed in the following steps.

[0147] Step 3: Obtain the reconstructed values ​​on both sides of the grid cell interface based on the mean value of the grid cells in the computational domain of the hypersonic vehicle.

[0148] The mean values ​​of the interfaces on both sides of the unit and the multidimensional reconstruction are both based on the general fifth-order WENO method, which is as follows.

[0149] Assume a uniform grid cell of one dimension is I i =(x i-1 / 2 ,x i+1 / 2 The grid spacing is Δx, and there exists an unknown function distribution f(x) on the grid cells. The integral average of f(x) over these cells is known to be... Using unit I i The five adjacent units centered on {I i-2 ,I i-1 ,I i ,I i+1 ,I i+2} is used as an interpolation template to reconstruct the fifth-order approximation of f(x) at any interpolation point in the unit V, and this reconstructed polynomial is denoted as Q(x).

[0150] In the general WENO method, the fourth-order reconstruction polynomial Q(x) with fifth-order accuracy should be represented as a combination of three quadratic polynomials to achieve high accuracy in smooth regions and suppress oscillations near discontinuities. Therefore, the above five-element interpolation template first needs to be divided into the following three sub-templates:

[0151] S0={I i-2 ,I i-1 ,I i},S1={I i-1 ,I i ,I i+1},S2={I i ,I i+1 ,I i+2}

[0152] Then, construct a quadratic polynomial p with third-order precision on each sub-template. j (x), j = 0, 1, 2, each polynomial has a similar expression p j (x)=a j x 2 +b j x+c jThe three interpolation polynomials must satisfy the following conditions:

[0153]

[0154]

[0155]

[0156] in Representation Unit I j The mean of the polynomials can be obtained by substituting their respective expressions into the above equation. For ease of calculation, the global coordinates of the computational domain are converted to local coordinates within the cell. Let the coordinates of the interpolation points be x. G =x i-1 / 2 +αΔx, where x i-1 / 2 Representation Unit I i The leftmost coordinate, α, is in cell I. i The range of values ​​for α is [0,1], and the interpolation polynomial after coordinate transformation is denoted as . The point x is obtained through calculation and transformation. G The fifth-order approximation at is:

[0157]

[0158]

[0159]

[0160] Similarly, a fourth-order interpolation polynomial Q(x) with fifth-order accuracy can be constructed on a smooth, continuous five-element template, which also needs to satisfy the following conditions:

[0161]

[0162] Similarly, by converting the coordinates to local coordinates, the polynomial can be obtained, denoted as .

[0163] Ideally, three quadratic polynomials The weighted combination should be equal to Therefore, ideal weights γ0, γ1, γ2 are introduced such that the interpolation polynomial satisfies the following within any grid cell:

[0164]

[0165] By comparing the polynomial coefficients on both sides of the equation, we can obtain the three ideal weighting coefficients:

[0166]

[0167]

[0168]

[0169] In non-ideal cases (near discontinuities), the weights of sub-templates containing discontinuities should be very small to avoid non-physical oscillations. Therefore, a method is needed to determine the smoothness of the interpolation sub-templates and calculate new weights for each template based on the ideal weights and the smoothness. These new weights can prevent non-physical oscillations in the weighted interpolation polynomial, and are therefore called oscillation-free weights. The parameter representing the smoothness of the sub-template is called the smoothness factor. The smoothness factors for the three sub-templates are denoted by IS0, IS1, and IS2, respectively, and their calculation method is as follows:

[0170]

[0171]

[0172]

[0173] Then, based on the smoothing factor and the ideal weights, the non-oscillating weights can be obtained:

[0174]

[0175]

[0176] τ5=|IS0-IS2|

[0177] These are all intermediate calculation variables and do not need to be given actual physical meaning. ε represents a small constant to avoid the denominator being zero, and can be taken as ε = 1 × 10. -16 Therefore, the interpolation polynomial in any case can be expressed as:

[0178]

[0179] The following is an explanation The specific implementation process is as follows: the reconstructed variables in other positions are calculated using the same method. For example... Figure 2 As shown, let the mean be The unit is V i,j The mean values ​​of the left and right sides of the interface i+1 / 2 are obtained by reconstructing along the x-direction using the unit mean values. and Then, the mean values ​​on both sides of the interface j+1 / 2 are obtained by reconstructing along the y-direction using the unit mean values. and beg At that time, three sub-templates are created along the x-direction as follows:

[0180] S0={V i-2,j V i-1,jV i,j},S1={V i-1,j V i,j V i+1,j},S2={V i,j V i+1,j V i+2,j}

[0181] The global coordinates of the computational domain are converted to local coordinates within the cell, and the coordinates of the interpolation point are denoted as x. G =x i-1 / 2 +αΔx, where x i-1 / 2,j Representation unit V i,j The horizontal axis on the left side of the interface, Δx, represents the unit V. i,j Width in the x-direction. α in element V i,j The range of values ​​for α is [0,1]. After coordinate transformation, the interpolation polynomial of each sub-template is denoted as... for Taking α = 1, the interpolation polynomials established by the three sub-templates are as follows:

[0182]

[0183]

[0184]

[0185] Calculate the smoothness factors IS0, IS1, IS2 for each polynomial:

[0186]

[0187]

[0188]

[0189] Calculate the ideal weights and the non-oscillating weights:

[0190]

[0191]

[0192]

[0193]

[0194]

[0195] τ5=|IS0-IS2|

[0196] Finally, calculate the reconstructed value.

[0197]

[0198] for Establish along the x-direction with V i+1,j Using a template centered at α = 0, we can obtain the following using the same method. for and Establish V along the y direction respectively i,j and V i,j+1 Using a template centered at the center, we can solve the problem by setting α = 1 and 0 respectively. This allows us to obtain the mean values ​​on both sides of all mesh interfaces.

[0199] Step 4: Calculate the reconstruction values ​​required by the multidimensional Riemann solver.

[0200] The multidimensional reconstruction is still based on the general fifth-order WENO method. However, since the reconstruction location is a Gaussian point, its ideal weights may be negative. In the one-dimensional case, α at the interpolation point is generally 0 or 1, and the ideal weights γ0, γ1, and γ2 are all positive. But when α takes other values, the ideal weights may be negative. This invention requires multidimensional reconstruction, where α at the interpolation point is not 0 or 1, thus requiring special handling. Furthermore, although based on the WENO method, certain constraints are still needed to prevent numerical oscillations.

[0201] When the ideal weights are negative, a simple separation method can be used to improve the interpolation polynomial. When min(γ0,γ1,γ2)<0, the ideal weights can be decomposed into... and Two parts: positive and negative

[0202]

[0203] θ is the weighting parameter, typically set to 3. Then, the ideal weights for the positive and negative parts are constructed.

[0204]

[0205] Where, σ ± The sum of the positive and negative parts from the previous decomposition yields two separate polynomials Q. ± (α):

[0206]

[0207] These two polynomials satisfy the following conditions:

[0208] Q(α)=σ + Q + (α)-σ - Q - (α)

[0209] Similarly, we need to consider the smoothness factor and the non-oscillating weights, and calculate them using the ideal weights of the positive and negative parts respectively:

[0210]

[0211] That is, the non-oscillating weights of the positive and negative parts. This represents an intermediate parameter, where ε represents a small constant, which can be taken as ε = 1 × 10⁻⁶. -16 This avoids the denominator being zero. Then, the polynomial Q with positive and negative splitting can be obtained. ± (α):

[0212]

[0213] The final WENO reconstruction formula is:

[0214] Q(α)=σ + Q + (α)-σ - Q - (α)

[0215] According to the fifth-order WENO method, when reconstructing the variables needed on the i+1 / 2 interface, we can first apply the fifth-order WENO formula along the x-direction to perform a one-dimensional reconstruction, obtaining the mean values ​​of the left and right sides of the interface. and Then, using the average values ​​of both sides of the interface obtained from the previous reconstruction, a one-dimensional reconstruction is performed along the y-direction to obtain the four inputs required for the multidimensional Riemann solver.

[0216] Theoretically, the above strategy can complete the reconstruction and meet the requirement of being essentially oscillatory-free. However, in complex flows, numerical oscillations still exist near discontinuities. Therefore, this invention introduces a limiter to impose a secondary constraint on the reconstructed value to avoid the emergence of new extrema.

[0217] Assuming the flow field variable q near the discontinuity satisfies q∈[m,M], then the reconstructed variable at the Gaussian point should also satisfy q. Gauss ∈[m,M]. When calculating the flux in the x-direction, q Gauss include The reconstruction obtained directly using the above-mentioned fifth-order WENO method generally cannot satisfy q. Gauss Since the integers are ∈[m,M], certain methods are needed to restrict them.

[0218] Assume unit V i,j The interpolation polynomial with fifth-order precision is p i,j (x,y), the restricted polynomial is expressed as The following methods can be used to impose restrictions:

[0219]

[0220] in,

[0221]

[0222]

[0223] In the above formula, M and m represent the extreme values ​​of the flow field solution, defined as the maximum and minimum values ​​of the target element and its neighboring elements. V and m V Representing the extrema of the reconstructed polynomial, the polynomial distribution within a cell is generally difficult to obtain; therefore, the reconstructed value is used for calculation.

[0224]

[0225]

[0226] The limiter only needs to be activated near the discontinuity, not in the smooth region. Therefore, a shock detector is needed to determine the location of the discontinuity in the computational domain. The discontinuity detector is defined as follows.

[0227]

[0228]

[0229] The above formula σ represents the mean and standard deviation of the flow field variables within a certain range. Generally, a 3×3 range is chosen for calculating the intermittent detector, such as... Figure 3 As shown, n = 9, and h represents the grid scale. Therefore, the limiter combined with the shock detector can be expressed as:

[0230]

[0231]

[0232] By combining the general WENO reconstruction method with a shock detector and a limiter, and employing a dimensional splitting strategy, a complete multidimensional fifth-order reconstruction strategy can be formed. Calculate f. i+1 / 2,j The multidimensional reconstruction strategy and flux calculation method are as follows:

[0233] Integral average using a two-dimensional grid A one-dimensional reconstruction is performed in the x-direction to obtain a higher-order approximation of the variable q at interface i+1 / 2 in the x-direction. Two values ​​can be reconstructed based on the left- and right-leaning templates, denoted as... and Let the mean of the grid integral be The mesh element is V i,j Based on a left-leaning template, i.e., V i,j A template centered on the right, or a template based on the V-shape. i+1,j The template is centered on this. For a one-dimensional Riemann solver, this is the input state quantity required for its computation.

[0234] The average value of the grid integral obtained by one-dimensional reconstruction in the x-direction respectively and Find the Gaussian point y along the y-direction jk A higher-order approximation on the surface. Here, it is assumed that the average value on the left side of the interface is... The unit is The average value on the right side of the interface is The unit is For the fifth-order scheme, there are 3 Gaussian points, namely... y j2 =y j+1 / 2 ,

[0235] The first and third points are two-dimensional Riemann problems, while the second point is actually a one-dimensional Riemann problem. It can be solved using the formula of a one-dimensional Riemann solver or by directly applying a two-dimensional Riemann solver.

[0236] For the first and third points and Use to The template reconstruction is centered on α, with the value of α being... and Using respectively and Reconstruct the template centered on α, where α takes the value of

[0237] For the second point, if a one-dimensional Riemann solver is used, only... and The two reconstructed values, if solved using a two-dimensional Riemannian solver, can be made... Similarly, the reconstruction value required on the right side of interface i+1 / 2 is reconstructed using the unit on the right side of the interface, and the value of α is the same.

[0238] After the reconstruction described above, a preliminary reconstructed value has been obtained at each Gaussian point, denoted as . Then a limiter is applied to each value to prevent numerical oscillations. For example, for The value after restriction The value of is determined by referring to the limiter calculation method described above.

[0239] Finally, the interface flux is calculated using the reconstructed values ​​after constraints. For any Gaussian point, the aforementioned two-dimensional Riemann solver is used. Solve for the flux at this point.

[0240] For interface i+1 / 2, the detailed implementation process is as follows; the method for reconstructing the grid interface in other positions is the same.

[0241] like Figure 4 As shown, assume that the left side of the interface i+1 / 2 has a mean value of The unit is The right side of the interface has a mean value. The unit is For the fifth-order scheme, the interface has three Gaussian points G1, G2, and G3, with ordinates of G1, G2, and G3 respectively. The second Gaussian point G2 can be considered a one-dimensional Riemann problem or a special two-dimensional Riemann problem. If we want to use a two-dimensional Riemann solver to calculate the flux, then let... These represent the initial reconstructed values ​​for the lower left, upper left, lower right, and upper right corners of the Gaussian point G2, respectively. Then, G1 and G3 are two-dimensional Riemannian problems, and the reconstruction variables for the four regions need to be calculated separately.

[0242] Reconstruct the variable at the lower left corner of Gaussian point G1 The reconstructed variable at the top left corner of Gaussian point G3 Create three sub-templates as follows:

[0243]

[0244]

[0245] The global coordinates of the computational domain are converted into local coordinates within the cell, and the coordinates of the interpolation point are denoted as y. G =y j-1 / 2 +αΔy, where y j-1 / 2 Representation unit The bottommost ordinate, Δy, represents the element. Width in the y-direction. α in the element. The range of values ​​for α is ∈ [0,1]. After coordinate transformation, the interpolation polynomial of each sub-template is denoted as... for for The reconstruction polynomial built on the sub-template is as follows:

[0246]

[0247]

[0248]

[0249] Calculate the smoothing factors IS0, IS1, IS2 for each sub-template:

[0250]

[0251]

[0252]

[0253] Calculate the ideal weights γ0, γ1, γ2:

[0254]

[0255]

[0256]

[0257] Determine whether special handling is needed for negative weights based on the sign of the ideal weights. If min(γ0,γ1,γ2)>0, then directly calculate the non-oscillating weights.

[0258]

[0259]

[0260] τ5=|IS0-IS2|

[0261] in, All are intermediate calculation variables, where ε represents a small constant, which can be taken as ε = 1 × 10 -16 This avoids zero denominators and directly calculates reconstructed variables based on non-oscillating weights and sub-template interpolation polynomials.

[0262]

[0263]

[0264] If min(γ0,γ1,γ2) < 0, then the ideal weights need to be divided into positive and negative parts.

[0265]

[0266] θ represents the weighting parameter, typically taken as 3, based on which the new positive and negative ideal weights are calculated.

[0267]

[0268] Where, σ ± The principle is to sum the positive and negative weights, and then calculate the non-oscillating weights of the positive and negative parts. The smoothing factor is the same as before:

[0269]

[0270] in, ε is an intermediate calculation parameter, representing a small constant, which can be taken as ε = 1 × 10. -16 This avoids the denominator being zero and utilizes two weights to form a reconstructed polynomial Q with positive and negative separation. ± (α):

[0271]

[0272] Finally, the reassembled values ​​at G1 and G3 are calculated using interpolation polynomials:

[0273]

[0274]

[0275] In addition, for the reconstructed variable at the top left corner of G1 The refactored variable at the bottom left of G3 It is necessary to establish them separately. and Reconstruct the template centered on α, where α takes the value of The required reconstruction values ​​on the right side of the interface are reconstructed using the mean value on the right side of the cell. Similarly, the reconstruction values ​​at all Gaussian points on cell interface j+1 / 2 are reconstructed along the x-direction using the mean values ​​on both sides of interface j+1 / 2. In this way, preliminary reconstruction values ​​at all Gaussian points on the interface can be obtained.

[0276] Use limiters to restrict the initial refactored variables:

[0277] for For example, this value is located in cell V. i,j Inside, firstly with V i,j Define a 3×3 grid cell range centered on the shock detector and calculate the indicator factor Detector:

[0278]

[0279]

[0280] in, σ represents the mean of the flow field variables within the 9 cells, and σ represents the variance of the variables within the 9 cells. The size of the detector determines whether a limiter should be applied.

[0281]

[0282]

[0283]

[0284]

[0285] in Indicates the limiter indicator factor, when When indicates that the limiter is not used, The time indicates the application limiter, M V and m V This represents the maximum and minimum values ​​of the reconstructed value, where M and m represent the maximum and minimum values ​​of the original variable within the 9 cells. The final reconstructed value after constraints is:

[0286]

[0287] Based on this method, the initial reconstructed values ​​at all Gaussian points are judged and restricted.

[0288] Step 5: Calculate the interface flux

[0289] To fully consider the information propagating along the interface normal and laterally and improve the calculation accuracy when solving for interface flux, a two-dimensional Riemann solver is used to calculate the interface flux.

[0290] A schematic diagram of the multidimensional wave propagation model is shown below. Figure 5 As shown, according to the multidimensional wave propagation model, the grid corners contain q. LD ,q LU ,q RD ,q RU The four discontinuous initial values ​​constitute a two-dimensional Riemann problem, from which the flux in the x-direction can be obtained as:

[0291]

[0292] In the above formula, f LD ,f LU ,f RD ,f RU The multidimensional wave propagation model represents the flux in the x-direction at the lower left, upper left, lower right, and upper right corners, respectively, where g is the flux. LD ,g LU ,g RD ,g RU S represents the flux in the y-direction at the lower left, upper left, lower right, and upper right corners of the multidimensional wave propagation model, respectively. R and S L S represents the maximum and minimum wave velocities along the x-direction in the multidimensional wave propagation model. U and S D These represent the maximum and minimum wave velocities propagating along the y-direction in the multidimensional wave propagation model. The expressions for these four limiting wave velocities are as follows:

[0293]

[0294] in, Represents state q RU Maximum wave velocity in the x-direction at a given location; Represents state q RU Minimum wave velocity in the x-direction; It means (q) LU ,q RU The maximum wave velocity along the x-direction in the average state of Roe between ) It means (q) LU ,q RU The minimum wave velocity of the Roe average state along the x-direction between )

[0295] Similarly, the flux in the y-direction can be expressed as:

[0296]

[0297] Conventional methods for solving interface flux consider the Riemann problem at the interface's midpoint and corners. The midpoint is treated as a one-dimensional Riemann problem, and the corners as a two-dimensional Riemann problem. Simpson's formula is used for weighting, and the flux f at interface i+1 / 2 is calculated. i+1 / 2,j The calculation formula is as follows:

[0298]

[0299] in, The corner flux is represented and calculated using a two-dimensional Riemann solver. The flux at the midpoint of the interface is represented and calculated using a one-dimensional Riemannian solver. The wave propagation model of the one-dimensional Riemannian solver is as follows: Figure 5 As shown, the calculation method is as follows:

[0300]

[0301] Among them, f L ,f R q represents the flux on the left and right sides. L ,q R This represents the flow field variables on the left and right sides. The superscript "m" represents the physical quantity related to the midpoint of the interface. and These represent the maximum wave speeds propagating to the left and right, respectively, and are defined as follows:

[0302]

[0303] In the above formula, 'a' represents the local speed of sound, and the superscript '~' indicates the Roe average.

[0304] However, the Simpson integral only has third-order accuracy, which cannot match the reconstruction scheme that achieves fifth-order accuracy. Therefore, this invention makes certain improvements to the interface flux calculation method, using a three-point Gaussian integral to achieve fifth-order accuracy, where f i+1 / 2,j and g i,j+1 / 2 Calculate using the following two formulas respectively:

[0305]

[0306] For the fifth-order scheme of this invention, we have β1 = 4 / 9, β2 = 5 / 18, β3 = 4 / 9. For interface i+1 / 2, the three Gaussian points are respectively For interface j+1 / 2, the three Gaussian points are respectively The second point is the midpoint of the interface. Since this point is equidistant from the interpolation templates on both sides, it can be considered a one-dimensional Riemann problem and calculated using a one-dimensional Riemann solver. However, for the sake of simplicity and consistency, the above formula is still expressed in two dimensions because when q LD =q LU ,q RD =q RU At that point, the two-dimensional formula degenerates into a one-dimensional Riemann solver.

[0307] The following describes in detail the calculation process of the flux at interface i+1 / 2; the calculation method is the same for other positions.

[0308] First, calculate the flux at the first Gaussian point G1 on the i+1 / 2 interface. After the multidimensional reconstruction in step 4, the reconstructed variables at the Gaussian point have been obtained. To simplify the formula, the following variable substitutions are made:

[0309]

[0310] First, calculate the wave speed:

[0311]

[0312] in, Represents state q RU Maximum wave velocity in the x-direction, Where u RU c represents the velocity in the x-direction. RU Indicates the speed of sound; Represents state q RU The minimum wave velocity in the x-direction at point X. It means (q) LU ,q RU The maximum wave velocity of the Roe average state along the x-direction between ) This represents the average velocity of Roe in the x-direction. Let Roe represent the average speed of sound. It means (q) LU ,q RU The minimum wave velocity along the x-direction of the Roe average state between )

[0313] The two-dimensional Riemann flux is then:

[0314]

[0315] The third Gaussian point is represented in the same way as above.

[0316] For the second Gaussian point, the flux is calculated using a one-dimensional Riemann solver.

[0317]

[0318] The superscript "m" represents the physical quantity related to the midpoint of the interface. and These represent the maximum wave speeds propagating to the left and right, respectively, and their values ​​are defined as follows:

[0319]

[0320] In the above formula, 'a' represents the local speed of sound, and the superscript '~' indicates the Roe average.

[0321] Finally, the interface flux is obtained by applying the Gaussian integral:

[0322]

[0323] Then, the flux through the interface j+1 / 2 along the y-direction can be calculated using the same method, but the calculation method for the multidimensional Riemann flux along the y-direction is slightly different from that in the x-direction, also based on... Figure 5 The multidimensional wave propagation model, flux at the first Gaussian point on interface j+1 / 2 The calculation method is as follows; the remaining points in the y-direction are calculated using the same method:

[0324]

[0325] The meanings of the symbols in the formula are the same as those in the text above, but q LD ,q LU ,q RD ,q RU We need to utilize the reconstructed value at the first Gaussian point on the interface y+1 / 2.

[0326] Step 6: Solve by time progression

[0327] The semi-discrete finite volume scheme is transformed into a fully discrete finite volume scheme in spacetime using the third-order Runge-Kutta discretization formula for the time variable:

[0328] q (1) =q n +ΔtL(q n )

[0329]

[0330]

[0331] L represents the residual, and the superscript "n" represents the time step. The flow field variables at the next time step are solved using a fully discrete finite volume scheme in spacetime. This process is repeated sequentially to obtain the numerical simulation results of the overall flow field stability of the hypersonic vehicle.

[0332] The following four implementation examples serve as specific implementation cases of the method disclosed in this invention.

[0333] Example 1: The Two-Dimensional Riemann Problem.

[0334] This problem describes the interaction between two shock waves and two contact discontinuities, with the following initial conditions:

[0335]

[0336] The computational grid is 1000×1000. Although this example is only a numerical example, it contains a relatively complex flow structure and can be used to demonstrate the high numerical accuracy of the present invention. Figure 7 The density distributions obtained at t=0.25 using the second-order reconstruction scheme and the fifth-order reconstruction scheme described in this invention are presented, with a total of 30 contour lines ranging from 0.54 to 1.70. It can be seen from the figure that the scheme described in this invention has higher shock wave resolution and captures more refined wave systems in the interference region compared to lower-order schemes.

[0337] Example 2: Radial Riemann Problem

[0338] The computational domain for this example is [0,1]×[0,1], the CFL number is 0.6, the computational grid is 200×200, and the initial conditions are as follows.

[0339]

[0340] Figure 8 Density contour plots obtained at t=0.13 using both the one-dimensional Riemann solver scheme and the scheme of this invention are presented. The figures show that, with the same mesh size, the two-dimensional Riemann solver of this invention, combined with a higher-order reconstruction scheme, achieves higher accuracy in capturing shock discontinuities and contact discontinuities.

[0341] Example 3: Supersonic Flow Around NACA0012

[0342] This example demonstrates supersonic flow around an airfoil, a crucial calculation required when designing supersonic and hypersonic aircraft. The incoming Mach number is 2.0, the angle of attack is 10°, and the mesh size is 300×300. Figure 9 The pressure cloud maps obtained by using the second-order reconstruction scheme and the scheme of the present invention are shown in the figure. It can be seen from the figure that when the number of grids is small, the high-order scheme of the present invention has higher accuracy in capturing the head and tail shock waves.

[0343] Example 4: Hypersonic Flow Around a Double Ellipse

[0344] This example illustrates hypersonic double-elliptical flow, which can represent the flow of hypersonic vehicles to some extent. The incoming Mach number is 8.15, the angle of attack is 30°, and the incoming flow density is 0.0231 kg / m³. 3 The incoming static pressure is 370.7 Pa and the incoming static temperature is 56 K. Figure 10 The figures show Mach number cloud maps obtained using the second-order reconstruction scheme and the fifth-order reconstruction scheme of this invention, respectively. It can be seen from the figures that the two schemes have comparable accuracy in simulating the detached shock wave at the aircraft nose, but the fifth-order scheme has higher accuracy in simulating the detached shock wave and secondary shock wave in front of the cockpit. Figure 11 The shock wave detector also shows that the fifth-order scheme provides a more precise solution for separating the shock wave, the secondary shock wave, and the nearby flow disturbances.

[0345] Although embodiments of the present invention have been shown and described above, it is understood that the above embodiments are exemplary and should not be construed as limiting the present invention. Those skilled in the art can make changes, modifications, substitutions and variations to the above embodiments within the scope of the present invention without departing from the principles and spirit of the present invention.

Claims

1. A high-precision numerical simulation method suitable for hypersonic vehicles, characterized in that: Includes the following steps: Step 1: Establish the geometric model and computational domain of the hypersonic vehicle, divide and read the computational mesh of the hypersonic vehicle, obtain mesh information, and assign initial values ​​according to the incoming flow conditions; Step 2: Discretize the governing equations to obtain a semi-discrete scheme in finite volume form; Step 3: Obtain the reconstructed values ​​on both sides of the grid cell interface based on the mean value of the grid cells in the computational domain of the hypersonic vehicle; Step 4: Using the reconstructed values ​​on both sides of all mesh element interfaces in Step 3, obtain the reconstructed values ​​at the Gaussian points of all mesh element interfaces in the hypersonic flow field for use in the multidimensional Riemann solver. Step 5: Use the multidimensional Riemann solver obtained in Step 4 to obtain the interface flux; Step 6: Determine the residual based on the interface flux obtained in Step 5, and perform time-progression solution to obtain the final hypersonic vehicle flow field.

2. The high-precision numerical simulation method for hypersonic vehicles according to claim 1, characterized in that: Step 2 includes the following: The differential form governing equations for the flow around the aircraft are as follows: in, t represents time, x and y represent the x and y coordinates in the computational domain of the hypersonic vehicle, respectively, q is a conserved variable, f and g represent the flux in the x and y directions, and ρ, u, v, p, and E represent the hypersonic flow field density, velocity in the x direction, velocity in the y direction, pressure, and energy, respectively. Spatial discretization of the flux term yields: Where i and j are the grid cell node numbers, and Δx and Δy represent the width of a single grid cell in the x and y directions, respectively; f i+1 / 2,j and g i,j+1 / 2 The interface fluxes in the x and y directions are respectively obtained by a higher-order reconstruction scheme and a multidimensional Riemann solver.

3. The high-precision numerical simulation method for hypersonic vehicles according to claim 1, characterized in that: In step 3, for the mesh element V with node numbers i and j i,j Its unit mean is Then, by reconstructing along the x-direction using the unit mean, the mean values ​​on both sides of the interface i+1 / 2 are obtained. and The mean values ​​on both sides of the interface j+1 / 2 are obtained by reconstructing along the y-direction using the unit mean values. and 4. The high-precision numerical simulation method for hypersonic vehicles according to claim 3, characterized in that: In step 3, for the mesh element V with node numbers i and j i,j ,calculate At that time, three sub-templates are created along the x-direction as follows: S0={V i-2,j ,V i-1,j ,V i,j },S1={V i-1,j ,V i,j ,V i+1,j },S2={V i,j ,V i+1,j ,V i+2,j } The global coordinates of the computational domain are converted into local coordinates within the cell, and the x-coordinate of the interpolation point is denoted as x. G =x i-1 / 2 +αΔx, where x i-1 / 2,j Representation unit V i,j The horizontal axis on the left side of the interface, Δx, represents the unit V. i,j Width in the x-direction; α in element V i,j The range of values ​​for α is [0,1]. After coordinate transformation, the interpolation polynomial of each sub-template is denoted as... for Taking α = 1, the interpolation polynomials are as follows: Calculate the smoothness factors IS0, IS1, IS2 for each polynomial: Calculate the ideal weights γ0, γ1, γ2, and the non-oscillating weights. τ5=|IS0-IS2| in, τ5 represents intermediate calculation variables with no particular physical meaning, ε represents a small constant used to avoid the denominator being zero, and finally the reconstructed value is obtained. calculate At that time, establish V along the x-direction i+1,j Using the template centered at α, and setting α = 0; calculate At that time, establish V along the y direction i,j Using the template centered at α, and setting α = 1; calculate At that time, establish V along the y direction i,j+1 The template centered on α is set to 0.

5. The high-precision numerical simulation method for hypersonic vehicles according to claim 1, characterized in that: Step 4 includes the following steps: Step 4.1: Obtain the preliminary reconstruction values ​​at the Gaussian points of the cell interface using the multidimensional fifth-order WENO method; Step 4.2: Use a limiter to restrict the initial reconstruction value obtained in step 4.1, thereby obtaining the reconstruction value at the Gaussian point of the element interface, and use it for the multidimensional Riemann solver.

6. The high-precision numerical simulation method for hypersonic vehicles according to claim 5, characterized in that: In step 4.1, for the grid interface i+1 / 2, the left side of the interface has a mean value of The unit is The right side of the interface has a mean value. The unit is For the fifth-order scheme, the interface has three Gaussian points G1, G2, and G3, with ordinates of G1, G2, and G3 respectively. For the second Gaussian point G2, the flux is calculated using a two-dimensional Riemann solver, let in These represent the initial reconstructed values ​​of the lower left, upper left, lower right, and upper right corners of the Gaussian point G2, respectively. For the first Gaussian point G1 and the third Gaussian point G3, calculate the preliminary reconstruction values ​​for the four regions respectively: Reconstruct the variable at the lower left corner of Gaussian point G1 The reconstructed variable at the top left corner of Gaussian point G3 Create three sub-templates as follows: The global coordinates of the computational domain are converted into local coordinates within the cell, and the coordinates of the interpolation point are denoted as y. G =y j-1 / 2 +αΔy, where y j-1 / 2 Representation unit The bottommost ordinate, Δy, represents the element. Width in the y-direction; α in the element The range of values ​​for α is [0,1]. After coordinate transformation, the interpolation polynomial of each sub-template is denoted as... for for The reconstruction polynomial built on the sub-template is as follows: Calculate the smoothing factors IS0, IS1, IS2 for each sub-template: Calculate the ideal weights γ0, γ1, γ2: Determine whether special handling is needed for negative weights based on the sign of the ideal weights. If min(γ0,γ1,γ2)>0, then directly calculate the non-oscillating weights. τ5=|IS0-IS2| in, τ5 are all intermediate calculation variables, and ε represents a small constant used to avoid zero denominators; the reconstructed values ​​are calculated based on non-oscillating weights and interpolation polynomials. and If min(γ0,γ1,γ2)<0, then the ideal weights are divided into positive and negative parts. θ represents the weighting parameter, from which the new positive and negative ideal weights are calculated. Where, σ ± The principle is to sum the positive and negative weights, and then calculate the non-oscillating weights of the positive and negative parts. in, Using intermediate calculation parameters, a reconstructed polynomial Q with positive and negative weights is formed. ± (α): Finally, the reassembled values ​​at G1 and G3 are calculated using interpolation polynomials: For the reconstructed variable at the top left corner of Gaussian point G1 The reconstructed variable at the lower left corner of Gaussian point G3 Establish separately and Reconstruct the template centered on α, where the values ​​of α are respectively The required reconstruction values ​​on the right side of the interface are reconstructed using the average value on the right side of the cell.

7. The high-precision numerical simulation method for hypersonic vehicles according to claim 5, characterized in that: In step 4.2, for This value is located in grid cell V. i,j Inside, firstly, with grid cell V i,j Define a 3×3 grid cell range centered on the shock detector and calculate the shock detector indicator factor (Detector): in, σ represents the mean of the flow field variables within the 9 cells, and σ represents the variance of the flow field variables within the 9 cells. The density in the flow field is used to calculate the shock detector indicator factor. The size of the detector determines whether a limiter should be applied. in Indicates the limiter indicator factor, when When indicates that the limiter is not used, when The time indicates that the limiter is activated, M V and m V The maximum and minimum values ​​of the reconstructed value are represented by M and m, respectively, where M and m represent the maximum and minimum values ​​of the original variables within the 9 cells. The final reconstructed value after constraints is:

8. The high-precision numerical simulation method for hypersonic vehicles according to claim 1, characterized in that: In step 5, the flux at all Gaussian points is calculated and combined using one-dimensional and two-dimensional Riemann solvers to obtain the interface flux.

9. The high-precision numerical simulation method for hypersonic vehicles according to claim 8, characterized in that: In step 5, for interface i+1 / 2, first calculate the flux at the first Gaussian point G1 of interface i+1 / 2; make the following variable substitutions: Where, q LD ,q LU ,q RD ,q RU To replace the reconstructed variables at the lower left, upper left, lower right, and upper right corners of the Gaussian point, calculate the wave velocity based on the reconstructed variables: Among them, S R and S L S represents the maximum and minimum wave velocities along the x-direction in the multidimensional wave propagation model. U and S D This represents the maximum and minimum wave velocities propagating along the y-direction in a multidimensional wave propagation model. Represents the state variable q RU Maximum wave velocity in the x-direction, Where u RU c represents the velocity in the x-direction. RU Indicates the speed of sound; Represents state q RU The minimum wave velocity in the x-direction at point X. It means (q) LU ,q RU The maximum wave velocity along the x-direction in the average state of Roe between ) This represents the average velocity of Roe in the x-direction. Let Roe represent the average speed of sound. It means (q) LU ,q RU The minimum wave velocity along the x-direction of the Roe average state between ) In the multidimensional wave propagation model, the two-dimensional Riemann flux is represented as: Among them, f LD ,f LU ,f RD ,f RU The multidimensional wave propagation model represents the flux in the x-direction at the lower left, upper left, lower right, and upper right corners, respectively, where g is the flux. LD ,g LU ,g RD ,g RU These represent the flux in the y-direction at the lower left, upper left, lower right, and upper right corners of the multidimensional wave propagation model, respectively. The flux at the third Gaussian point G3 It was obtained using the same method as Gaussian point G1; For the second Gaussian point G2, the flux is calculated using a one-dimensional Riemann solver: Among them, f L ,f R These represent the flux on the left and right sides, respectively. These represent the reconstructed values ​​on the left and right sides of the Gaussian point after constraints, respectively. The superscript "m" represents the physical quantity related to the midpoint of the interface. and The maximum wave speeds representing leftward and rightward propagation in the respective wave propagation models are defined as follows: In the formula, 'a' represents the local speed of sound, and the superscript "~" indicates the Roe average. Finally, the interface flux f is obtained by applying Gaussian integral. i+1 / 2,j :

10. The high-precision numerical simulation method for hypersonic vehicles according to claim 1, characterized in that: Step 6 includes the following steps: The semi-discrete finite volume scheme is transformed into a fully discrete finite volume scheme in spacetime using the third-order Runge-Kutta discretization formula for the time variable: q (1) =q n +ΔtL(q n ) L represents the residual, and the superscript "n" represents the time step. The flow field variable values ​​at the next time step are solved using a spatiotemporally discrete finite volume scheme. This process is repeated to obtain the numerical simulation results of the full flow field stability of the hypersonic vehicle.

Images