Aerodynamic thermal environment rapid prediction method based on computational fluid mechanics

Abstract

The invention belongs to the technical field of aerospace engineering, and particularly relates to a computational fluid mechanics-based aerodynamic thermal environment rapid prediction method, which comprises the following steps of: generating a surface grid according to an aerodynamic configuration; generating a pneumatic grid based on the surface grid, and then carrying out non-viscous flow field calculation based on a CFD method to obtain a flow field variable; searching all surface units according to the characteristic that the stationary point speed is 0, and determining the position of a stationary point; coordinates (x, y, z) are used as independent variables, and an epsilon-curve is obtained through calculation around stationary points; using coordinates (x, y, z) as independent variables, adopting a Gaussian kernel function to realize surface streamline tracking until the surface streamline is intersected with the epsilon-curve, and then calculating along the streamline to obtain a shape factor; and based on an axial symmetry comparison idea, predicting the aerodynamic thermal environment along the surface streamline by using an aerodynamic thermal environment calculation formula. In the method, due to the fact that the hybrid grid and the polyhedral grid are higher in calculation efficiency, the calculation cost can be remarkably reduced when aerodynamic thermal environment prediction is carried out based on the method.

Description

Technical Field

[0001] The present invention belongs to the technical field of aerospace engineering, and in particular relates to a method for rapidly predicting an aerodynamic thermal environment based on computational fluid dynamics. Background Art

[0002] The paper "Approximate method for computing convective heating on hypersonic vehicles using unstructured grids," Journal of Spacecraft and Rockets, 51 (2014), pp. 1288-1305, uses computational fluid dynamics (CFD) methods based on axisymmetric analogy to obtain inviscid flow field variables. By calculating shape factors along surface streamlines for triangular surface mesh elements, the proposed method effectively predicts the aerodynamic thermal environment of the vehicle surface. However, the proposed method not only requires constant adjustment of independent variables for different locations, which is not conducive to algorithm implementation, but also only works with triangular surface mesh elements and fails to account for the effects of sideslip angle. However, for practical engineering problems, the following three issues must be considered: 1. The surface mesh must use different mesh element types, i.e., a hybrid mesh, for different locations to achieve satisfactory solution accuracy; 2. Using polyhedral meshes to obtain flow field variables for CFD calculations is highly efficient; and 3. The sideslip angle significantly affects the aerodynamic thermal environment. In summary, how to quickly and accurately predict the aerodynamic thermal environment based on surface hybrid grids and polyhedral grid units, and taking into account the influence of sideslip angle, thereby improving the efficiency and accuracy of the aircraft design stage, is a technical problem to be solved in this field. Summary of the Invention

[0003] (1) Technical issues to be solved

[0004] The technical problem to be solved by the present invention is: in order to solve the prominent problem that the method in the literature is only applicable to triangular surface grid units and cannot consider the influence of sideslip angle, and at the same time break through the limitation of the literature method that requires continuous transformation of independent variables, which is not conducive to algorithm implementation, the present invention aims to propose a new aerodynamic thermal environment rapid prediction method based on computational fluid dynamics, which can not only overcome the limitation of the literature method in algorithm implementation caused by the need to continuously transform independent variables, but more importantly, it can be applicable to surface mixed grids and polyhedral grid units, and can be applied to working conditions with sideslip angle, while obtaining satisfactory prediction accuracy, reducing the computational cost required for aerodynamic thermal environment prediction, and improving the engineering applicability of the method.

[0005] (2) Technical solution

[0006] In order to solve the above technical problems, the present invention provides a method for rapid prediction of aerodynamic thermal environment based on computational fluid dynamics. Figure 1 The following is a flow chart of the implementation of the method of the present invention: the method comprises the following steps:

[0007] Step 1: Generate surface mesh based on aerodynamic shape;

[0008] Step 2: Generate an aerodynamic mesh based on the surface mesh in step 1, and then perform inviscid flow field calculations based on the CFD method to obtain flow field variables;

[0009] Step 3: Based on the characteristic that the stagnation point velocity is 0, search all surface units to determine the location of the stagnation point;

