Simulation parameter optimization method for numerical simulation of hypersonic velocity external flow field

Abstract

The invention discloses a simulation parameter optimization method for numerical simulation of a hypersonic velocity external flow field, which comprises the following steps: selecting a sharp trailing edge NACA0012 airfoil profile, establishing an airfoil profile mathematical model, and generating an airfoil profile geometric model; carrying out grid division by adopting ICEM CFD, and setting the height of a first layer of grid units to enable y + to be less than or equal to 1; setting an expansion layer and a grid growth rate by adopting an expansion layer grid strategy; setting a Reynolds number of a near shock wave surface grid unit, and determining a grid division scheme; carrying out hypersonic flow calculation by adopting a density basis solver; carrying out numerical calculation by using an ROE flux splitting method; performing flow field variable interpolation calculation by adopting a second-order windward format, and performing gradient calculation based on a least square method; and performing error evaluation on a simulation result through wind tunnel test data. According to the method, a definition formula modeling method is adopted, direct use of a database for interpolation modeling is avoided, and therefore the smoothness and accuracy of the geometric curve of the airfoil profile are improved.

Description

Technical Field

[0001] The invention relates to the technical field of hypersonic external flow field sampling, and in particular to a simulation parameter optimization method for numerical simulation of a hypersonic external flow field. Background Art

[0002] During the flight of a hypersonic vehicle, the airflow around it presents complex phenomena such as shock wave-boundary layer interaction, strong nonlinear turbulence, aerodynamic thermal effects, and flow separation. In order to study the aerodynamic characteristics of hypersonic vehicles, wind tunnel experiments are usually combined with computational fluid dynamics (CFD) numerical simulations for analysis. However, due to the high equipment cost, short test time, and difficulty in measuring complex flow field parameters in hypersonic wind tunnel experiments, the use of high-precision CFD numerical simulation methods to predict hypersonic flow phenomena has become an important direction in the study of hypersonic aerodynamics.

[0003] At present, the main difficulties in numerical simulation of hypersonic external flow fields are:

[0004] 1. Insufficient accuracy of airfoil modeling. Existing NACA airfoil modeling is usually based on database or interpolation methods, which may lead to insufficient curve smoothness, thus affecting the accuracy of numerical calculations. In addition, too few or too many data points will affect the numerical stability of the airfoil geometry, thereby affecting the simulation error.

[0005] 2. Insufficient optimization of the far-field boundary conditions in the computational domain and unreasonable setting of the far-field boundary distance may lead to boundary reflection errors and affect the stability of the flow field calculation. Summary of the invention

[0006] The purpose of this section is to summarize some aspects of embodiments of the present invention and briefly introduce some preferred embodiments. Some simplifications or omissions may be made in this section and the specification abstract and the invention title of this application to avoid blurring the purpose of this section, the specification abstract and the invention title, and such simplifications or omissions cannot be used to limit the scope of the present invention.

[0007] In view of the above-mentioned and / or existing problems in the simulation parameter optimization method for numerical simulation of hypersonic outer flow field, the present invention is proposed.

[0008] Therefore, the problem to be solved by the present invention is how to provide a simulation parameter optimization method for numerical simulation of hypersonic external flow field.

[0009] In order to solve the above technical problems, the present invention provides the following technical solutions: a simulation parameter optimization method for numerical simulation of hypersonic external flow field, which comprises: selecting a sharp trailing edge NACA0012 airfoil, establishing an airfoil mathematical model, and generating an airfoil geometric model; setting the calculation domain far field distance to 16 times the airfoil chord length, using the Pressure farfield condition for the input boundary, the Pressure far field condition for the output boundary, and the No-slip isothermal wall condition for the airfoil wall boundary; using ICEMCFD for grid division, setting the height of the first layer of grid cells, so that y + ≤1, adopt the expansion layer grid strategy, set the expansion layer and grid growth rate; set the Reynolds number of the near shock wave surface grid unit, and determine the grid division scheme; use the density-based solver for hypersonic flow calculation; select the ROE flux splitting method for numerical calculation; use the second-order upwind scheme for flow field variable interpolation calculation, and use the least squares method for gradient calculation; use the Spalart-Allmaras turbulence model for turbulence calculation, and set the turbulent viscosity ratio at the input boundary; and use the wind tunnel test data to evaluate the error of the simulation results.

[0010] As a preferred solution of the simulation parameter optimization method for numerical simulation of hypersonic external flow field described in the present invention, wherein: the establishment of the airfoil mathematical model refers to calculating the airfoil curve using the NACA4 definition formula, which is expressed as

[0011] y=±0.5947[0.2983x 1 / 2 -0.1271x-0.3579x 2 +0.292x 3 -0.1052x 4 ]

[0012] Where y is the airfoil thickness distribution and x is the position of the airfoil surface along the chord length.

[0013] As a preferred solution of the simulation parameter optimization method for numerical simulation of hypersonic external flow field described in the present invention, data points are connected through CAD software or Python / Matlab code, an airfoil surface is created, and an airfoil geometric model is generated.

[0014] As a preferred solution of the simulation parameter optimization method for numerical simulation of hypersonic external flow field described in the present invention, when setting the height of the first layer of grid cells, the following formula is used for calculation:

[0015]

[0016] y H =2y p

[0017] In the formula, y H is the height of the first layer grid unit, y + is the dimensionless wall distance parameter, u τ is the friction velocity, ρ is the atmospheric density, and μ is the free stream kinematic viscosity.

[0018] As a preferred solution of the simulation parameter optimization method for numerical simulation of hypersonic external flow field described in the present invention, the expansion layer is set to 30 and the grid growth rate is set to 1.05.

[0019] As a preferred scheme of the simulation parameter optimization method for numerical simulation of hypersonic external flow field described in the present invention, for the sharp trailing edge airfoil, the grid unit Reynolds number is set to 0.469, 0.505 and 0.493 respectively; for the blunt trailing edge airfoil, the minimum value of the grid orthogonality is 0.454, 0.495 and 0.472 respectively.

