Unsteady Shock Wave Interaction Flow Field and Aerodynamic Heat Prediction Method under High Mach Number Conditions

Abstract

The present invention relates to a method for predicting the unsteady shock interference flow field and aerodynamic heat under high Mach number conditions for an internal-external flow integrated aircraft. To solve the problem of accurately predicting the unsteady flow field and aerodynamic heat of an internal-external flow integrated aircraft, a numerical simulation analysis method is formed by combining an unsteady time discretization method with a low-dissipation spatial discretization scheme: 1. The delayed detached eddy simulation (IDDES) method is used to numerically simulate the characteristics of the interference unsteady flow field of the internal-external flow integrated aircraft, and flow field parameter results with higher accuracy than the time-averaged method are obtained; 2. The second-order accurate low-dissipation flux splitting spatial discretization scheme (LDFSS-M) is used to numerically simulate the aerodynamic heat trend of the internal-external flow integrated aircraft, and a more accurate surface heat flux distribution trend than the ASUM+ and Roe schemes is obtained.

Description

Technical Field

[0001] The present invention relates to a method for predicting shock-stabilized aerodynamic heat, belonging to the technical field of aerodynamic heat design of air-breathing combined power aircraft, and is mainly used to improve the accuracy of predicting the flow field and aerodynamic heat in the local shock interference area of the existing numerical simulation method. Background Technique

[0002] To meet the requirements of lift augmentation and drag reduction, air-breathing combined power aircraft mostly adopt an integrated internal and external flow design. There are various complex shock / shock and shock / boundary layer interferences in the flow field, and the unsteady flow characteristics are significant. The unsteady flow leads to a sharp increase in pressure and temperature, which is the largest load faced by the aircraft thermal protection system.

[0003] Currently, the research on predicting the interference flow field and aerodynamic heat between shock waves and between shock waves and boundary layers mainly focuses on the high-precision discretization format of the NS equation, turbulence models, etc. The analysis of weak shock interference and the wave system after reflection in the boundary layer can be completed by theoretical analysis and high-order accurate numerical methods of Reynolds averaging. The simulation results are also relatively consistent with the wind tunnel test results. However, the shock interference of the integrated internal and external flow aircraft is unstable, and the flow field parameters change violently over time. The unsteady flow field contains pulsations of various parameters such as pressure, temperature, and density, resulting in a large gradient of the thermal load distribution in the interference area, and it is very difficult to accurately predict the flow field parameters and heat flux in the interference area. Summary of the Invention

[0004] The technical problem to be solved by the present invention is: to solve the problem of accurately predicting the unsteady flow field and aerodynamic heat of the integrated internal and external flow aircraft, a method for predicting the unsteady shock interference flow field and aerodynamic heat under high Mach number conditions of the integrated internal and external flow aircraft is formed by combining the unsteady time discretization method with a low-dissipation spatial discretization format.

[0005] The technical solution of the present invention is: a method for predicting the unsteady shock interference flow field and aerodynamic heat under high Mach number conditions of an integrated internal and external flow aircraft, including:

[0006] Based on the numerical simulation finite volume method, a calculation framework of the NS equation suitable for predicting unsteady shock interference is built;

[0007] According to the Delayed Detached Eddy Simulation (IDDES) method, by solving the scale function L IDDES , the near-wall and far-wall regions are divided, and the calculation methods of the NS equation calculation framework in the two regions are determined;

[0008] The flux splitting spatial discretization format LDFSS-M is used to spatially discretize the NS equation calculation framework. During the spatial discretization process, the problem of too small numerical viscosity near the shock wave and too large numerical viscosity in the boundary layer during shock capture is solved by redefining the numerical sound speed, and shock capture is realized.

[0009] Using the completed discrete NS equation calculation framework to predict the unsteady shock interference flow field and aerodynamic heat under high Mach number conditions of the aircraft.