[0010] Step 4: Using the coordinates (x, y, z) as independent variables, calculate the ε-curve around the stationary point;

[0011] Step 5: Also using the coordinates (x, y, z) as independent variables, based on the radial basis function (RBF) interpolation method, the Gaussian kernel function is used to trace the surface streamlines until they intersect with the ε-curve, and then the shape factor is calculated along the streamlines;

[0012] Step 6: Based on the axisymmetric analogy, the aerodynamic thermal environment calculation formula is used to predict the aerodynamic thermal environment along the surface streamlines.

[0013] Wherein, in step 1, the surface grid unit is any type such as triangle, quadrilateral or polygon.

[0014] Wherein, in said step 3, according to the characteristic that the stagnation point velocity is 0, all surface units are searched to determine the position of the stagnation point;

[0015] The aerodynamic shape surface is expressed by formula (1)

[0016] F(x,y,z)=0 (1)

[0017] In formula (1), (x, y, z) are coordinate variables in the Cartesian coordinate system; (u, v, w) are defined as the velocity variables of the air in the (x, y, z) directions in the Cartesian coordinate system, and it is known that the velocity at the stagnation point satisfies (u, v, w) = 0; for hypersonic air flow, (v, w) is small compared to u, so in each unit, let (v, w) = 0, and use (y, z) as the independent variable to search the grid cells within a reasonable range of the head to determine the stagnation point position.

[0018] Wherein, in said step 4, the coordinates (x, y, z) are used as independent variables to calculate the ε-curve around the stationary point;

[0019] like Figure 2As shown, a streamline coordinate system (ξ,η,τ) is defined on the aerodynamic surface; where ξ is the streamline direction, η is tangent to the surface and perpendicular to ξ, τ is the direction of the surface normal, and u is the velocity of the air in the x direction in the Cartesian coordinate system in step 3. At the same time, are the unit vectors in the (ξ, η, τ) directions in the streamline coordinate system, so:

[0020] In the Cartesian coordinate system, the position vector of any point on the surface is The coordinates (x, y, z) of the point are expressed as follows

[0021]

[0022] In formula (2), are the unit vectors in the (x, y, z) directions in the Cartesian coordinate system;

[0023] In the streamline coordinate system, the position vector of any point on the surface is Expressed as

[0024]

[0025] Unit vector Calculated using formula (4)

[0026]

[0027] Air velocity vector In the Cartesian coordinate system, it is expressed as formula (5)

[0028]

[0029] Air velocity vector In the streamline coordinate system, it is expressed as formula (6)

[0030]

[0031] In formula (6), V is the absolute velocity of air, and its value is

[0032] Therefore, combining formula (5) and formula (6) we can get Expression

[0033]

[0034] Definition (F x ,F y ,F z ) is the derivative of the aerodynamic surface function F(x,y,z) in the (x,y,z) direction, For (F x,F y ,F z ) but Expressed as follows

[0035]

[0036] Combining formulas (4), (7) and (8), we get The expression is as follows

[0037]

[0038] At this point, the ε-curve around the stationary point is calculated using the following formula

[0039]

[0040] In formula (10), s η is the lower edge of the streamline coordinate system The ε-curve can be obtained by integrating formula (10) around the stationary point.

[0041] Wherein, in step 4, the formula (10) is applicable to both the conditions with attack angle and the conditions with sideslip angle.

[0042] In step 5, the coordinates (x, y, z) are also used as independent variables. Based on the radial basis function (RBF) interpolation method, a Gaussian kernel function is used to track the surface streamlines until they intersect with the ε-curve, and then the shape factor is calculated along the streamlines.

[0043] According to the definition of streamline, in the Cartesian coordinate system, the streamline can be traced by integrating the time t according to formula (11);

[0044]

[0045] According to formula (2), the unit vector Expressed as

[0046]

[0047] In formula (12), They are The partial derivative with respect to the streamline coordinate η is, Corresponding to the partial derivative The amplitude of

[0048] Combining formula (9) and formula (12) we can get

[0049]

[0050] definition is the shape factor variable along the streamline; Formula (13) is rewritten as

[0051]

[0052] And there are

[0053]

[0054] Using the chain derivation rule of composite functions, we can obtain the following relationship:

[0055]

[0056] In formula (16), are the partial derivatives of the air velocity variables (u, v, w) with respect to the coordinates (x, y, z); by giving the initial value h0, the initial value can be calculated according to formula (14) Then, using formula (16) to integrate along the streamline, and combining with formula (15), the shape factor variable h along the streamline can be calculated;

[0057] In this process, the 45th-order adaptive Longo Kutta method is used to integrate formulas (11) and (16). It is not difficult to find that in the integration process, it is necessary to calculate the integral of each integral time step. and The value of , and the corresponding value of each grid unit also need to be calculated In addition, it is also necessary to interpolate the air velocity variables (u, v, w) and other flow field variables at each integration point, including pressure p and density ρ;

[0058] In order to be applicable to triangular meshes, hybrid meshes, and polyhedral meshes, this step uses the RBF interpolation method to implement variable interpolation, and derives the partial derivative expressions of all variables; in order to avoid the zero division phenomenon, the Gaussian kernel function is selected as the interpolation kernel function; the partial derivative of the air velocity variable u in the x direction is given below expression;

[0059] In any unit e, u is expressed as follows

[0060]

[0061] In formula (17), is the RBF interpolation kernel function, d = ‖xx m ‖2 is the Euclidean distance between two points; n is the number of interpolation nodes, c is the interpolation coefficient, and the subscript m represents the mth point;

[0062] By formula (17) It is expressed as follows

[0063]

[0064] Define σ as an adjustable parameter and set the Gaussian kernel function Substituting the expression into formula (18), we get

[0065]

[0066] Comparing formula (18) and formula (19), we find that by using the Gaussian kernel function, we can avoid the zero division phenomenon when d→0; similarly, the other partial derivatives The same method is used to calculate it; by integrating formula (11) and formula (16) along time and combining formula (19) and other derived expressions, the aerodynamic thermal environment can be predicted.

[0067] Wherein, in said step 6, based on the axisymmetric analogy idea, the aerodynamic thermal environment is predicted along the surface streamline using the aerodynamic thermal environment calculation formula;

[0068] For laminar flow, the heat flux Calculate using the following formula

[0069]

[0070] Where ρ, u, and μ are the density, velocity, and viscosity coefficient of air, respectively. The superscript * indicates that it is calculated based on the Eckert reference enthalpy method. The subscript e indicates the boundary layer outer edge parameter. H w ,H e and H r are the wall enthalpy, boundary layer outer edge enthalpy and recovery enthalpy of air respectively; is the momentum thickness Reynolds number, and the momentum thickness θ along the streamline coordinate system is calculated using the following formula

[0071]

[0072] In formula (21), s is the lower edge of the streamline coordinate system Distance of direction;

[0073] The Eckert reference enthalpy relationship is shown below

[0074] H * =0.19H r +0.23H e +0.58H w (twenty two)

[0075] In formula (22), Pr w is the Prandtl number corresponding to the aerodynamic surface, and Pr w=0.71.

[0076] (3) Beneficial effects

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

[0078] (1) The method proposed in the present invention can realize aerodynamic thermal environment prediction for triangular grid cells, triangular / quadrilateral hybrid grid cells, and polyhedral grid cells, and has stronger engineering applicability;

[0079] (2) The method of the present invention can predict the aerodynamic thermal environment for the side slip angle condition;

[0080] (3) Since hybrid grids and polyhedron grids have higher computational efficiency, the aerodynamic thermal environment prediction based on the method of the present invention can significantly reduce the computational cost. BRIEF DESCRIPTION OF THE DRAWINGS

[0081] Figure 1 It is an implementation flow chart of the method of the present invention.

[0082] Figure 2 It is a schematic diagram of the streamline coordinate system.

[0083] Figure 3 This is a diagram of the aerodynamic thermal environment results obtained by applying the method of the present invention to a hypersonic blunt double-cone example.

[0084] Figure 4 This is a comparison diagram of the aerodynamic thermal environment results obtained by applying the method of the present invention to a hypersonic blunt double-cone example.

[0085] Figure 5 This is a streamline tracking result diagram obtained by applying the method of the present invention to a hypersonic double ellipsoid example.