[0020] As a preferred solution of the simulation parameter optimization method for numerical simulation of hypersonic external flow field described in the present invention, wherein: the grid division scheme is to divide the number of grids into 610,000.

[0021] The beneficial effects of the present invention are as follows: the definition formula modeling method is adopted to avoid the direct use of database interpolation modeling, thereby improving the smoothness and accuracy of the airfoil geometry curve. NACA4 data points (200) are selected to optimize the distribution of airfoil curve data points, improve the geometric accuracy of simulation calculations, and reduce numerical errors. By comparing the effects of 12L, 16L and 20L far field distances on calculation accuracy, it is found that the 16L far field distance can effectively reduce boundary reflection errors and improve calculation stability. The far field boundary condition uses Pressure far field to ensure that the flow field boundary conditions are consistent with the actual flow and improve the accuracy of numerical simulation. Improve the grid division strategy and improve the turbulence analysis capability. BRIEF DESCRIPTION OF THE DRAWINGS

[0022] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the following briefly introduces the drawings required for describing the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without creative labor. Among them:

[0023] Figure 1 The distribution of data points is modeled for six NACA0012 airfoils.

[0024] Figure 2 Schematic diagram of the correlation between three far-field distances and numerical accuracy in the incompressible external flow field.

[0025] Figure 3Schematic diagram of the flow field outside the calculation domain of the blunt trailing edge.

[0026] Figure 4 Schematic diagram of the flow field outside the calculation domain of the sharp trailing edge.

[0027] Figure 5 This is a steady-state heat diffusion analysis diagram.

[0028] Figure 6 This is a flow analysis diagram with a 90° turn.

[0029] Figure 7 Analytical diagram for mesh orthogonality analysis.

[0030] Figure 8 This is the node design scheme diagram at the grid level of the pointed trailing edge airfoil.

[0031] Fig. 9 This is a diagram of the node design scheme at the grid level for the blunt trailing edge airfoil.

[0032] Fig.10 is the numerical error rate of the sharp trailing edge airfoil designed based on NACA4 at three far-field distances.

[0033] Fig.11 Numerical error rate of the blunt trailing edge airfoil designed based on Airfoil tools at three far-field distances.

[0034] Fig.12 Numerical error rates of the sharp trailing edge airfoil designed based on the defined formula at three far-field distances.

[0035] Fig.13 It is a linear upwind difference algorithm.

[0036] Fig.14 It is the least squares difference algorithm.

[0037] Fig.15 This is a comparison chart of the numerical error rate under four grid aspect ratios.

[0038] Fig.16 Optimum amplitude plots for the optimal total mean error rate and the associated numerical accuracy for six NACA0012 airfoils.

[0039] Fig.17 Optimal total mean error rates for far-field distances, turbulence models, and flux types, as well as associated numerical optimization amplitudes. DETAILED DESCRIPTION

[0040] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the specific implementation methods of the present invention are described in detail below in conjunction with the accompanying drawings.

[0041] In the following description, many specific details are set forth to facilitate a full understanding of the present invention, but the present invention may also be implemented in other ways different from those described herein, and those skilled in the art may make similar generalizations without violating the connotation of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below.

[0042] Secondly, the term "one embodiment" or "embodiment" as used herein refers to a specific feature, structure, or characteristic that may be included in at least one implementation of the present invention. The term "in one embodiment" that appears in different places in this specification does not necessarily refer to the same embodiment, nor does it refer to a separate or selective embodiment that is mutually exclusive with other embodiments.

[0043] Example 1

[0044] Reference Figure 1 to Figure 17 , which is the first embodiment of the present invention, provides a simulation parameter optimization method for numerical simulation of hypersonic external flow field, which comprises the following steps:

[0045] Select the NACA0012 airfoil with a sharp trailing edge, establish the airfoil mathematical model, and generate the airfoil geometric model;

[0046] The far field distance of the computational domain is set to 16 times the chord length of the airfoil, the input boundary adopts the Pressure far field condition, the output boundary adopts the Pressure far field condition, and the airfoil wall boundary adopts the No-slip isothermal wall condition;

[0047] ICEM CFD is used for meshing, and the height of the first layer of mesh cells is set so that y + ≤1;

[0048] Adopt the expansion layer grid strategy and set the expansion layer and grid growth rate;

[0049] Set the Reynolds number of the grid unit near the shock wave surface and determine the grid division scheme;

[0050] The density-based solver is used for hypersonic flow calculations;

[0051] The ROE flux splitting method is used for numerical calculation;

[0052] The second-order upwind scheme is used for flow field variable interpolation calculation, and the least square method is used for gradient calculation;

[0053] The Spalart-Allmaras turbulence model is used for turbulence calculation, and the turbulent viscosity ratio is set at the input boundary;

[0054] The error of simulation results is evaluated through wind tunnel test data.

[0055] Specifically, there are two types of trailing edge shapes for NACA0012 airfoils: blunt trailing edge and sharp trailing edge. The design can be achieved through three main tools, namely Airfoil tools, NACA4 Airfoil Generator and definition formula. Formula (1) is the definition formula of NACA0012. The modeling data point on the X-axis is represented by x. By inputting the value of x, the corresponding Y-axis modeling data point y can be calculated. Specifically, Airfoiltools provides 132 X-axis modeling data points, while NACA4Airfoiltools provides 200. Using the X-axis data points provided by these two tools and combining them with the definition formula, the corresponding Y-axis data points can be calculated to design the required NACA0012 airfoil.

[0056] y=±0.5947[0.2983x 1 / 2 -0.1271x-0.3579x 2 +0.292x 3 -0.1052x 4 ]

[0057] Where y is the airfoil thickness distribution and x is the position of the airfoil surface along the chord length.