[0010] Preferably, the prediction results are verified through simulation examples and wind tunnel test results, and the prediction process is corrected when there are deviations.

[0011] Preferably, the simulation examples include NACA0012 airfoil and HYFLEX aircraft simulation examples, and the Ma6 shock wind tunnel is used for the wind tunnel test. [[ID=[]]

[0012] Preferably, the NS equation calculation framework processes the NS equation using a third-order accurate WENO reconstruction method, and the time discretization format of the NS equation calculation framework adopts a third-order TVD Runge-Kutta format.

[0013] Preferably, a non-linear weight is used to replace the smooth factor in the WENO reconstruction method to avoid the reduction of accuracy caused by oscillation in the non-smooth region; the formula for the non-linear weight is as follows:

[0014]

[0015] In the above formula, γ i，j is the linear weight, ω i，j is the non-linear weight, IS i，j is the smooth factor, and ε is a fixed parameter introduced to prevent the denominator from being zero in the WENO reconstruction method; i and j represent the rows and columns of the coefficient matrix in the WENO reconstruction method.

[0016] Preferably, an adaptive function is used to replace the fixed parameter ε to improve the accuracy of the extreme points; the adaptive function ε is a function of the smooth factor, and the expression is:

[0017]

[0018] In the formula, δ is a minimum value taken to prevent the denominator from being zero.

[0019] Preferably, during the construction of the NS equation calculation framework, since there is excessive numerical dissipation when simulating low Mach number flows, the lower interface pressure term of the NS equation calculation framework is processed by scaling the velocity diffusion term to avoid excessive numerical dissipation; the processed lower interface pressure term p s is:

[0020] f(M) is a scaling function, and its definition is f(M) = min(1.0, max(M, M ref)); In the above formula, the subscripts "L" and "R" respectively indicate that the variable is obtained by interpolation through the limiter on the left and right sides of the unit interface, D + represents the forward difference value at the unit interface, D - represents the backward difference value at the unit interface; M is the local Mach number, and its value is related to the Mach number |M L | on the left side of the interface and the Mach number |M R | on the right side. M = min(1.0, max(|M L |, |M R |)), M ref is the Mach reference value, obtained by solving the NS equation; p is the unit pressure value.

[0021] Preferably, the RANS method is used for calculation in the near-wall region, and the LES method is used for calculation in the region far from the wall, so that a smooth transition can be achieved in the shock interference region during calculation, and higher-precision flow field parameters can be obtained.

[0022] Preferably, the numerical speed of sound a redefined by the AUSMPW+ format i+1 / 2 :

[0023]

[0024]

[0025]

[0026] where U is the integral average value of the compressible NS equation in the calculation cell, V is the velocity of the calculation cell, and the subscripts "L" and "R" respectively indicate that the variable is on the left and right sides of the unit interface; H total,L and H total,R are the total specific enthalpies obtained by interpolation through the limiter on the left and right sides of the calculation cell interface respectively, and γ represents the specific heat of the gas.

[0027] Preferably, the shock anomaly phenomenon that occurs when calculating hypersonic blunt bodies is solved by introducing multi-dimensional numerical dissipation, that is, the multi-dimensional numerical dissipation SVM is used to replace the dissipation SV in LDFSS-M; the formula for the multi-dimensional numerical dissipation SVM is as follows:

[0028]

[0029]

[0030]

[0031] where ΔU represents the difference in the integral average value dissipation before and after in the calculation cell, Δu and Δv are the differences in the dissipation of the X-direction and Y-direction velocity components before and after respectively, n x and n yThey are the X - and Y - component vectors of the unit vector n respectively, a represents the sound speed in the element, ρ represents the density in the element, and p Lk , p Rk Interpolate through the limiter to calculate the pressures at the boundaries where the unit interfaces with the k - th adjacent unit on the left and right sides of the unit interface.

[0032] The advantages of the present invention compared with the prior art are as follows:

[0033] The present invention uses the delayed detached eddy simulation (IDDES) method to numerically simulate the unsteady flow field characteristics of the interference of an internal - external flow integrated aircraft, obtaining more accurate flow field parameter results compared with the time - average method; uses the second - order accurate low - dissipation flux - splitting spatial discretization scheme (LDFSS - M) to numerically simulate the aerodynamic heat trend of the internal - external flow integrated aircraft, obtaining a more accurate surface heat flux distribution trend compared with the ASUM + and Roe schemes. Specifically:

[0034] (1) The present invention innovatively proposes and implements a new aerodynamic heat prediction method that can meet the requirements of high - Mach shock stability and boundary - layer solution. After careful comparison, the IDDES method has a pressure prediction error of about 25% in the shock oscillation region of the interference of the internal - external flow integrated aircraft, and the RANS method has a pressure prediction error of about 50% in the shock oscillation region. The IDDES method has better accuracy in the simulation of unsteady shock oscillation flow field parameters;

[0035] (2) By using the LDFSS - M scheme to correct the aerodynamic heat prediction of the internal - external flow integrated aircraft, under the same grid conditions, the LDFSS - M scheme has a shock oscillation heat flux prediction error of 26.99% in the interference region of the internal - external flow integrated aircraft, and the Roe's FDS scheme has a heat flux prediction error of 45%. The LDFSS - M scheme has better simulation accuracy in the heat flux prediction of the interference region. Brief Description of the Drawings

[0036] Figure 1 It is a comparison diagram of heat flux results before and after the improved method;

[0037] Figure 2 It is the simulation result of the wall heat flux on the windward side of the HYFLEX aircraft by the method of the present invention. Detailed Embodiment

[0038] A method for predicting the unsteady shock interference flow field and aerodynamic heat of an internal - external flow integrated aircraft under high - Mach conditions includes the following steps:

[0039] S1, based on the existing numerical simulation WENO finite - volume method, build an NS equation calculation framework suitable for unsteady shock interference prediction.

[0040] (1) NS calculation framework

[0041] Ignoring the effects of heat transfer, viscosity, and chemical reactions, the NS equations for compressible fluids are the two-dimensional Euler equations, as follows:

[0042]

[0043] where

[0044]

[0045]

[0046] In Equation (2-2), ρ, u, v, p, and E represent the density of the fluid, the velocity component in the x-direction, the velocity component in the y-direction, the pressure, and the total energy per unit volume, respectively. To close the governing equations, the equation of state of the gas needs to be supplemented. Assuming the gas satisfies the ideal gas assumption, the expression for the equation of state of the gas is

[0047]

[0048] γ represents the specific heat of the gas;

[0049] Integrating the NS equations for compressible fluids over the control volume gives the following integral form:

[0050]

[0051] In the above equation, Ω i is the area of the i-th control volume, n is the unit outer normal vector, F = (F, G), and U i is the average value of the conserved variables within the i-th control volume:

[0052]

[0053] The Gaussian quadrature method is used to calculate the fluxes on the control volume boundary. On the k-th edge:

[0054]

[0055] In the above equation, n k represents the unit outer normal vector of the k-th edge, U G represents the physical quantity at the Gaussian point, q represents the number of Gaussian integration points, and |Γ k | represents the length of the k-th edge; when reconstructing a quadratic polynomial, q is equal to 2; since the reconstructed values on the left and right sides of the boundary are inconsistent, a numerical flux is needed:

[0056]

[0057] Thus, the NS equation calculation framework is as follows:

[0058]

[0059]

[0060] The time discretization format of the NS equation calculation framework adopts the third-order TVD Runge-Kutta format, as follows:

[0061]

[0062] (2) WENO reconstruction method

[0063] The present invention uses the third-order precision WENO method to reconstruct the physical quantity at the Gaussian point. The reconstruction process involves 10 units including the current unit. First, a quadratic polynomial Q(x,y) is constructed, which is in the form

[0064]

[0065]

[0066] Where (x0, y0) is the coordinate of the current cell's centroid. To solve the coefficients of the quadratic polynomial, the integral of Q(x, y) over all cells must be equal to the average value of that cell. Subtracting the average value of all cells in the template from the average value of the current cell yields a set of equations. For example, for cells i and O, the following equation is used:

[0067]

[0068]

[0069] Writing all nine equations in matrix form gives:

[0070]

[0071] Formula (2-14) has 5 unknowns and 9 equations, which can be solved by transposing the least squares method:

[0072] A T Aφ=A T ψ (2-16)

[0073] in,

[0074]

[0075] To increase the stability of the algorithm, this paper multiplies each row of the coefficient matrix by a weight. The weight is equal to the inverse of the distance from the centroid of the calculation unit to the centroid of the central unit, and its specific form is

[0076]

[0077] After adopting the weighted least squares, the original equation can be written as

[0078]

[0079] where the expression of the coefficient matrix is

[0080]

[0081] After reconstructing the quadratic polynomial Q(x, y) using the k-exact method, and using the idea of the WENO method, the entire template is divided into 9 sub-templates: S1 = {O, i, j}, S2 = {O, i, k}, S3 = {O, j, k}, S4 = {O, i, i1}, S5 = {O, i, i2}, S6 = {O, j, j1}, S7 = {O, j, j2}, S8 = {O, k, k1}, S9 = {O, k, k2}. Within each template, a linear polynomial in the following form is defined:

[0082] p(x, y) = U0 + β1ξ + β2η (2-21)

[0083] For each Gauss point, find a set of non-negative linear weights γ j such that the 9 linear polynomials satisfy

[0084]

[0085]

[0086] Using the element U j Equating the coefficients on both sides of the equation gives a system of equations for γ j

[0087]

[0088] In equation (2-23), the elements a i,j in the coefficient matrix and λ i on the right-hand side of the equal sign are both related to the grid, and a i,j = f(β i ). The linear weights of the sub-template can be obtained through equations (2-22) and (2-23). The last 8 equations of (2-23) can be written in the following form:

[0089]

[0090] To ensure computational stability, this paper derives the algebraic expression of the inverse matrix of the coefficient matrix of equation (2-24):

[0091]

[0092] ​Among them, the expression of B is:

[0093]

[0094] This avoids the operation of matrix inversion in the program and improves the robustness and computational efficiency of the algorithm.

[0095] (3) Nonlinear Weight Suppression of Discontinuous Oscillation Method

[0096] In the smooth region, the nonlinear weight approaches the linear weight, and third-order accuracy can be achieved; in the region containing discontinuities, in order to ensure no oscillation, the WENO reconstruction method will convert the third-order accuracy to second-order accuracy. In order to suppress oscillation near the discontinuity in the present invention, it is necessary to use the smooth factor in the nonlinear weight replacement method. The process of solving the nonlinear weight using the linear weight and the smooth factor is as follows:

[0097]

[0098] In the above formula, γ i，j is the linear weight, ω i，j is the nonlinear weight, IS i，j is the smooth factor, ε is a fixed parameter introduced to prevent the denominator from being zero in the WENO reconstruction method; i, j represent the rows and columns of the coefficient matrix in the WENO reconstruction method. In the literature, ε is generally taken as a positive number between 10 -2 and 10 -6 .

[0099] At the same time, in order to avoid that the third-order accuracy WENO format can only reach the second-order accuracy at the extreme points, this paper proposes to use an adaptive function to replace the fixed small parameter method to improve the accuracy of the extreme points. That is, in formula (3-1), an adaptive function is used to replace the original fixed small parameter ε. The adaptive function is a function of the smooth factor, and its expression is:

[0100]

[0101] Among them, δ is taken as 10 -40 to prevent the denominator from being zero and does not change with the problem being calculated.

[0102] (4) There is excessive numerical dissipation when simulating low Mach number flows, and the velocity diffusion term is scaled

[0103] The interface pressure term and the expanded formula under the NS equation framework are as follows:

[0104]

[0105]

[0106] To avoid excessive numerical dissipation, the high-order terms of Ma are ignored under subsonic conditions:

[0107]

[0108] where f(M) is a scaling function, defined as

[0109] f(M) = min(1.0, max(M, M ref )) M = min(1.0, max(|M L |, |M R |)) (4-4)

[0110] S2. According to the delayed detached eddy simulation IDDES method, the near-wall and far-wall regions are divided by solving the scale function L IDDES . Based on the NS equation calculation framework, RANS calculation is used in the boundary layer of the near-wall region, and LES calculation is used in the far-wall region.

[0111] The IDDES method is achieved by modifying the length scale function in the destruction term of the RANS turbulent kinetic energy transport equation. In the present invention, the SST k-ω model is adopted, and its turbulent kinetic energy equation is:

[0112]

[0113] In the formula, σ k is a constant; the destruction term D k = β * ρωk, β * = 0.09. Since the scale function in the RANS equation Therefore, it can be obtained

[0114]

[0115] In the present invention, the IDDES method is adopted, and the assignment takes the form of:

[0116]

[0117] L RANS = k 12 / (β * ω), L LES = C DES Δ

[0118]

[0119]

[0120]

[0121] where vt represents the turbulent kinematic viscosity coefficient, К is the Karman constant, taking the value of 0.41. C dt is an empirical constant, taking the value of 20 in the SST model of the present invention. f B is defined as follows:

[0122] f B = min{2exp(-9α 2 ), 1.0}

[0123] where α = 0.25 - d w / h max .

[0124] Δ is the grid scale and is defined as follows:

[0125] Δ = min{max[C w d w , C w h max , h min , h max}

[0126] h max = max(Δ x , Δ y , Δ z )

[0127] h min = min(Δ x , Δ y , Δ z )

[0128] Here, d w is the distance from the wall, and C w is an empirical constant, taking the value of 0.15 in the present invention.

[0129] When performing unsteady shock interference calculations, within the boundary layer Therefore, L IDDES = L RANS , and the flow field calculation uses the RANS calculation; in the region far from the wall r dt << 1, f dt ≈ 1, Only in the region where d w / h max < 0.5, L IDDES = L RANS , and the flow field calculation uses the RANS calculation; while in the region where d w / h max > 1.0, L IDDES completely transitions to L LES , and the flow field calculation uses the LES calculation.

[0130] S3. Based on the calculation framework in step S2, a LDFSS-M numerical format is proposed to solve the problem of low shock-capturing resolution in existing prediction methods.

[0131] (1) Establishment of the LDFSS-M format:

[0132] The numerical sound speed definition method adopted by the LDFSS format is arithmetic mean, which makes LDFSS unable to capture shocks acutely. Therefore, a new numerical sound speed a is introduced in the LDFSS-M format i+1 / 2 Calculation formula

[0133]

[0134]

[0135]

[0136] where U is the integral average value of the compressible NS equation in the computational cell, V is the velocity, and γ represents the specific heat of the gas. H total,L and H total,R are the total specific enthalpies interpolated through the limiter on the left and right sides of the computational cell interface respectively, where C p,L and C p,R are the pressure coefficients interpolated through the limiter on the left and right sides of the computational cell interface respectively, T L and T R are the temperatures interpolated through the limiter on the left and right sides of the computational cell interface respectively, u L and u R are the X-direction velocity components interpolated through the limiter on the left and right sides of the computational cell interface respectively.

[0137] (2) When calculating hypersonic blunt bodies, shock anomalies are likely to occur, and multi-dimensional numerical dissipation needs to be introduced.

[0138] The multi-dimensional numerical dissipation SVM formula is as follows

[0139]

[0140]

[0141]

[0142] where ΔU represents the difference in the integral average value in the computational cell before and after dissipation, Δu and Δv are the differences in the velocity components in the X and Y directions before and after dissipation, n x and n yThey are the X and Y components of the unit vector n respectively. a represents the sound speed in the element, ρ represents the density in the element, and the pressure-sensitive basis function g is as shown in the above formula. Suppose this element is connected to k adjacent elements, then its boundary can be divided into k boundaries connected to the adjacent elements, p Lk , p Rk Interpolate through the limiter to calculate the pressures at the boundaries where this element is connected to the k-th adjacent element on the left and right sides of the element interface.

[0143] Using the discrete NS equation calculation framework obtained after the above process and combining with specific regional calculation methods, that is, it constitutes a method for predicting the unsteady shock interference flow field and aerodynamic heat under high Mach number conditions for an integrated internal and external flow aircraft. The prediction method is verified and corrected through standard models such as typical airfoils and high-speed aircraft and the test results of the Ma6 shock tunnel;

[0144] Among them:

[0145] (1) The NACA0012 airfoil verification model verifies the suppression of oscillations near the shock discontinuity by the prediction method of the present invention; Figure 1 In the left figure, the heat flux calculation result obtained without using the prediction method of the present invention shows obvious anomalies, while the oscillations near the shock discontinuity are suppressed by the prediction method of the present invention, and a smooth heat flux distribution is obtained ( Figure 1 right figure).

[0146] (2) The HYFLEX aircraft verification model verifies the capture of shock discontinuities and the elimination of shock instability phenomena by the prediction method of the present invention; The original format generates oscillations at the head heat flux, showing a non-physical heat flux distribution. And as Figure 2 shown, the LDFSS-M format of the present invention obtains a reasonably distributed heat flux in the high-heat area at the head.

[0147] (3) For the unsteady flow region dominated by complex wave systems such as shock / shock and shock / boundary layer interference, by developing the IDDES method, while improving the shock capture resolution, the advantage of calculation efficiency is obtained;

[0148] (4) For the prediction of local peak heat fluxes in local complex flows, such as irregular shapes and V-shaped overflow ports, develop the LDFSS-M format, and effectively avoid shock anomalies by strengthening the shock capture stability and high-resolution characteristics of the boundary layer.

[0149] By combining theory with actual engineering applications, finally, the method for predicting the unsteady shock interference flow field and aerodynamic heat under high Mach number conditions for an integrated internal and external flow aircraft is applied to the prediction of the shock interference flow field of high-speed aircraft.

[0150] The present invention relates to a method for predicting the unsteady shock interference flow field and aerodynamic heat under high Mach number conditions for an internal-external flow integrated aircraft, which provides accurate distribution characteristics of the aerodynamic heat environment in the complex flow region of shock interference and technical support for hypersonic aircraft in the conceptual design stage. At the same time, this shock prediction method can be applied to the thermal environment design of similar aircraft, providing a basis and technical support for the shape optimization of hypersonic aircraft with complex shapes, the design of thermal protection system solutions, and material selection.

[0151] The above is the best specific implementation mode of the present invention, but the protection scope of the present invention is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed by the present invention should be covered within the protection scope of the present invention.

[0152] The content not described in detail in the specification of the present invention belongs to the well-known technology of those skilled in the art.

Claims

1. A method for predicting the unsteady shock interference flow field and aerodynamic heat under high Mach number conditions of an internal-external flow integrated aircraft, characterized in that Including: Based on the numerical simulation finite volume method, build a NS equation calculation framework suitable for predicting unsteady shock wave interference; According to the IDDES method of delayed detached eddy simulation, by solving the scale function L IDDES , the near-wall region and the far-wall region are divided, and the calculation methods of the NS equation calculation framework in the two regions are determined; Use the flux splitting spatial discretization format LDFSS-M to discretize the NS equation calculation framework. During the spatial discretization process, solve the problem that the numerical viscosity near the shock wave is too small during shock wave capture and too large within the boundary layer by redefining the numerical sound speed, and achieve shock wave capture; Use the NS equation calculation framework that has completed the above discretization to predict the unsteady shock wave interference flow field and aerodynamic heat under high Mach number conditions of the aircraft; The NS equation calculation framework is processed by using the third-order accurate WENO reconstruction method for the NS equation, and the time discretization format of the NS equation calculation framework adopts the third-order TVD Runge-Kutta format; Use non-linear weights to replace the smooth factor in the WENO reconstruction method to avoid the reduction of accuracy caused by oscillation in the non-smooth area; the formula for the non-linear weights is as follows: In the above formula, γ i，j is the linear weight, ω i，j is the non-linear weight, IS i，j is the smoothing factor, and ε is a fixed parameter introduced to prevent the denominator from being zero in the WENO reconstruction method; i and j represent the rows and columns of the coefficient matrix in the WENO reconstruction method; Replace the fixed parameter ε with an adaptive function to improve the accuracy of the extreme point; the adaptive function ε is a function of the smooth factor, and the expression is: In the formula, δ is a minimum value taken to prevent the denominator from being zero.

2. The method according to claim 1, wherein: Verify the prediction results through simulation examples and wind tunnel test results, and complete the correction of the prediction process when there are deviations.

3. The method according to claim 2, characterized in that: The simulation examples include the NACA0012 airfoil and HYFLEX aircraft simulation examples, and the Ma6 shock tunnel is used for the wind tunnel test.

4. The method according to claim 1, wherein: During the process of building the NS equation calculation framework, due to excessive numerical dissipation when simulating low Mach number flows, the interface pressure term in the NS equation calculation framework is processed by scaling the velocity diffusion term to avoid excessive numerical dissipation; the processed lower interface pressure term p s is as follows: f(M) is a scaling function, which is defined as f(M) = min(1.0, max(M, M ref )); in the above formula, the subscripts "L" and "R" respectively indicate that the variable is obtained by interpolating through the limiter on the left and right sides of the cell interface, D + represents the forward difference value at the cell interface, D - represents the backward difference value at the cell interface; M is the local Mach number, and its value is related to the Mach number |M L | on the left side of the interface and the Mach number |M R | on the right side. M = min(1.0, max(|M L |, |M R |)), M ref is the Mach reference value, obtained by solving the NS equation; p is the cell pressure value.

5. The method according to claim 1, wherein: Use the RANS method to calculate in the near-wall region and the LES method to calculate in the region far from the wall, so that smooth transition can be achieved in the shock wave interference area during calculation, and higher-precision flow field parameters can be obtained.

6. The method according to claim 1, characterized in that: The numerical sound speed a redefined using the AUSMPW+ format i+1 / 2 : Where, U is the integral average value of the compressible NS equation in the computational cell, V is the velocity of the computational cell, and the subscripts "L" and "R" respectively indicate that the variable is on the left and right sides of the cell interface; H total,L and H total,R are respectively the specific total enthalpies obtained by limiter interpolation on the left and right sides of the computational cell interface, and γ represents the specific heat of the gas.

7. The method according to claim 1, wherein: Solve the abnormal shock wave phenomenon that appears when calculating hypersonic blunt bodies by introducing multi-dimensional numerical dissipation, that is, use the multi-dimensional numerical dissipation SVM to replace the dissipation SV in LDFSS-M; the formula for the multi-dimensional numerical dissipation SVM is as follows: Among them, ΔU represents the difference in the integral average dissipation before and after in the calculation unit, Δu and Δv are the differences in the X-direction and Y-direction velocity components before and after dissipation respectively, n x and n y are the X-direction and Y-direction components of the unit vector n respectively, a represents the sound speed in the unit, ρ represents the density in the unit, p Lk , p Rk The left and right sides of the interface of the calculation unit obtain the pressures at the boundaries where the unit is connected to the k-th adjacent unit through limiter interpolation.

Images