[0086] Figure 6 This is a comparison diagram of the aerodynamic thermal environment results obtained by applying the method of the present invention to the hypersonic double ellipsoid calculation example.

[0087] Figure 7 This is a diagram of the aerodynamic thermal environment results obtained by applying the method of the present invention to a hypersonic double ellipsoid example at different sideslip angles.

[0088] Figure 8 This is a comparison chart of the mesh size generated based on different surface meshes and the time required for CFD inviscid flow field calculation for the hypersonic double ellipsoid example. DETAILED DESCRIPTION

[0089] In order to make the purpose, content, and advantages of the present invention more clear, the specific implementation methods of the present invention are further described in detail below with reference to the accompanying drawings and examples.

[0090] In order to solve the above technical problems, the present invention provides a method for rapid prediction of aerodynamic thermal environment based on computational fluid dynamics. Figure 1 The following is a flow chart of the implementation of the method of the present invention: the method comprises the following steps:

[0091] Step 1: Generate surface mesh based on aerodynamic shape;

[0092] Step 2: Generate an aerodynamic mesh based on the surface mesh in step 1, and then perform inviscid flow field calculations based on the CFD method to obtain flow field variables;

[0093] Step 3: Based on the characteristic that the stagnation point velocity is 0, search all surface units to determine the location of the stagnation point;

[0094] Step 4: Using the coordinates (x, y, z) as independent variables, calculate the ε-curve around the stationary point;

[0095] Step 5: Also using the coordinates (x, y, z) as independent variables, based on the radial basis function (RBF) interpolation method, the Gaussian kernel function is used to trace the surface streamlines until they intersect with the ε-curve, and then the shape factor is calculated along the streamlines;

[0096] Step 6: Based on the axisymmetric analogy, the aerodynamic thermal environment calculation formula is used to predict the aerodynamic thermal environment along the surface streamlines.

[0097] Wherein, in step 1, the surface grid unit is any type such as triangle, quadrilateral or polygon.

[0098] Wherein, in said step 3, according to the characteristic that the stagnation point velocity is 0, all surface units are searched to determine the position of the stagnation point;

[0099] The aerodynamic shape surface is expressed by formula (1)

[0100] F(x,y,z)=0 (1)

[0101] In formula (1), (x, y, z) are coordinate variables in the Cartesian coordinate system; (u, v, w) are defined as the velocity variables of the air in the (x, y, z) directions in the Cartesian coordinate system, and it is known that the velocity at the stagnation point satisfies (u, v, w) = 0; for hypersonic air flow, (v, w) is small compared to u, so in each unit, let (v, w) = 0, and use (y, z) as the independent variable to search the grid cells within a reasonable range of the head to determine the stagnation point position.

[0102] Wherein, in said step 4, the coordinates (x, y, z) are used as independent variables to calculate the ε-curve around the stationary point;

[0103] like Figure 2As shown, a streamline coordinate system (ξ,η,τ) is defined on the aerodynamic surface; where ξ is the streamline direction, η is tangent to the surface and perpendicular to ξ, τ is the direction of the surface normal, and u is the velocity of the air in the x direction in the Cartesian coordinate system in step 3. At the same time, are the unit vectors in the (ξ, η, τ) directions in the streamline coordinate system, so:

[0104] In the Cartesian coordinate system, the position vector of any point on the surface is The coordinates (x, y, z) of the point are expressed as follows

[0105]

[0106] In formula (2), are the unit vectors in the (x, y, z) directions in the Cartesian coordinate system;

[0107] In the streamline coordinate system, the position vector of any point on the surface is Expressed as

[0108]

[0109] Unit vector Calculated using formula (4)

[0110]

[0111] Air velocity vector In the Cartesian coordinate system, it is expressed as formula (5)

[0112]

[0113] Air velocity vector In the streamline coordinate system, it is expressed as formula (6)

[0114]

[0115] In formula (6), V is the absolute velocity of air, and its value is

[0116] Therefore, combining formula (5) and formula (6) we can get Expression

[0117]