[0058] In summary, the NACA0012 airfoil has two trailing edge shapes, three modeling methods, and two modeling data point sources. The present invention uses the modeling tool DesignModeler provided by Ansys to establish six NACA0012 models with different configurations. The detailed features of these models, including the shape of their trailing edges, the selected modeling methods, the sources of data points, and the specific number of data points, are all shown in detail in Table 1. By comparing and exploring the specific effects of the trailing edge morphology and its modeling method, the source of data points used in the definition formula, and the number of data points on the accuracy of numerical calculations. The distribution of modeling data points for these six NACA0012 airfoils is shown in the figure below. Figure 1 As shown, L is 1m, and it can be observed that there are obvious position differences between different modeling data points.

[0059] Table 1 Six NACA0012 airfoils established

[0060] NACA0012 airfoil Modeling approach Number of data points for modeling Blunt trailing edge airfoil Airfoiltools 132 data points Sharp trailing edge airfoil NACA4 airfoil generator 200 data points Sharp trailing edge airfoil Defining formulas Using 132 data points provided by Airfoiltools Sharp trailing edge airfoil Defining formulas Double the 132 data points to 264 data points Sharp trailing edge airfoil Defining formulas Using 200 data points provided by NACA4 Sharp trailing edge airfoil Defining formulas Double the 200 data points to 400 data points

[0061] After the airfoil model is established, it is necessary to construct an external flow field in the computational domain suitable for numerical simulation based on the determined far-field distance. Reasonable far-field distance values ​​are very important for numerical simulation. For the NACA0012 airfoil, the official recommendation of Ansys is to set the far-field distance to be in the range of 12 to 20 times the airfoil chord length L

[23] . However, the official only provides a rough reference range, neither specifying an exact value nor analyzing the possible impact of different far-field distances on the accuracy of numerical calculations. In order to explore this factor in more detail, based on three far-field distances of 12L / 16L / 20L, the correlation between the far-field distance and numerical accuracy in the incompressible external flow field is studied. Figure 2 As shown in the figure, the results show that there is an obvious correlation between the far field distance and the numerical accuracy under incompressible environmental conditions. Based on the NACA0012 airfoil, the present invention also uses 12 times, 16 times and 20 times L as the far field distance, and constructs Figure 3 and Figure 4 The flow field outside the computational domain is shown. The correlation between different far-field distances and numerical accuracy under hypersonic external flow field conditions is analyzed, and the similarities and differences with the research conclusions of incompressible external flow fields are compared. The black square in the figure marks the origin of the coordinate system. INLET stands for inputboundary, that is, input boundary. And OUTLET stands for outputboundary, that is, output boundary. Both use the Pressurefarfield boundary condition used to simulate the free flow at infinity. AIRFOIL is the wall boundary (wallboundary), which adopts no-slip and isothermal wall conditions. Considering that the mesh division at the wall of NACA0012 airfoil has a great influence on the numerical calculation, the blocks near the wall of the airfoil are divided twice.

[0062] There is a thin layer called boundary layer near the NACA0012 airfoil, which contains all the viscous effects. The viscous effect makes the fluid layer close to the airfoil surface have no relative sliding with the airfoil surface, resulting in a sharp change in the velocity and temperature of the fluid in a very small range near the airfoil wall, from a large value far away from the wall to a value consistent with the wall, thus forming a significant gradient change in the normal direction of the wall. This requires accurate calculation of the appropriate first layer grid unit height (y H ) value, so that the first layer of grid cells is located within the sub-viscous layer to ensure the accuracy and reliability of numerical calculations. The determination of this value is affected by multiple parameters, including y + , R e , the velocity U of the far-field free flow t , speed of sound C air , the friction coefficient of the wall C f , shear stress τW , friction speed u τ , and the specific distance y between the airfoil wall and the center of the first layer of grid cells p The calculation method is shown in formulas (2) to (8):

[0063]

[0064] U t =C air ×M a (3)

[0065] C air =20.05(T t ) 0.5 (4)

[0066] C f =[2log 10 (R e )-0.65] -2.3 (5)

[0067] τ W =0.5×ρU t 2 C f (6)

[0068]

[0069] Calculate y using formulas (2)-(8) p , and then calculate y according to formula (9) H .

[0070] y H =2y p (9)

[0071] T t =81.2K,P t =576Pa, atmospheric density ρ = 0.0247kg / m3, C air =180.6m / s. Because M a =10, so U t =1806m / s. e =10×10 6 , ρ, L and U t Substitute R e The formula is μ=4.46082×10 -5 Pa.s.y + The initial value of is 1

[24] , and R e ,ρ,U t and + Substituting into formulas (5) to (9), we can calculate yH The initial value is 4.6×10 -5 m. In the calculation process, due to the use of empirical formula, the y H is an estimated value. In order to ensure that the near-wall mesh element y + The value is kept within a reasonable range (not greater than 1) during the whole simulation process, and a series of adjustments are required. + =0.3,y H =1.4×10 -5 m can meet the requirements.

[0072] (2) Number of expansion layers (N) and grid height growth rate (r)

[0073] The boundary layer thickness (δ) is calculated as follows:

[0074]

[0075] Thickness of expansion layer (y T ) is calculated as follows:

[0076]

[0077] is a geometric series formula, which can be rewritten using the following identity:

[0078]

[0079] Therefore, the error between the boundary layer and the expansion layer thickness is:

[0080]

[0081] The value of N is related to the turbulence model used. The present invention uses the Reynolds average turbulence model (RANS), and it is generally recommended that the value of N range from 15 to 30. Therefore, after the recommended value of N is given, formula (13) can be regarded as a function related only to r:

[0082]

[0083] Formula (14) aims to make the boundary layer thickness equal to the expansion layer thickness, that is, to solve f(r) = 0. In order to find this solution using the bisection root-finding algorithm, a reasonable initial guess value for r is needed. Experiments have shown that r values ​​of 1.05 or 1.1 have little effect on the numerical calculation results. However, if the value is increased to 1.2, it will be obvious that the numerical calculation results have changed significantly. Furthermore, when the value of r increases to 1.4, this change becomes more obvious. Therefore, in order to ensure the stability and accuracy of the numerical calculation, the initial value of r should be carefully controlled below 1.1. Finally, N = 30 and r = 1.05 are set. The above parameters can ensure the stable execution of numerical simulation.