[0118] Definition (F x ,F y ,F z ) is the derivative of the aerodynamic surface function F(x,y,z) in the (x,y,z) direction, For (F x,F y ,F z ) but Expressed as follows

[0119]

[0120] Combining formulas (4), (7) and (8), we get The expression is as follows

[0121]

[0122] At this point, the ε-curve around the stationary point is calculated using the following formula

[0123]

[0124] In formula (10), s η is the lower edge of the streamline coordinate system The ε-curve can be obtained by integrating formula (10) around the stationary point.

[0125] Wherein, in step 4, the formula (10) is applicable to both the conditions with attack angle and the conditions with sideslip angle.

[0126] In step 5, the coordinates (x, y, z) are also used as independent variables. Based on the radial basis function (RBF) interpolation method, a Gaussian kernel function is used to track the surface streamlines until they intersect with the ε-curve, and then the shape factor is calculated along the streamlines.

[0127] According to the definition of streamline, in the Cartesian coordinate system, the streamline can be traced by integrating the time t according to formula (11);

[0128]

[0129] According to formula (2), the unit vector Expressed as

[0130]

[0131] In formula (12), They are The partial derivative with respect to the streamline coordinate η is, Corresponding to the partial derivative The amplitude of

[0132] Combining formula (9) and formula (12) we can get

[0133]

[0134] definition is the shape factor variable along the streamline; Formula (13) is rewritten as

[0135]

[0136] And there are

[0137]

[0138] Using the chain derivation rule of composite functions, we can obtain the following relationship:

[0139]

[0140] In formula (16), are the partial derivatives of the air velocity variables (u, v, w) with respect to the coordinates (x, y, z); by giving the initial value h0, the initial value can be calculated according to formula (14) Then, using formula (16) to integrate along the streamline, and combining with formula (15), the shape factor variable h along the streamline can be calculated;

[0141] In this process, the 45th-order adaptive Longo Kutta method is used to integrate formulas (11) and (16). It is not difficult to find that in the integration process, it is necessary to calculate the integral of each integral time step. and The value of , and the corresponding value of each grid unit also need to be calculated In addition, it is also necessary to interpolate the air velocity variables (u, v, w) and other flow field variables at each integration point, including pressure p and density ρ;

[0142] In order to be applicable to triangular meshes, hybrid meshes, and polyhedral meshes, this step uses the RBF interpolation method to implement variable interpolation, and derives the partial derivative expressions of all variables; in order to avoid the zero division phenomenon, the Gaussian kernel function is selected as the interpolation kernel function; the partial derivative of the air velocity variable u in the x direction is given below expression;

[0143] In any unit e, u is expressed as follows

[0144]

[0145] In formula (17), is the RBF interpolation kernel function, d = ‖xx m ‖2 is the Euclidean distance between two points; n is the number of interpolation nodes, c is the interpolation coefficient, and the subscript m represents the mth point;

[0146] By formula (17) It is expressed as follows

[0147]

[0148] Define σ as an adjustable parameter and set the Gaussian kernel function Substituting the expression into formula (18), we get

[0149]

[0150] Comparing formula (18) and formula (19), we find that by using the Gaussian kernel function, we can avoid the zero division phenomenon when d→0; similarly, the other partial derivatives The same method is used to calculate it; by integrating formula (11) and formula (16) along time and combining formula (19) and other derived expressions, the aerodynamic thermal environment can be predicted.

[0151] Wherein, in said step 6, based on the axisymmetric analogy idea, the aerodynamic thermal environment is predicted along the surface streamline using the aerodynamic thermal environment calculation formula;

[0152] For laminar flow, the heat flux Calculate using the following formula

[0153]

[0154] Where ρ, u, and μ are the density, velocity, and viscosity coefficient of air, respectively. The superscript * indicates that it is calculated based on the Eckert reference enthalpy method. The subscript e indicates the boundary layer outer edge parameter. H w ,H e and H r are the wall enthalpy, boundary layer outer edge enthalpy and recovery enthalpy of air respectively; is the momentum thickness Reynolds number, and the momentum thickness θ along the streamline coordinate system is calculated using the following formula

[0155]

[0156] In formula (21), s is the lower edge of the streamline coordinate system Distance of direction;