[0084] The value of the grid aspect ratio has an impact on both steady-state and transient numerical calculations. First, consider the steady-state case, such as Figure 5 As shown, the black dot represents the center point of the grid. The left side and the bottom are both adiabatic walls. To simplify the analysis, only the heat diffusion effect is considered. 1 、x 2 and x 3 is the coordinate of the center point of the grid, T 1 , T 2 and T 3 is the temperature of the corresponding grid cell. Heat flux Q 12 , Q 13 In relation to the temperature gradient, each grid unit is energy-conserving. According to Fourier's law of heat conduction, the following equation can be obtained, where k is the thermal conductivity, A is the contact surface area of ​​the grid unit and A 13 >A 12 , S is the volume heat source term.

[0085]

[0086] Q 12 +Q 13 =S(17)

[0087] Substituting formula (15) and (16) into formula (17), we can get:

[0088]

[0089] Assumptions Formula (18) can be organized as:

[0090]

[0091] Therefore, the calculation result of temperature change is directly related to the value of the grid aspect ratio. As the value of the grid aspect ratio increases further, the contact area A decreases and the grid spacing increases, which means that the effect of the related diffusion effect is further reduced. Both the steady-state momentum and thermal energy equations contain diffusion terms, and the diffusion terms are related to the value of aspectratio. When the value of aspectratio is large, the value of the related diffusion term will be very small, so the value of aspectratio has a certain influence on the steady-state numerical calculation. The transient situation is discussed below. The transient numerical simulation is limited by the Courant number. Generally, it is necessary to ensure that the Courant number does not exceed 1 to maintain the stability and accuracy of the numerical calculation, that is, within a time step Δt, the movement of the fluid should be smaller than the width of the grid unit. For the boundary layer grid, because y + Due to the limitation of the value, the height of the boundary layer grid unit is very small, resulting in a large corresponding aspectratio value. However, according to the properties of the boundary layer, the fluid velocity near the wall is often very small. In the area outside the boundary layer, as the height of the grid unit increases, the aspectratio value decreases, but at the same time the fluid velocity also increases. Therefore, although there is a large aspectratio in the boundary layer, because the flow velocity is very small, it will not have much impact on the stability of the transient calculation. But if there is such a Figure 6 As shown in the flow with a 90° turn, there is a large aspectratio in this case, which may cause divergence in numerical calculations. Therefore, for transient simulations, if the grid has a large aspectratio value, it is necessary to analyze the position of these units in the grid to determine whether it affects the numerical calculation. This also explains why high aspectratio is a warning message rather than an error message. In summary, for steady-state flow, the value of aspectratio has a direct impact on the diffusion effect. The larger the aspectratio value, the smaller the contact surface, which leads to a lower diffusion effect. In particular, it has a significant impact on the calculation of the pressure correction equation for incompressible flow. For transient flow, if the value of the grid aspectratio is large, a smaller time step Δt may be required, but this is not absolute and depends on the specific position of the high aspectratio grid. The present invention discusses the maximum aspectratio values ​​of the sharp trailing edge airfoil and the blunt trailing edge airfoil at three far-field distances, which are 2100, 2870, 3630 and 2760, 3680, 4770, respectively, which can ensure the stability of numerical simulation.

[0092] There are two cases in grid orthogonality analysis, such as Figure 7As shown in the figure, d is the vector connecting the center points of two mesh units, n is the surface normal vector of the contact surface of two mesh units, c is the vector connecting the center point of the current mesh unit and the center point of the contact surface, and θ is the orthogonal angle for evaluating the orthogonality of the mesh. Therefore, there are two ways to define the orthogonal angle, as shown in formulas (20) and (21):

[0093]

[0094] The maximum value is taken as the index to judge the orthogonality of the grid. When θ = 90°, the orthogonality quality of the grid is 0, which is the worst quality; when θ = 0°, the orthogonality quality of the grid is 1, which is the best quality. The quality of grid orthogonality and the diffusion term of the flow field control equation The discretization of is closely related to the discretization of . According to the divergence theorem, the volume integral is converted into the surface integral, as shown in formula (22):

[0095]

[0096] v fi is the value of the kinematic viscosity of the corresponding grid unit surface, which can be calculated by the linear upwind difference scheme, near the shock surface k:

[0097] n=Δ+k (23)

[0098] Substituting formula (23) into formula (22), we can get

[0099]

[0100] Formula (24) contains the part of the orthogonal quantity Among them U P Represents the velocity variable of the current grid cell, U N represents the velocity variable of the adjacent grid unit. Both of them are unknown quantities, so the part containing orthogonal quantities is also called implicit quantity. Take the numerical calculation of the last iteration to calculate the non-orthogonal quantity Therefore, it is also called the display quantity. The orthogonal angle θ is positively correlated with the non-orthogonal quantity k. The larger θ is, the larger the display quantity including the non-orthogonal quantity will be, which will reduce the stability of the numerical calculation and increase the possibility of divergence. If the mesh quality is poor, Ansysfluent will use non-orthogonal corrector loops for correction, and perform an inner loop of pressure value correction to ensure the convergence of pressure values. This will undoubtedly consume more computing resources, so a better mesh should be divided as much as possible. According to the simulation test results, for the pointed trailing edge airfoil, the minimum values ​​of the mesh orthogonality at the three far-field distances are 0.469, 0.505 and 0.493 respectively. For the blunt trailing edge airfoil, the minimum values ​​of the mesh orthogonality are 0.454, 0.495 and 0.472 respectively. The above values ​​can ensure the stable operation of the numerical simulation.

[0101] To ensure the credibility of the numerical simulation results, it is necessary to refine the grid to different degrees and verify whether the effect of grid refinement on the results is independent, that is, grid independence verification. This means that it is necessary to compare the error rate between the simulation results and the experimental data, and at the same time comprehensively consider the currently available computing resources to determine the most appropriate number of grids. In the present invention, three grid independence analyses of different fineness are implemented, and the specific numbers of grids are about 440,000, 610,000 and 850,000 respectively. The node designs of the three grid levels are as follows: Figure 8 and Fig. 9 As shown in the figure, the number of nodes in each direction increases by about 20% each time the refinement is made. Taking the sharp trailing edge airfoil designed based on NACA4, combined with the SSTk-omega turbulence model and ROE flux type as an example, the error rate between the simulation data and the wind tunnel test data at the three far-field distances of the above three grid levels is calculated. Fig.10 As shown, P / P at 12L far field distance and 3 grid levels t 、T / T t and U / U t The average numerical error rates between the sampling locations and the wind tunnel test data are (3.094% 5.608% 2.188%), (2.277% 3.964% 1.034%) and (2.248% 3.936% 1.025%), respectively. t T / T t U / U t ) are 3.630%, 2.425%, and 2.403% respectively. Similarly, at 16L far field distance, P / P t 、T / T t and U / U t The average numerical error rate and the total average error rate are (11.570% 5.638% 3.038%), (6.229% 3.251% 2.356%), (6.212% 3.218% 2.338%) and 6.749%, 3.945%, 3.923% respectively. At 20L far field distance, P / P at three grid levels t 、T / T t and U / U tThe average numerical error rate and the total average error rate are (10.198% 5.949% 2.835%), (3.151% 4.140% 2.277%), (3.141% 4.106% 2.249%) and 6.327%, 3.189%, 3.165% respectively. At the three far-field distances, as the number of grids increases from 610,000 to 850,000, there is almost no further optimization of the numerical error rate. The numerical error rates of the airfoils designed based on the other two Airfoiltools and defined formula modeling methods at the three far-field distances are as follows Fig.11 and 12 As shown, P / P at three grid levels t 、T / T t and U / U t The numerical error rate also has a similar situation. Under the three far-field distances, the further increase in the number of grids hardly brings about obvious optimization of calculation accuracy, and the numerical calculation converges. After balancing the time required for numerical calculation and the available hardware resources, the present invention selects a partitioning scheme containing 610,000 grid cells.

[0102] AnsysFluent provides users with two flow field solver options: pressure-based solver and density-based solver. a For low velocity incompressible / incompressible flows with M < 3, a pressure-based solver can be used. In contrast, a density-based solver is better at handling M a ≥3 high-speed compressible flow. In view of the fact that the present invention focuses on the discussion of hypersonic external flow fields, a density-based solver is selected for calculation. ROE and AUSM are the two main differential algorithms used for inviscid flux vectors. The present invention takes the NACA0012 airfoil as the feature object and also uses the above two flux types to study the correlation between the above two flux types and the prediction accuracy in the hypersonic external flow field. In the hypersonic external flow field, it is necessary to solve the mass, momentum and energy equations, and calculate the values ​​of the flow field variables at the center of the grid surface. Ansysfluent stores the flow field variables at the center of the grid unit by default, and there are three methods: linear / center difference, upwind difference and linear upwind difference. Linear difference has good second-order accuracy, but it may also produce unbounded solutions and non-physical oscillations, which in turn lead to instability problems in numerical calculations. The upwind difference method introduces the direction of mass flux to avoid the generation of unbounded solutions and oscillations and ensure the stability of numerical calculations, but this method has only first-order accuracy. Linear upwind difference is based on upwind difference and introduces gradient to improve the accuracy of the interpolation algorithm.

[0103] The variable varies linearly between the center of the grid cell and the center of the surface, so this method has second-order accuracy and is also called second-order upwind in Ansysfluent. Fig.13 As shown, where φ 1 ,φ2 and φ f Represent the variables of the current grid cell, the adjacent grid cell and the center of the contact surface, respectively. 1 , X 2 and X f are the corresponding coordinates, and represents the gradient at the center of the grid. The calculation equation is shown in formula (25), where r is the distance vector between the center of the grid unit and the center of the surface, and F f is the mass flux direction. However, this method may lead to local maxima or minima, so a gradient limiter is needed to limit the gradient, as shown in formula (26), where is a parameter used to constrain the gradient, Gradient at the center of the grid cell There are three main algorithms: the grid-based Green-Gauss method and the grid-based least squares method. The first two are derived based on the divergence theorem, and the only difference is how to calculate the value of the center of the grid surface, while the least squares method does not use the value of the center of the grid surface. Fig.14 As shown, red represents the grid cell φ where the gradient is to be calculated 1 , green represents the neighboring grid cell φ 2 ~φ 5 , d represents the distance between the grids, as shown in formula (27). This equation can be rearranged as form, where d 1 , and φ 1 As shown in formula (28), the grid unit has 4 faces, so d 1 is a 4*3 matrix. Suppose the grid unit has M faces, then d 1 is an M*3 matrix, is a 3*1 matrix, φ 1 is an M*1 matrix. Generally speaking, M>3 means that the solution There is only an approximate solution, and the least squares method is used to minimize the sum of squares of the errors, so this method is called the least-squares cell-based interpolation method. Its solution is shown in formula (29), where φ c Represents the value at the center of the grid cell to be calculated, φ s Represents the value at the center of the surrounding grid cells. T d is a 3*3 matrix, so the inverse operation is very simple.

[0104] But there is a problem here. Fig.14 The boundary layer grid is shown, and its distance perpendicular to the wall is much smaller than the distance along the streamline direction, that is,