[0157] The Eckert reference enthalpy relationship is shown below

[0158] H * =0.19H r +0.23H e +0.58H w (twenty two)

[0159] In formula (22), Pr w is the Prandtl number corresponding to the aerodynamic surface, and Pr w=0.71.

[0160] Example 1

[0161] Based on the present invention's method, simulations were conducted on the blunt double-cone and hypersonic double-ellipsoid examples at Mach 9.86 from NASA TP-2334. The present invention is described in detail with reference to the experimental results. Obviously, the described embodiments are only a portion of the present invention, not all of it. All other embodiments derived by persons of ordinary skill in the art based on the embodiments described herein without inventive effort are also within the scope of protection of the present invention.

[0162] 1. NASA_TP_2334 blunt double cone example

[0163] Figure 3 The aerodynamic thermal environment results obtained by applying the method of the present invention to the hypersonic blunt double cone calculation example are shown in the figure. Qw is used to represent the heat flux density in formula (22). The calculation conditions are Mach number 9.86, a high angle of attack α = 16°, a sideslip angle β = 0°, and laminar flow. A hybrid mesh of triangles and quadrilaterals is used for the surface mesh. As can be seen, the proposed method can produce reasonably distributed and smooth heat flux predictions, demonstrating the successful application of the proposed method to hybrid meshes. Figure 4 The comparison diagram of the predicted aerodynamic thermal environment results on the symmetry plane and the test results. In the figure, q is used to represent the heat flux density in formula (22). Windward represents the direction into the wind, Leeward represents the direction from the leeward, and the prefix Exp indicates that the data point is experimental data. As can be seen from the figure, the thermal environment predicted by the method of the present invention is consistent with the experimental results, indicating that the method of the present invention can also achieve satisfactory calculation accuracy.

[0164] 2. Hypersonic double ellipsoid example

[0165] The calculation conditions for the hypersonic double ellipsoid example are Mach number 8.04, angle of attack α = 0°, sideslip angle β = 0°, and laminar flow. Polyhedral elements are used to mesh the complex surface, which is more computationally efficient. Figure 5 The following figure shows the streamline tracking effect of the method of the present invention in the case of hypersonic double ellipsoid. The red line is obtained by using the commercial software Tecplot, and the green point is the data point obtained by the method of the present invention. It can be seen from the figure that the method of the present invention can achieve a good streamline tracking. Figure 6 As shown in the figure, q is used to represent the heat flux density in formula (22) q_ref represents the reference heat flux used for normalization, polyhedral grid represents a polyhedral surface grid, and tet grid represents a triangular surface grid. It can be seen that the results predicted by the method of the present invention agree well with the experimental results in most areas. The only area where the results agree poorly with the experimental results is the boundary between the two ellipsoids. This is due to the complex shock wave environment in this area, a phenomenon also found in many other literatures. In summary, the method of the present invention can not only be applied to polyhedral grids, but also predict thermal environments that meet engineering precision requirements, and the mesh type has strong applicability.

[0166] In order to show the applicability of the method of the present invention to the side slip angle condition, Figure 7 The comparison of the predicted thermal environment results under the calculation conditions of sideslip angle β = 0° and β = 5° is given. In the figure, Qw is used to represent the heat flux density in formula (22). It can be seen that the method of the present invention successfully predicts the aerodynamic thermal environment of the side slip angle working condition, which greatly expands the scope of engineering application. In addition, the number of grids generated by using triangular grid units and polyhedral grid units in the hypersonic double ellipsoid calculation example is compared, and the calculation time required for CFD calculation of the inviscid flow field is compared. Figure 8 As shown in the figure, it can be seen that: 1. On the basis of achieving the same calculation accuracy, the mesh size generated based on the polyhedral surface mesh is 126,000, compared with 288,000 meshes generated based on the triangular surface mesh, which is a 56% reduction in mesh size; 2. The CFD calculation based on the polyhedral mesh takes 49 seconds per 50 steps, while the CFD calculation based on the tetrahedral mesh takes 66 seconds per 50 steps, which is a 25% reduction in calculation time. Considering that the mesh size is often larger for actual engineering problems, the aerodynamic thermal environment prediction based on the method of the present invention can significantly reduce the calculation cost.