[0105] Therefore, the gradient result based on the least squares algorithm is mainly composed of the values ​​along the streamline direction (φ 4 ,φ 5 ). However, boundary layer theory tells us that the rate of change of variables in the direction perpendicular to the wall is much greater than the rate of change of variables in the streamline direction. Therefore, formula (29) will produce a large error when calculating the boundary layer grid gradient. Therefore, the weighted coefficient matrix ω related to the distance is introduced into formula (29), and the gradient equation after introducing ω is shown in formula (30), where ω is an M*M matrix, as shown in formula (31). The least squares approximate solution after introducing ω is shown in formula (32). This method does not use the center value of the grid surface to calculate the gradient, so there is no problem of skewness. This method also does not need to calculate the node position value, and the calculation amount is concentrated on the matrix related to the distance (d T ω T ωd) -1 d T ω T ω on. Among them d T ω T ωd is a 3*3 matrix, and its inverse matrix is ​​easy to solve. And most numerical calculations use static grids, which means that (d T ω T ωd) -1 d T ω T The ω matrix only needs to be calculated once, so the computational efficiency of this method is relatively high. In summary, the variables of the flow field control equations use the second-order upwind interpolation method, and the gradient uses the grid-based least squares interpolation method.

[0106]

[0107]

[0108] High-speed flow often involves turbulence, which is an unsteady, random motion phenomenon observed at medium / high Reynolds numbers, and its essence can be described by the Navier-Stokes (NS) equations. However, it is extremely difficult to solve turbulence problems by direct numerical simulation (DNS) in actual operation, so it is usually necessary to perform some form of averaging on the NS equations to eliminate or weaken the influence of the turbulent component. The Reynolds-averaged Navier-Stokes (RANS) equations are a method currently widely used in turbulence simulation, which simplifies the problem by averaging the turbulent fluctuation time terms in the NS equations. In the numerical simulation of the present invention, two relatively accurate and widely used RANS turbulence models, SSTk-omega and Spalart-Allmaras, were selected. For the SSTk-omega turbulence model, the intensity and viscosity ratio are used as the turbulence specification method on the input boundary, where the turbulence intensity is set to 1% and the viscosity ratio is set to 1. For the Spalart-Allmaras turbulence model, the turbulence specification method for the input boundary selects the turbulence viscosity ratio, which takes a value of 1. The hypersonic external flow field needs to consider the compression effect of the fluid, so the viscosity μ is not a constant and Sutherland Law is used. Because the experimental Mach number reaches 10, it should be considered whether there is a real gas effect and whether the atmospheric density ρ can use the ideal gas model. The ideal gas approximation is applicable to the low-pressure gas area of ​​the PT diagram and PV diagram. If the fluid conditions meet P / P c <<1, where P c =3.77MPa, which is the critical pressure of air fluid, so the ideal gas model can be used. The initial pressure of the far-field free flow is 576Pa, and the maximum pressure value of the numerical calculation is about 73728Pa, that is, P / P c The maximum value of is about 0.019, which satisfies the above inequality. Therefore, the ideal gas model can be used for the hypersonic external flow field.

[0109] Based on 3 far-field distances (12L / 16L / 20L), 2 RANS turbulence models (SSTk-omega / Spalart-Allmaras) and 2 flux types (ROE / AUSM), the following numerical simulation configurations are obtained: SST+12L (ROE / AUSM), SA+12L (ROE / AUSM), SST+16L (ROE / AUSM), SA+16L (ROE / AUSM), SST+20L (ROE / AUSM) and SA+20L (ROE / AUSM). Based on the grid strategy and numerical method parameters, 12 sets of numerical simulations were performed for each NACA0012 airfoil. In this way, a total of 72 sets of numerical simulations were performed. The P / P of the airfoil surface was calculated. t 、T / T t and U / U tThe numerical results at the sampling position and the error rate with the wind tunnel test data. The error rates of the six NACA0012 airfoils under different simulation configurations are shown in Table 2, and the optimal solution of the numerical calculation and its corresponding simulation parameters can be determined.

[0110] Table 2 Different parameter configurations and their numerical error rates

[0111]

[0112]

[0113] The Reynolds number of the grid element near the shock surface (R cell ) is a key parameter that affects the numerical error rate. Taking a blunt cylinder as the characteristic object, the prior art points out that R cell The value of should be no less than 8. In addition, another conclusion in the prior art shows that the value of the aspect ratio of the grid near the shock wave surface is also an important factor affecting the numerical accuracy. Based on the above existing research conclusions, the present invention analyzes the influence of the Reynolds number and the aspect ratio of the grid unit near the shock wave surface of the NACA0012 airfoil on the numerical accuracy under hypersonic conditions. cell The calculation of is shown in formula (24), where ρ, U t , μ are the density, velocity and viscosity of the far-field free flow, respectively, y H is the height of the first layer grid unit, so R cell With y + and H First, under the premise that the total number of grids remains unchanged, based on the optimal simulation parameters obtained from the numerical results in Section 2, three kinds of Reynolds numbers of grids near the shock surface are selected to perform numerical simulation. cell ,y + and H The values ​​are (16, 0.3, 1.4e-5), (8, 0.15, 7e-6) and (4, 0.08, 3.5e-6). Tests show that the above three values ​​can ensure the airfoil y + The maximum value of does not exceed 1. As shown in Table 3, the corresponding average error rates under the three grid Reynolds numbers are (2.54%, 1.86%, 1.74%), (2.62%, 2.09%, 2.01%) and (6.38%, 4.77%, 2.47%), respectively. + =0.3, R cell =16, the optimal error rate is obtained. Then, the influence of changing the aspect ratio of the grid near the shock wave surface on the numerical accuracy is numerically studied under the following conditions: (1) the total number of grids remains unchanged; (2) the Reynolds number of the grid near the shock wave surface remains unchanged; (3) only the aspect ratio of the airfoil wall grid in a small range near the shock wave surface is changed. The grid aspect ratios are 760, 380, 190 and 95, respectively, and the corresponding numerical error rates are as follows: Fig.15 As shown in the figure. When the aspect ratio of the near shock wave surface grid is small, the corresponding total average error is large. As the aspect ratio increases, the total average error rate decreases. When the aspect ratio is 380, the optimal simulation result is obtained. Compared with the other three aspect ratios, the numerical accuracy is improved by 63.97%, 46.75% and 65.37%, respectively. As the aspect ratio increases further, the error rate also increases.

[0114]

[0115] In summary, unlike the existing studies that use blunt cylinders as feature objects, the recommended value of the unit Reynolds number under the condition of NACA0012 airfoil as feature object should be no less than 16. Reducing this value will reduce the numerical accuracy. The grid aspect ratio also has a similar situation. The smaller the value, the lower the numerical accuracy. The recommended value is 380.

[0116] Table 3 Average error rate under three grid unit Reynolds numbers

[0117] Unit Reynolds number <![CDATA[P / P t ]]> <![CDATA[T / T t ]]> <![CDATA[U / U t ]]> Total average error rate <![CDATA[16(y + 0.3)]]> 2.54％ 1.86％ 1.74％ 2.05％ <![CDATA[8(y + 0.15)]]> 2.62％ 2.09％ 2.01％ 2.24％ <![CDATA[4(y + 0.08)]]> 6.38％ 4.77％ 2.47％ 4.54％

[0118] Fig.16The optimal average error rate of six NACA0012 airfoils and their numerical prediction accuracy variation are shown. Regarding the numerical performance of the two trailing edge airfoils, it is found that except for the sharp trailing edge airfoil designed based on the definition formula of 264 data points, the blunt trailing edge airfoil designed by Airfoiltools is inferior to the other sharp trailing edge airfoils in prediction accuracy, and its accuracy decreases by 11.41%, 2.14%, 31.73% and 9.88% respectively. Comparing the three modeling methods of Airfoiltools, NACA4 and definition formula, it can be observed that: (1) the optimal simulation error rate of the Airfoiltools method is 2.70%, which is relatively high; (2) the NACA4 method performs better, with an optimal simulation error rate of 2.42%; (3) the definition formula method can achieve the minimum numerical error rate of 2.05% under certain circumstances (based on 200 data points). However, in other cases, the numerical performance of the airfoil designed by NACA4 is still better than that of the definition formula. Further exploration of the effect of the number and source of data points on simulation accuracy revealed the following findings: When using the data provided by Airfoil tools, the increase from 132 data points to 264 actually resulted in a 4.55% drop in accuracy. Similarly, when using the data provided by NACA4, the increase in data points from 200 to 400 also resulted in a 19.71% drop in numerical accuracy. Further analysis of the definition formulas for four different numbers of data points showed that the optimal numerical error rates were 2.64%, 2.76%, 2.05%, and 2.454%, respectively. The definition formulas using 200 and 400 data points performed better than the definition formulas using 132 and 264 data points. Based on the above analysis, the following conclusions can be drawn: First, the selection of the trailing edge shape of the NACA0012 airfoil has a great influence on the simulation accuracy. Improper selection of the trailing edge shape may lead to a significant increase in the simulation error, with the maximum amplitude reaching 50.67%. It is recommended to use a sharp trailing edge airfoil, and the source of the data points used in modeling determines which sharp trailing edge modeling method is used. Secondly, the performance of the modeling data points provided by NACA4 is better than that provided by Airfoil tools, and there is no simple positive correlation between the number of data points and the simulation accuracy, that is, more data points do not necessarily mean improved accuracy. Finally, in response to the simulation requirements of the hypersonic external flow field, it is recommended to use the 200 data points provided by NACA4 to establish the sharp trailing edge NACA0012 airfoil.

[0119] Fig.17The comparison of the far-field distance, turbulence model and flux type in terms of the optimal numerical error rate is shown. Regarding the influence of the far-field distance, it can be observed that as the far-field distance increases from 12L to 16L and even further to 20L, the accuracy of the numerical calculation first significantly increases by 15.42%, and then declines by 19.88%. For the turbulence model, the numerical performance of SST k-omega under the far-field distance conditions of 12L and 20L is better than that of the SA model, and the calculation accuracy is improved by 14.3% and 17.54% respectively. However, under the far-field distance of 16L, the minimum numerical error rate can be achieved based on the SA model. For the flux type, compared with AUSM, ROE has better numerical performance, with the maximum accuracy improvement of 37.3% and the minimum accuracy of 22.68%. In summary, first, there is no positive correlation between the far-field distance and the numerical accuracy, and it is recommended to use the 16L far-field distance to design the flow field outside the calculation domain. Secondly, the choice of turbulence model depends on the far-field distance value. According to the optimal far-field distance of 16L, the SA model should be given priority. Finally, the flux should be of ROE type.

[0120] Table 4 lists in detail the optimal error rate and accuracy change of key parameters. The configuration of SST+12L combined with ROE flux is suitable for blunt / sharp trailing edge airfoils designed using Airfoil tools and NACA4. For the sharp trailing edge airfoil designed according to the defined formula, when using the 200 data points provided by NACA4, the best numerical error rate can be achieved by combining the SA+16L+ROE configuration. It should be emphasized that if the numerical parameters are not selected properly, the simulation accuracy may be significantly reduced. The flux type will bring the largest error rate reduction, followed by the sharp trailing edge airfoil modeling method, and then the shape of the trailing edge airfoil. The reduction caused by the turbulence model and the far field distance is close and relatively small.

[0121] Table 4 Optimal error rate and accuracy variation of key parameters

[0122]

[0123] The present invention provides a simulation parameter optimization method for numerical simulation of hypersonic external flow field, which achieves the following beneficial effects by improving key links such as airfoil modeling, calculation domain setting, grid division, numerical calculation method, and turbulence model optimization:

[0124] 1. Improve the accuracy of airfoil modeling and reduce geometric errors. Use the defined formula modeling method to avoid direct use of database interpolation modeling, thereby improving the smoothness and accuracy of the airfoil geometry curve. Select NACA4 data points (200), optimize the distribution of airfoil curve data points, improve the geometric accuracy of simulation calculations, and reduce numerical errors.