[0167] The above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the technical principles of the present invention. These improvements and modifications should also be regarded as the scope of protection of the present invention.

Claims

1. A method for rapid prediction of aerodynamic thermal environment based on computational fluid dynamics, characterized in that: The method comprises the following steps: Step 1: Generate surface mesh based on aerodynamic shape; Step 2: Generate an aerodynamic mesh based on the surface mesh in step 1, and then perform inviscid flow field calculations based on the CFD method to obtain flow field variables; Step 3: Based on the characteristic that the stagnation point velocity is 0, search all surface units to determine the location of the stagnation point; Step 4: Using the coordinates (x, y, z) as independent variables, calculate the ε-curve around the stationary point; Step 5: Also using the coordinates (x, y, z) as independent variables, based on the radial basis function interpolation method, the Gaussian kernel function is used to trace the surface streamlines until they intersect with the ε-curve, and then the shape factor is calculated along the streamlines; Step 6: Based on the axisymmetric analogy, the aerodynamic thermal environment calculation formula is used to predict the aerodynamic thermal environment along the surface streamlines.

2. The method for rapid prediction of aerodynamic thermal environment based on computational fluid dynamics according to claim 1, characterized in that: In step 1, the surface grid units are triangles, quadrilaterals or polygons.

3. The method for rapid prediction of aerodynamic thermal environment based on computational fluid dynamics according to claim 2, characterized in that: In step 3, based on the characteristic that the stagnation point velocity is 0, all surface units are searched to determine the position of the stagnation point; The aerodynamic shape surface is expressed by formula (1) F(x,y,z)=0 (1) In formula (1), (x, y, z) are coordinate variables in the Cartesian coordinate system; (u, v, w) are defined as the velocity variables of the air in the (x, y, z) directions in the Cartesian coordinate system, and it is known that the velocity at the stagnation point satisfies (u, v, w) = 0; for hypersonic air flow, (v, w) is small compared to u, so in each unit, let (v, w) = 0, and use (y, z) as the independent variable to search the grid cells within a reasonable range of the head to determine the stagnation point position.

4. The method for rapid prediction of aerodynamic thermal environment based on computational fluid dynamics according to claim 3, characterized in that: In step 4, the coordinates (x, y, z) are used as independent variables to calculate the ε-curve around the stationary point; Define the streamline coordinate system (ξ,η,τ) on the aerodynamic surface; where ξ is the streamline direction, η is tangent to the surface and perpendicular to ξ, τ is the direction of the surface normal, and u is the velocity of the air in the x direction in the Cartesian coordinate system in step 3. At the same time, define are the unit vectors in the (ξ, η, τ) directions in the streamline coordinate system, so: In the Cartesian coordinate system, the position vector of any point on the surface is The coordinates (x, y, z) of the point are expressed as follows In formula (2), are the unit vectors in the (x, y, z) directions in the Cartesian coordinate system; In the streamline coordinate system, the position vector of any point on the surface is Expressed as Unit vector Calculated using formula (4) Air velocity vector In the Cartesian coordinate system, it is expressed as formula (5) Air velocity vector In the streamline coordinate system, it is expressed as formula (6) In formula (6), V is the absolute velocity of air, and its value is Therefore, combining formula (5) and formula (6) we can get Expression Definition (F x ,F y ,F z ) is the derivative of the aerodynamic surface function F(x,y,z) in the (x,y,z) direction, For (F x ,F y ,F z ) but Expressed as follows Combining formulas (4), (7) and (8), we get The expression is as follows At this point, the ε-curve around the stationary point is calculated using the following formula In formula (10), s η is the lower edge of the streamline coordinate system The ε-curve can be obtained by integrating formula (10) around the stationary point.

5. The method for rapid prediction of aerodynamic thermal environment based on computational fluid dynamics according to claim 4, characterized in that: In step 4, the formula (10) is applicable to both the conditions with attack angle and the conditions with sideslip angle.