[0125] 2. Optimize the far field distance of the calculation domain to improve the calculation stability. By comparing the effects of 12L, 16L and 20L far field distances on the calculation accuracy, it is found that the 16L far field distance can effectively reduce the boundary reflection error and improve the calculation stability. The far field boundary condition uses Pressure far field to ensure that the flow field boundary conditions are consistent with the actual flow and improve the accuracy of the numerical simulation.

[0126] 3. Improve the grid division strategy to improve the turbulence analysis capability. Set the height of the first layer of grid cells, optimize the wall shear stress calculation, and improve the turbulence analysis accuracy. Use the expansion layer grid strategy to optimize the boundary layer analysis capability. Set the Reynolds number of the grid cells near the shock wave surface to ensure the accuracy of shock wave structure analysis and reduce shock wave diffusion errors. Control the grid aspect ratio, optimize the grid quality, and improve the stability of numerical calculations.

[0127] 4. Improve the numerical calculation method and optimize the shock wave capture accuracy. The density-based solver is suitable for hypersonic compressible flow calculations to improve the solution efficiency. The ROE flux splitting method is used, which can reduce the numerical dissipation error and improve the shock wave capture accuracy compared with the traditional AUSM method. The second-order upwind scheme is used for flow field variable interpolation calculation, combined with the gradient calculation method based on the least squares method, to improve the stability of the numerical solution and reduce the calculation error.

[0128] 5. Optimize the turbulence model to improve calculation accuracy and convergence. Use the Spalart-Allmaras turbulence model to optimize the turbulence simulation capability of the hypersonic flow field and improve calculation accuracy. Set the turbulent viscosity ratio at the input boundary to ensure calculation stability, reduce turbulence errors, and improve simulation convergence.

[0129] 6. Error control optimization to improve simulation reliability. By comparing the simulation results with wind tunnel test data, under the optimized parameter settings, the error rates of pressure, temperature and flow rate are controlled within 2.05%, which reduces the numerical error by more than 50% compared with traditional methods.

[0130] This optimization method can be widely used in the fields of hypersonic aircraft aerodynamic calculation, external flow field optimization design, aerodynamic thermal simulation analysis, etc., to improve the reliability and calculation accuracy of hypersonic aircraft design.

[0131] The present invention achieves improved numerical calculation accuracy, enhanced turbulence analysis capability, and improved calculation stability through system optimization of airfoil modeling, calculation domain setting, grid division, numerical calculation method, and error control. It can ultimately be used for external flow field simulation and optimization design of hypersonic aircraft, and has the beneficial effects of high calculation accuracy, small error, and good numerical stability.

[0132] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention rather than to limit it. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that the technical solutions of the present invention may be modified or replaced by equivalents without departing from the spirit and scope of the technical solutions of the present invention, which should all be included in the scope of the claims of the present invention.

Claims

1. A simulation parameter optimization method for numerical simulation of hypersonic external flow field, characterized in that: include, Select the NACA0012 airfoil with a sharp trailing edge, establish the airfoil mathematical model, and generate the airfoil geometric model; The far field distance of the computational domain is set to 16 times the chord length of the airfoil, the input boundary adopts the Pressure far field condition, the output boundary adopts the Pressure far field condition, and the airfoil wall boundary adopts the No-slip isothermal wall condition; ICEM CFD is used for meshing, and the height of the first layer of mesh cells is set so that y + ≤1; Adopt the expansion layer grid strategy and set the expansion layer and grid growth rate; Set the Reynolds number of the grid unit near the shock wave surface and determine the grid division scheme; The density-based solver is used for hypersonic flow calculations; The ROE flux splitting method is used for numerical calculation; The second-order upwind scheme is used for flow field variable interpolation calculation, and the least square method is used for gradient calculation; The Spalart-Allmaras turbulence model is used for turbulence calculation, and the turbulent viscosity ratio is set at the input boundary; The error of simulation results is evaluated through wind tunnel test data.

2. The method for optimizing simulation parameters for numerical simulation of hypersonic external flow field according to claim 1, characterized in that: The establishment of the airfoil mathematical model refers to the use of the NACA4 definition formula to calculate the airfoil curve, which is expressed as y=±0.5947[0.2983x 1 / 2 -0.1271x-0.3579x 2 +0.292x 3 -0.1052x 4 ] Where y is the airfoil thickness distribution and x is the position of the airfoil surface along the chord length.

3. The method for optimizing simulation parameters for numerical simulation of hypersonic external flow field according to claim 1 or 2, characterized in that: Connect data points, create airfoil surfaces, and generate airfoil geometry models using CAD software or Python / Matlab code.

4. The method for optimizing simulation parameters for numerical simulation of hypersonic external flow field according to claim 3, characterized in that: When setting the height of the first layer of grid cells, the following formula is used for calculation: <h2 style=";text-align:left;direction:ltr">y<h2 style=";text-align:left;direction:ltr"> H <h2 style=";text-align:left;direction:ltr"> <2y<h2 style=";text-align:left;direction:ltr"> p In the formula, y H is the height of the first layer grid unit, y + is the dimensionless wall distance parameter, u τ is the friction velocity, ρ is the atmospheric density, and μ is the free stream kinematic viscosity.

5. The method for optimizing simulation parameters for numerical simulation of hypersonic external flow field according to any one of claims 1, 2 and 4, characterized in that: The expansion layer is set to 30 and the mesh growth rate is set to 1.

05.

6. The method for optimizing simulation parameters for numerical simulation of hypersonic external flow field according to claim 5, characterized in that: For the sharp trailing edge airfoil, the grid unit Reynolds numbers are set to 0.469, 0.505 and 0.493 respectively; for the blunt trailing edge airfoil, the minimum values ​​of grid orthogonality are 0.454, 0.495 and 0.472 respectively.

7. The method for optimizing simulation parameters for numerical simulation of hypersonic external flow field according to claim 6, characterized in that: The grid division scheme is to divide the grid into 610,000 grids.

Images