6. The method for rapid prediction of aerodynamic thermal environment based on computational fluid dynamics according to claim 5, characterized in that: In step 5, the coordinates (x, y, z) are also used as independent variables. Based on the radial basis function interpolation method, the Gaussian kernel function is used to track the surface streamlines until they intersect with the ε-curve, and then the shape factor is calculated along the streamlines. According to the definition of streamline, in the Cartesian coordinate system, the streamline can be traced by integrating the time t according to formula (11); According to formula (2), the unit vector Expressed as In formula (12), They are The partial derivatives of x, y, and z with respect to the streamline coordinate η, Corresponding to the partial derivative The amplitude of Combining formula (9) and formula (12) we can get definition is the shape factor variable along the streamline; Formula (13) is rewritten as And there are Using the chain derivation rule of composite functions, we can obtain the following relationship: In formula (16), are the partial derivatives of the air velocity variables (u, v, w) with respect to the coordinates (x, y, z); by giving the initial value h0, the initial value can be calculated according to formula (14) Then, using formula (16) to integrate along the streamline, and combining with formula (15), the shape factor variable h along the streamline can be calculated; In this process, the 45th-order adaptive Longo Kutta method is used to integrate formulas (11) and (16). It is found that in the integration process, it is necessary to calculate the integral of each integral time step. and The value of F corresponding to each grid unit also needs to be calculated x ,F y ,F z , In addition, it is also necessary to interpolate the air velocity variables (u, v, w) and other flow field variables at each integration point, including pressure p and density ρ; In order to be applicable to triangular meshes, hybrid meshes, and polyhedral meshes, this step uses a radial basis function interpolation method to implement variable interpolation, and derives the partial derivative expressions of all variables; in order to avoid the division by zero phenomenon, the Gaussian kernel function is selected as the interpolation kernel function; the partial derivative of the air velocity variable u in the x direction is given below expression; In any unit e, u is expressed as follows In formula (17), is the radial basis function interpolation kernel function, d = ‖xx m ‖2 is the Euclidean distance between two points; n is the number of interpolation nodes, c is the interpolation coefficient, and the subscript m represents the mth point; By formula (17) It is expressed as follows Define σ as an adjustable parameter and set the Gaussian kernel function Substituting the expression into formula (18), we get Comparing formula (18) and formula (19), we find that by using the Gaussian kernel function, the zero division phenomenon when d→0 can be avoided; Similarly, other partial derivatives The same method is used to calculate it; by integrating formula (11) and formula (16) along time and combining formula (19) and other derived expressions, the aerodynamic thermal environment can be predicted.

7. The method for rapid prediction of aerodynamic thermal environment based on computational fluid dynamics according to claim 6, characterized in that: In step 6, based on the axisymmetric analogy concept, the aerodynamic thermal environment calculation formula is used to predict the aerodynamic thermal environment along the surface streamlines; For laminar flow, the heat flux Calculate using the following formula Where ρ, u, and μ are the density, velocity, and viscosity coefficient of air, respectively. The superscript * indicates that it is calculated based on the Eckert reference enthalpy method. The subscript e indicates the boundary layer outer edge parameter. H w ,H e and H r are the wall enthalpy, boundary layer outer edge enthalpy and recovery enthalpy of air respectively; is the momentum thickness Reynolds number, and the momentum thickness θ along the streamline coordinate system is calculated using the following formula In formula (21), s is the lower edge of the streamline coordinate system Distance of direction; The Eckert reference enthalpy relationship is shown below H * =0.19H r +0.23H e +0.58H w (22) In formula (22), Pr w is the Prandtl number corresponding to the aerodynamic shape surface.

8. The method for rapid prediction of aerodynamic thermal environment based on computational fluid dynamics according to claim 7, characterized in that: In step 6, take Pr w =0.

71.

9. The method for rapid prediction of aerodynamic thermal environment based on computational fluid dynamics according to claim 7, characterized in that: The method can realize aerodynamic thermal environment prediction for triangular grid units, triangular / quadrilateral mixed grid units and polyhedron grid units, and has strong engineering applicability.

10. The method for rapid prediction of aerodynamic thermal environment based on computational fluid dynamics according to claim 7, characterized in that: The method can predict the aerodynamic thermal environment for a side slip angle operating condition.

Images