High-speed aircraft thermal protection plate thermodynamic elasticity rapid calculation method based on aerodynamic force / thermal agent model

Abstract

The invention discloses an aerodynamic force / thermal agent model-based rapid calculation method for the thermo-dynamic elasticity of a thermal protection plate of a high-speed aircraft, and the method comprises the steps: building an unsteady aerodynamic force / thermal agent model of a thermal protection wall plate in a high-speed flow based on a Kriging agent model method; a structural dynamics method based on a finite element method is used for establishing an aeroelastic module, a transient heat conduction equation is used for establishing an aerothermal module, and a thermal protection wallboard thermodynamic elasticity analysis frame is built. Finally, the aerodynamic force / thermal agent model is used for carrying out thermal protection wall plate thermodynamic elasticity calculation and analysis. Compared with a traditional high-order method, the calculation efficiency is improved by four orders of magnitude.

Description

Technical Field

[0001] The present invention belongs to the technical field of aircraft, and particularly relates to a rapid calculation method for thermo-aeroelasticity of a thermal protection panel of a high-speed aircraft based on an aerodynamic / thermal surrogate model. Background Art

[0002] Due to the severe aerodynamic heating phenomenon existing in high-speed flight, it will lead to material degradation, thermal stress and thermal bending moment under temperature change, and finally bring about the thermo-aero-elastic multi-physical field coupling problem of the aircraft structure. Especially for the thin-walled structures widely existing in high-speed aircraft such as Figure 1 , the thermo-aeroelastic problem is particularly prominent. In addition, the air flow state and physical properties in the high-speed flow change violently, and the calculation of aerodynamic force / heat becomes the part with the largest computational effort in the research. Therefore, the calculation of aerodynamic force and heat is the key point and difficulty in the research.

[0003] Typical aerodynamic force / heat calculation methods include experimental measurement, engineering algorithms and numerical calculation. However, the engineering algorithms are efficient but lack accuracy, and the numerical calculation has high accuracy but requires a large amount of computing resources.

[0004] Aerodynamic engineering algorithms include piston theory, Van Dyke theory, unified lifting surface theory, etc. With the development of computer technology, using numerical methods to solve the CFD method enables high-fidelity and high-resolution simulation of complex flow fields. Many researchers have carried out aeroelasticity and thermo-aeroelasticity research using aerodynamic engineering algorithms and CFD methods, indicating that there are errors between aerodynamic engineering algorithms and CFD methods.

[0005] Aerodynamic heat engineering algorithms are established based on a large amount of existing data, empirical formulas and combined with theoretical analysis. There are mainly Fay-Riddle formula, Lees formula, etc. Currently, the widely used methods are reference temperature method and Eckert reference enthalpy method. The numerical calculation method of aerodynamic heat directly solves the flow field control equation and has high solution accuracy. However, the numerical format, grid strategy and turbulence model in the calculation will all affect the solution.

[0006] To solve the contradiction between the calculation accuracy and efficiency of aerodynamic force / heat, model reduction methods have been developed to establish aerodynamic force and aerodynamic heat models with both accuracy and computational efficiency. Typical aerodynamic force reduction methods such as Volterra series method, artificial neural network model, modular model, etc. Typical aerodynamic heat reduction methods such as global aerodynamic heat model directly take the aerodynamic heat distribution as the output.

[0007] However, the above aerodynamic force / heat reduction models have problems such as difficult parameter adjustment and too large model dimension. Summary of the Invention

[0008] To overcome the deficiencies of the prior art, the present invention provides a rapid computational method for thermo-aeroelasticity of the thermal protection panel of a high-speed aircraft based on an aerodynamic / thermal surrogate model. Based on the Kriging surrogate model method, an unsteady aerodynamic / thermal surrogate model of the thermal protection wall panel in high-speed flow is established. The aerodynamic elastic module is established using the structural dynamics method based on the finite element method, and the aerothermal module is established using the transient heat conduction equation to build the thermo-aeroelastic analysis framework of the thermal protection wall panel. Finally, the aerodynamic / thermal surrogate model is used to carry out the thermo-aeroelastic calculation and analysis of the thermal protection wall panel. Compared with the traditional high-order method, the computational efficiency of the present invention is increased by four orders of magnitude.

[0009] The technical solution adopted by the present invention to solve its technical problems is as follows:

[0010] Step 1: Set the research object;

[0011] Take the TPS wall panel on the surface of the forebody-inlet of a high-speed aircraft as the research object. The TPS wall panel is composed of a thermal protection layer and a bottom metal plate. Among them, the thermal protection layer includes a radiation layer and a heat insulation layer, which are respectively used to radiate heat to the environment and block the downward conduction of heat; the two ends of the bottom metal plate are simply supported; when the aircraft is in high-speed oncoming flow, an oblique shock wave is generated at the leading edge of the sharp wedge, and the airflow behind the wave flows parallel to the surface of the sharp wedge. The gas parameters at the leading edge of the elastic wall panel are the same as those in the airflow area behind the wave;

[0012] Step 2: Based on the Kriging surrogate model method, establish an unsteady aerodynamic / thermal surrogate model of the thermal protection wall panel in high-speed flow;

[0013] Step 2-1: Determine the state space of the system: Parametrize the input and output variables of the system to obtain the input and output parameters expected by the surrogate model, and determine the parameter boundaries;

[0014] Step 2-2: Sample the parameter state space: Use the Latin hypercube sampling method to select a certain number of state points in the state space as samples and for verification;

[0015] Step 2-3: Response of the sample points: Use the high-order method to calculate the working conditions corresponding to the sample points to generate training data;

[0016] Step 2-4: Generate the surrogate model: Based on the input and output of the sample state points and the response results, establish the surrogate model;

[0017] Step 2-5: Accuracy evaluation: Evaluate the accuracy of the surrogate model at the evaluation points according to the cross-validation idea;

[0018] Step 2-6: Apply the surrogate model for prediction: If the requirements are met, the surrogate model can be used.

[0019] Step 3: Establish a thermo-aeroelastic coupling model;

[0020] The thermo-pneumoelastic coupling analysis method is used to decompose the three-field coupling system into two modules, namely two coupled subsystems: the thermal module and the pneumoelastic module. Among them, the thermal module is a coupled system of aerothermal and structural transient heat conduction. The aerothermal calculation obtains the heat flux and transmits it to the heat conduction module. The heat conduction calculation obtains the transient temperature distribution of the structure and then transmits it to the aerothermal calculation to update the heat flux at the next time step. The pneumoelastic module calculates the unsteady aerodynamic force by the aerodynamic force model and transmits it to the structural dynamics module. The elastic deformation is calculated by the structural dynamics model and then transmitted to the aerodynamic force calculation to update the aerodynamic force at the next time step.

[0021] Preferably, the Kriging model is specifically expressed as follows:

[0022] Step 4-1: The mathematical expression of the Kriging model is written as a linear superposition of the global estimate and the local difference:

[0023] f(x) = R(x) + C(x) (1)

[0024] In the formula, x is an n-dimensional vector, and its dimension represents the number of input parameters, that is, the number of design variables; R(x) is the regression function representing the global approximate estimate, and R(x) takes a constant value β; C(x) is a correlation function with a mean of 0, a variance of σ 2 , a non-zero covariance, representing the local difference; therefore, the above formula can be rewritten as:

[0025] y(x) = β + C(x) (2)

[0026] Step 4-2: The local difference in the prediction of the unknown quantity x is represented by a stochastic process. The Gaussian random function is used as the correlation function to interpolate the sample points to estimate the trend of the stochastic process; C(x i ) and C(x j ) The correlation between them is related to the distance between the two points x i and x j , i = 1, 2,..., m, j = 1, 2,..., m, and m is the number of sample points used to generate the surrogate model; in the Kriging model, the weighted distance is used:

[0027]

[0028] Among them, are the k-th elements in the i / j-th sample input parameter vector respectively;

[0029] Combined with the Gaussian random function, the correlation function R between x i and x j is expressed as:

[0030]

[0031] Wherein, n is the number of input parameters, and θ, 0 ≤ θ ≤ ∞, is a vector of relevant parameters;

[0032] Step 4-3: The predicted output of Kriging is written as:

[0033]

[0034] Wherein, is the estimated value of β; R is an m×m dimensional correlation matrix, and the matrix element R (i,j) = R(x i , x j ); r is an m-dimensional correlation vector, r i (x) = R(x, x i ); y is an m-dimensional column vector composed of the response values of the sampling sample points, as shown in the following formula:

[0035] y = [y(x1),..., y(x m )] (6)

[0036] The covariance matrix of C(x) is expressed as:

[0037] Cov[C(x i ), C(x j )] = σ 2 R[R(x i , x j )] (7)

[0038] Solve for the unknown parameters, β and σ 2 The estimates in the sense of least squares are obtained from the following two formulas:

[0039]

[0040] The solution of the relevant parameter θ is a one-dimensional optimization problem, as shown in the following formula:

[0041]

[0042] Solve for the value of the relevant parameter θ by maximum likelihood estimation.

[0043] Preferably, the specific content of the step 2-1 is:

[0044] The input of the aerodynamic / thermal surrogate model is the incoming flow Mach number Ma1, the wall panel deformation, and the wall panel surface temperature; the Mach number is directly used as the model input parameter as a scalar, and the wall panel deformation and surface temperature are parameterized through spatial coordinates;

[0045] Specifically, the wall panel deformation w(x) is represented by the superposition of its first six-order modal functions, as shown in formula (11), Φ iis a modal function, h is the thickness of the panel; the surface temperature distribution T(x) is expressed as an nth-order polynomial of the coordinate x in the chordwise direction of the panel, as shown in Equation (12):

[0046]

[0047] T(x) = T1 + T1x + T3x 2 + T4x 3 (12)

[0048] Finally, the panel deformation is parameterized by the dimensionless modal displacement amplitude parameters, and the surface temperature of the panel is parameterized by the polynomial coefficient T i parameters; in addition, displacement amplitude and surface temperature extreme value constraints are set: the dimensionless deformation amplitude of the panel ≤ 10, 200K ≤ the surface temperature extreme value of the panel ≤ 2000K.

[0049] Preferably, the specific content of step 2-2 is as follows:

[0050] According to the parameter boundary and extreme value constraints, use the constrained Latin hypercube sampling method to sample in the parameter space, and select sample points and evaluation points;

[0051] Preferably, the specific content of step 2-3 is as follows:

[0052] Referring to Equation (1), the aerodynamic / thermal surrogate model based on Kriging is expressed as the sum of the local deviation C(d, X) and the global estimate R(d, X), as shown in the following equation:

[0053] y(d) = R(d, X) + C(d, X) (13)

[0054] where y(d) is the prediction of the surrogate model for the physical quantity distribution on the panel surface at a certain state point in the parameter space, that is, under a certain Mach number, a certain deformation, and a certain surface temperature distribution; d is the vector composed of the input parameters corresponding to a certain state point, and X is the input parameter matrix of the training sample points. The elements in d and X are shown in Equations (14) and (15), where n is the number of sample points used to establish the model;

[0055]

[0056] Preferably, the specific content of step 2-4 is as follows:

[0057] Based on the sample point input parameter matrix and the corresponding sample point response data matrix as training data, establish an aerodynamic / thermal surrogate model; use the dimensionless parameter aerodynamic pressure coefficient C P and the Stanton number St as the output of the model, and obtain the aerodynamic pressure and aerodynamic heat flux values through Equations (16) and (17):

[0058]

[0059] Among them, y CP (d) is the predicted value of the wall panel surface pressure coefficient, and ρ3 and U3 are the density and velocity in the leading edge region of the elastic wall panel;

[0060] Q aero = c p U3ρ3(T s (x) - T 0,3 )y St (d) (17)

[0061] Among them, y St (d) is the predicted value of the wall panel surface Stanton number, and c p is the specific heat at constant pressure.

[0062] Preferably, the specific steps of step 2-5 are as follows:

[0063] Based on the idea of cross-validation, 300 evaluation points with non-overlapping sample points used to construct the surrogate model are selected, and the model accuracy is quantitatively evaluated by comparing the prediction results of the surrogate model with the CFD calculation results;

[0064] Root mean square error and maximum error L ∞ , such as:

[0065]

[0066] Among them, i represents the i-th node on the plate, j represents the j-th evaluation point, and p is the total number of nodes on the plate; ROM is the predicted value of the surrogate model, and Full is the value obtained by full-order CFD calculation; based on the above indicators, four values are calculated as the evaluation indicators for quantitatively evaluating the entire surrogate model, namely: average root mean square error Average maximum error Maximum maximum error L ∞max And L ∞ The proportion per of the number of evaluation points where it is less than 10% L∞≤10％ , and the calculation formula is as follows:

[0067]

[0068] Among them, n is the number of evaluation points;

[0069] At the same time, the above four evaluation indicators are tested to verify the convergence of the model as the number of sample points changes.

[0070] Preferably, the specific steps of step 2-6 are as follows:

[0071] Calculate the distribution of the average relative error between the predicted value and the CFD method calculation result in the chordwise position of the wall panel:

[0072]

[0073] For the prediction of unsteady aerodynamic forces, based on the prediction of steady-state aerodynamic forces by the surrogate model, a correction term considering unsteady effects based on the third-order piston theory is combined to establish an unsteady aerodynamic force surrogate model for the thermal protection wall panel;

[0074] The calculation formula for the unsteady pressure coefficient based on the surrogate model is as follows:

[0075]

[0076] where is the steady-state pressure coefficient predicted by the surrogate model, C p,vel and are the aerodynamic damping terms based on the third-order piston theory, and γ is the specific heat ratio of air;

[0077] Finally, the calculation formula for the unsteady aerodynamic pressure on the wall panel surface is:

[0078]

[0079] Preferably, the specific step 3 is as follows:

[0080] Use finite element to solve structural dynamics. The basic motion equation in transient dynamics analysis is:

[0081]

[0082] where M is the mass matrix, C is the damping matrix, K is the stiffness matrix, is the nodal acceleration vector, is the nodal velocity vector, u is the nodal displacement vector, and F is the load matrix varying with time t;

[0083] In the aerodynamic heat module, due to the existence of the radiation layer, the upper surface is a convective boundary, and there is aerodynamic heat flux Q aero input and radiation heat flux Q rad output. Therefore, the net heat flux on the upper surface of the thermal protection wall panel is Q = Q aero -Q rad ;

[0084] The radiation heat flux is related to the temperature difference between the wall surface temperature T wall and the ambient temperature T env The calculation formula is:

[0085]

[0086] where σ is the Stefan-Boltzmann constant, taking σ = 5.669×10 -8 / (W / m2 / T 4 ), the surface of the structure is non - bold, and the radiation emissivity ε is taken as 0.7;

[0087] Thermal conduction is considered inside the structure. The transient heat conduction partial differential equation of the two - dimensional panel is as follows:

[0088]

[0089] Among them, x is the chordwise direction of the panel, z is the thickness direction of the panel, ρ is the material density, c is the specific heat ratio of the material, T is the temperature, k is the thermal conductivity of the material, and t is the time; Central difference is used for spatial discretization, and forward difference is used for time marching. The discrete algebraic equation:

[0090]

[0091] Among them, i and j are the node numbers in the x and z directions respectively, numbered from left to right and from top to bottom; The leading and trailing edge end faces and the lower surface are adiabatic boundaries, and their temperatures are equal to the corresponding solid boundary temperatures.

[0092] The beneficial effects of the present invention are as follows:

[0093] 1) By simplifying the thin - wall structure on the surface of a typical high - speed aircraft, using a three - layer thermal protection panel as the physical model, establishing the aeroelasticity, aerothermal, and heat conduction models, and performing coupled modeling between multiple physical fields, aero - thermo - elastic analysis technology for the thermal protection panel structure is formed;

[0094] 2) To solve the contradiction between the calculation accuracy and efficiency of aerodynamic force / heat, the Kriging model is used to form the aerodynamic force / heat surrogate model technology for the panel in high - speed flow;

[0095] 3) The aerodynamic force / heat surrogate model established based on the Kriging method has a certain generalization ability and still has accuracy for predictions outside the state space;

[0096] 4) The calculation efficiency of the aerodynamic force / heat surrogate model is 4 orders of magnitude higher than that of traditional high - order methods (such as CFD). Brief Description of the Drawings

[0097] Figure 1 It is a schematic diagram of the structure of a typical high - speed aircraft;

[0098] Figure 2 It is the TPS thermal protection panel structure on the inclined surface of a scramjet engine;

[0099] Figure 3 It is the flow chart for establishing the surrogate model;

[0100] Figure 4 It is the flow chart for evaluating the model accuracy;

[0101] Figure 5 Schematic diagram of the thermo-pneumatoelastic tetrahedron model;

[0102] Figure 6 Schematic diagram of the coupling strategy for thermo-pneumatoelastic analysis;

[0103] Figure 7 Schematic diagram of the basic idea of the surrogate model;

[0104] Figure 8 Schematic diagram of the high-speed thermal protection wall panel;

[0105] Figure 9 Flow chart of the thermo-pneumatoelastic analysis of the thermal protection wall panel based on the aerodynamic / thermal surrogate model;

[0106] Figure 10 Schematic diagram of the variation of the evaluation index of the aerodynamic surrogate model with the number of samples;

[0107] Figure 11 Schematic diagram of the position distribution of the average relative error of the aerodynamic surrogate model;

[0108] Figure 12 Schematic diagram of the unsteady generalized aerodynamic force of the moving wall panel;

[0109] Figure 13 Schematic diagram of the variation of the aerothermal evaluation index with the number of samples;

[0110] Figure 14 Schematic diagram of the position distribution of the average relative error of the aerothermal surrogate model;

[0111] Figure 15 Time history diagram;

[0112] Figure 16 Schematic diagram of the wall panel displacement at the moment of response stability;

[0113] Figure 17 Schematic diagram of the thermal deformation of the thermal protection wall panel caused by the temperature field;

[0114] Figure 18 Diagram of the limit cycle motion response, (a) time history diagram, (b) phase plane diagram;

[0115] Figure 19 Diagram of the period-doubling motion response, (a) time history diagram, (b) frequency spectrum diagram, (c) phase plane diagram; Figure 20 Diagram of the quasi-periodic motion response, (a) time history diagram, (b) frequency spectrum diagram, (c) phase plane diagram;

[0116] Figure 21 Diagram of the chaotic motion response, (a) time history diagram, (b) frequency spectrum diagram, (c) phase plane diagram. Detailed implementation manners

[0117] The present invention will be further described below in conjunction with the accompanying drawings and embodiments.

[0118] The surface of an aircraft flying at high speed will be subjected to a complex aerodynamic / thermal environment, resulting in thermal stress in the structure and material degradation, and thermo-aeroelastic problems occurring under the coupling relationship of the three physical fields of fluid-structure-thermal. On the other hand, for the aerodynamic / thermal calculation method under high-speed flow, its engineering algorithm has high calculation efficiency but insufficient accuracy, while the numerical simulation method improves the calculation accuracy but has high calculation costs. In view of the contradiction between the aerodynamic / thermal calculation accuracy and efficiency, the present invention has developed a fast thermo-aeroelastic calculation technology for the thermal protection panel of a high-speed aircraft based on an aerodynamic / thermal surrogate model. It realizes accurate and fast prediction of high-speed aerodynamic / thermal and fast calculation of the structural thermo-aeroelastic response. It can further provide theoretical guidance for the thermal protection structure design of high-speed aircraft, the ballistic design optimization under thermal constraints, flight safety assessment, etc. The surrogate model idea is as Figure 7 shown.

[0119] The present invention takes the TPS wall panel on the surface of the forebody-inlet of a general high-speed aircraft as the research object ( Figure 2 ), which is composed of a thermal protection layer and a bottom metal plate. The structural schematic diagram of the thermal protection wall panel is as Figure 8 shown. Among them, the thermal protection layer includes a radiation layer and a heat insulation layer, which are respectively used to radiate heat to the environment and block the downward conduction of heat. The two ends of the bottom metal plate are simply supported. When the aircraft is under a high-speed oncoming flow (zone 1), an oblique shock wave is generated at the leading edge of the sharp wedge, and the airflow behind the wave (zone 2) flows parallel to the surface of the sharp wedge. The gas parameters at the leading edge of the elastic wall panel (zone 3) are the same as those in zone 2. The airflow in zone 2 can be obtained through the oblique shock wave theory and its relational expressions.

[0120] Based on the Kriging surrogate model method, an unsteady aerodynamic / thermal surrogate model of the thermal protection wall panel in high-speed flow is established. The aerodynamic elastic module is established using the structural dynamics method based on the finite element method, and the aerothermal module is established using the transient heat conduction equation to build a thermo-aeroelastic analysis framework for the thermal protection wall panel. Finally, the aerodynamic / thermal surrogate model is used to carry out the thermo-aeroelastic calculation and analysis of the thermal protection wall panel.

[0121] The most direct method to obtain the input-output relationship of the system is to reproduce physical phenomena or conduct experiments, but the complexity of the problem or the limitations of experimental conditions often limit the acquisition of data. Therefore, by establishing a mathematical model of the system and solving the mathematical model to describe the input-output relationship of the original system has become a more commonly used means. However, the accuracy of the mathematical model depends on physical assumptions and prior knowledge. At the same time, the complexity and high cost of mathematical solution reduce the engineering practicability.

[0122] In contrast, the surrogate model does not require prior knowledge. Based on mathematical theory, a mathematical model similar to the output of the original system or the mathematical model is constructed, thus "replacing" or acting as a surrogate for the original system or the mathematical model. Moreover, it usually has a simple form and strong timeliness, can be well combined with other discipline models, and significantly reduces the input of computing resources, with remarkable advantages.

[0123] The present invention uses the Kriging model to establish a surrogate model. Compared with machine learning methods, etc., the parameters of this method are all obtained by mathematical regression, without parameter adjustment and long-time data training. The mathematical expression of the Kriging model can be written as a linear superposition of the global estimate and the local difference:

[0124] f(x) = R(x) + C(x) (1)

[0125] In the formula, x is an n-dimensional vector, and its dimension represents the number of input parameters, that is, the number of design variables; R(x) is the regression function representing the global approximate estimate, and generally R(x) takes a constant value β; C(x) is a correlation function with a mean of 0, a variance of σ 2 , a non-zero covariance, representing the local difference. Therefore, the above formula can be written as:

[0126] y(x) = β + C(x) (2)

[0127] In the model, the local difference in predicting the unknown quantity x is represented by a stochastic process. A Gaussian random function is used as the correlation function to interpolate the sample points, so as to estimate the trend of the stochastic process. The correlation between C(x i ) and C(x j ) is related to the distance between the two points x i and x j , i = 1, 2,..., m, j = 1, 2,..., m, and m is the number of sample points used to generate the surrogate model. Since the Euclidean distance has the same weight for all sample points, in the Kriging model, a special weighted distance is used:

[0128]

[0129] Combined with the Gaussian random function, the correlation function R between x i and x j is expressed as:

[0130]

[0131] In the formula, n is the number of input parameters, and θ (0 ≤ θ ≤ ∞) is the correlation parameter vector.

[0132] The predicted output of Kriging can be written as:

[0133]

[0134] In the formula, is the estimated value of β; R is an m×m dimensional correlation matrix, and the matrix element R (i,j) = R(x i , x j ); r is an m-dimensional correlation vector, r i (x) = R(x, x i ); y is an m-dimensional column vector composed of the response values of the sampling sample points, as shown in the following formula:

[0135] y = [y(x1),..., y(x m )] (6)

[0136] The covariance matrix of C(x) can be expressed as:

[0137] Cov[C(x i ), C(x j )] = σ 2 R[R(x i , x j )] (7)

[0138] Solving for the unknown parameters, the estimates of β and σ 2 in the sense of least squares are obtained from the following two formulas:

[0139]

[0140] The solution of the relevant parameter θ is a one-dimensional optimization problem, as shown in the following formula:

[0141]

[0142] The value of the relevant parameter θ is solved by Maximum Likelihood Estimation.

[0143] 1. Establish an aerodynamic / thermal surrogate model elastic model

[0144] Figure 3 shows the general process of establishing a surrogate model, including:

[0145] ① Determine the state space of the system: Parametrize the input and output variables of the system to obtain the expected input and output parameters of the surrogate model, and determine the parameter boundaries;

[0146] ② Sample the parameter state space: Generally use the Latin hypercube sampling method to select a certain number of state points in the state space as samples and validations;

[0147] ③ Response of sample points: Use a high-order method to calculate the working conditions corresponding to the sample points to generate training data;

[0148] ④ Generate a surrogate model: Based on the input-output of the sample state points and response results, establish a surrogate model;

[0149] ⑤ Precision evaluation: Evaluate the precision of the surrogate model at the evaluation points according to the cross-validation idea.

[0150] ⑥ Apply the surrogate model for prediction: If the requirements are met, the surrogate model can be used.

[0151] The inputs of the aerodynamic / thermal surrogate model are the selected incoming flow Mach number Ma1, the wall panel deformation, and the wall panel surface temperature. The Mach number, as a scalar, can be directly used as an input parameter of the model, and the wall panel deformation and surface temperature are parameterized through spatial coordinates.

[0152] Specifically, the wall panel deformation w(x) is represented by the superposition of its first six-order modal functions, as shown in Equation (11), where Φ i is the modal function and h is the wall panel thickness; the surface temperature distribution T(x) is expressed by an n-order polynomial of the coordinate x in the chord direction of the wall panel, as shown in Equation (12).

[0153]

[0154] T(x) = T1 + T1x + T3x 2 + T4x 3 (12)

[0155] Finally, the wall panel deformation is parameterized by the dimensionless modal displacement amplitude parameter, and the wall panel surface temperature is parameterized by the polynomial coefficient T i parameter. The parameter space boundaries are shown in Table 1. In addition, set the extreme value constraints for the displacement amplitude and surface temperature: the dimensionless deformation amplitude of the wall panel ≤ 10, 200K ≤ the extreme value of the wall panel surface temperature ≤ 2000K.

[0156] Table 1 Surrogate model parameter space

[0157]

[0158] According to the parameter boundaries and extreme value constraints, use the constrained Latin hypercube sampling method to sample in the parameter space and select the sample points and evaluation points.

[0159] Referring to Equation (1), the aerodynamic / thermal surrogate model based on Kriging can be expressed as the sum of the local deviation C(d, X) and the global estimate R(d, X), as shown in the following Equations (14) and (15):

[0160] y(d) = R(d, X) + C(d, X) (13)

[0161] Among them, y(d) is the prediction of the physical quantity distribution on the panel surface by the surrogate model at a certain state point in the parameter space, that is, under a certain Mach number, a certain deformation, and a certain surface temperature distribution. d is the vector composed of the input parameters corresponding to a certain state point, and X is the input parameter matrix of the training sample points. The elements in d and X are shown in the following two equations, where n is the number of sample points used to establish the model.

[0162]

[0163] Based on the sample point input parameter matrix and the corresponding sample point response data matrix as training data, an aerodynamic / thermal surrogate model is established. To reduce the interference of the data magnitude on the model, the dimensionless parameters aerodynamic pressure coefficient C P and Stanton number St are used as the outputs of the model. The numerical values of aerodynamic pressure and aerodynamic heat flux are obtained through equations (16) and (17).

[0164]

[0165] Among them, y CP (d) is the predicted value of the pressure coefficient on the panel surface, and ρ3 and U3 are the density and velocity in zone 3.

[0166] Q aero =c p U3ρ3(T s (x)-T 0,3 )y St (d) (17)

[0167] Among them, y St (d) is the predicted value of the Stanton number on the panel surface, ρ3 and U3 are the density and velocity in zone 3, and c p is the specific heat at constant pressure.

[0168] Based on the idea of cross-validation, 300 evaluation points with non-overlapping sample points used to construct the surrogate model are selected. By comparing the prediction results of the surrogate model with the CFD calculation results, the accuracy of the model is quantitatively evaluated. The surrogate model accuracy evaluation process is as Figure 4 shown.

[0169] The root mean square error (Normalized Root Mean Square Error, NRMSE) and the maximum error L ∞ , such as:

[0170]

[0171] Among them, i represents the i-th node on the panel, j represents the j-th evaluation point, and p is the total number of nodes on the panel; ROM is the predicted value of the surrogate model, and Full is the value obtained by full-order CFD calculation. Based on the above indicators, four values are calculated as evaluation indicators for quantitatively evaluating the entire surrogate model, namely: root mean square error Mean maximum error Maximum maximum error (L ∞max ) and the proportion of the number of evaluation points where L ∞ is less than 10% (per L∞≤10％ ), and the calculation formula is as follows:.

[0172]

[0173]

[0174] Among them, n is the number of evaluation points. At the same time, the above four evaluation indicators are tested to verify the convergence of the model as the number of sample points changes.

[0175] In addition, to illustrate the prediction accuracy of the surrogate model for different positions of the panel, the distribution of the average relative error between the predicted value and the CFD calculation result in the chordwise position of the panel is calculated by the following formula:

[0176]

[0177] For the prediction of unsteady aerodynamic forces, based on the prediction of steady aerodynamic forces by the surrogate model, a correction term considering unsteady effects based on the third-order piston theory is combined to establish an unsteady aerodynamic force surrogate model for the thermal protection panel.

[0178] The calculation formula for the unsteady pressure coefficient based on the surrogate model is as follows:

[0179]

[0180] Among them is the steady pressure coefficient predicted by the surrogate model, C p,vel and are the aerodynamic damping terms based on the third-order piston theory, and γ is the specific heat ratio of air.

[0181] Finally, the calculation formula for the unsteady aerodynamic pressure on the panel surface is:

[0182]

[0183] 2. Establish a thermo-aeroelastic coupling model;

[0184] The thermo-aeroelastic problem is a typical fluid-structure-thermal coupling problem, including the coupling problems of aerodynamics, heat, and structure, and simultaneously involves the coupled interactions of inertial forces, elastic forces, aerodynamic forces, and thermal effects.Figure 5 The force tetrahedron model for the thermo-pneumatoelastic problem is given.

[0185] The thermo-pneumatoelastic coupling analysis method decomposes the three-field coupling system into two modules, namely two coupled subsystems, the thermal module and the pneumatoelastic module. Among them, the thermal module is a coupled system of aerothermal and structural transient heat conduction. The aerothermal calculation obtains the heat flux and transmits it to the heat conduction module. The heat conduction calculation obtains the transient temperature distribution of the structure and then transmits it to the aerothermal calculation to update the heat flux at the next time step; the pneumatoelastic module calculates the unsteady aerodynamic force by the aerodynamic force model and transmits it to the structural dynamics module. The elastic deformation is calculated by the structural dynamics model and then transmitted to the aerodynamic force calculation to update the aerodynamic force at the next time step.

[0186] The aerothermal-pneumatoelastic coupling strategy of the present invention is that the structural temperature field obtained by the thermal module is transmitted to the pneumatoelastic module. The coupling strategy process is as Figure 6 shown.

[0187] In the aeroelastic module, the finite element method is used for structural dynamics solution. The basic motion equation in transient dynamics analysis is:

[0188]

[0189] where M is the mass matrix, C is the damping matrix, K is the stiffness matrix, is the nodal acceleration vector, is the nodal velocity vector, u is the nodal displacement vector, and F is the load matrix that varies with time t.

[0190] In the aerothermal module, due to the presence of the radiation layer, the upper surface is a convective boundary, and there is an aerothermal flux Q aero input and a radiation heat flux Q rad output. Therefore, the net heat flux on the upper surface of the thermal protection wall panel is Q = Q aero -Q rad .

[0191] The radiation heat flux is related to the temperature difference between the wall surface temperature T wall and the ambient temperature T env . The calculation formula is:

[0192]

[0193] where σ is the Stefan-Boltzmann constant, taking σ = 5.669×10 -8 / (W / m 2 / T 4 ), the structure surface is non-blackbody, and the radiation emissivity ε is taken as 0.7.

[0194] Considering heat conduction inside the structure, the transient heat conduction partial differential equation of the two-dimensional wall panel is as follows:

[0195]

[0196] Among them, x is the chordwise direction of the panel, z is the thickness direction of the panel, ρ is the material density, c is the specific heat ratio of the material, T is the temperature, k is the thermal conductivity of the material, and t is the time. Central difference is used for spatial discretization, and forward difference is used for time marching. The discrete algebraic equations are:

[0197]

[0198] Among them, i and j are in the x and z directions respectively, and are the node numbers from left to right and from top to bottom. The leading and trailing edge end faces and the lower surface are adiabatic boundaries, and their temperatures are equal to the corresponding solid boundary temperatures.

[0199] The flowchart of the thermo-aeroelastic analysis of the thermal protection panel based on the aerodynamic / thermal surrogate model is as Figure 9 shown.

[0200] Example 1: Aerodynamic surrogate model of high-speed flow panel

[0201] Table 2 gives the accuracy evaluation indexes of the aerodynamic surrogate models established using different numbers of sample points.

[0202] Table 2 Evaluation indexes of aerodynamic surrogate models with different sample numbers

[0203]

[0204] Figure 10 Shows the changes of each with the increase of the sample number. It can be seen that as the sample number gradually increases from 200 to 1500, the accuracy evaluation indexes of the surrogate model no longer improve, verifying the convergence of the sample number used to generate the surrogate model. Figure 11 Shows the distribution of the relative error between the calculated results of the aerodynamic pressure coefficient of the aerodynamic surrogate model generated using 1500 sample points and the CFD results in the chordwise position of the panel.

[0205] Use the aerodynamic surrogate model to predict the unsteady generalized aerodynamic force (GAF) of the panel with a specified motion form, and verify the accuracy of the unsteady aerodynamic model based on the surrogate model.

[0206] The calculation formula of GAF is:

[0207] GAF(t) = ∫(Φ(x)q a (x,t))dx (33)

[0208] Among them, Φ(x) is the modal shape. The vibration equation of the panel is as shown in Equation (34), and the frequency f = 140Hz.

[0209]

[0210] The first-order mode is selected as the motion form of the panel.

[0211] Calculation conditions: incoming flow Mach 10.0, altitude 30 km, and the surface temperature of the panel is uniformly distributed at 1200 K.

[0212] Figure 12 The results of calculating the unsteady generalized aerodynamic force of the moving panel by various methods are shown. The results show that the results of calculating the unsteady generalized aerodynamic force by the surrogate model are consistent with the CFD unsteady calculation results, verifying the accuracy of calculating the unsteady aerodynamic force based on the surrogate model.

[0213] It is noted that the maximum non-dimensional amplitude of the panel vibration is 12, exceeding the state space of the surrogate model. And its unsteady GAF prediction results are consistent with the CFD unsteady calculation results, indicating that the established unsteady aerodynamic force model based on the surrogate model has a certain generalization ability, that is, when predicting the state points outside the state space, it still has a certain prediction accuracy.

[0214] Figure 12 The calculation results using the unsteady third-order piston theory are also shown, which have obvious differences from the CFD unsteady calculation results, demonstrating the importance of establishing a surrogate model with high-precision prediction once again.

[0215] In terms of calculation efficiency, when using the third-order unsteady piston theory and the unsteady aerodynamic force surrogate model, the time consumed for calculating the verification conditions is within 0.5 s. While using CFD for unsteady aerodynamic force calculation, the calculation time is more than 8 hours (CPU configuration: Gold 6152 2.10 GHz * 12), which is 60,000 times that of the surrogate model. Therefore, compared with the CFD method, the surrogate model method improves the calculation efficiency by four orders of magnitude.

[0216] Example 2: Aerodynamic heat surrogate model of high-speed flow panel;

[0217] Table 3 gives the accuracy evaluation indexes of the aerodynamic heat surrogate models established using different numbers of sample points.

[0218] Table 3 Evaluation indexes of aerodynamic heat surrogate models with different sample numbers

[0219]

[0220] Figure 13 The variations with the increase in the number of sample points are shown. It can be seen that as the number of sample points gradually increases from 200 to 2400, the accuracy evaluation indexes of the surrogate model no longer improve, verifying the convergence of the number of sample points used to generate the surrogate model.

[0221] Figure 14 Shows the distribution of the relative error between the Stanton number calculation results of the aerodynamic force surrogate model generated using 2200 sample points and the CFD results in the chordwise position of the panel.

[0222] Example 3: Thermal Aeroelastic Analysis of High-Speed Flow Panels;

[0223] Using the present invention, two configurations, namely the panel deformation snapshots under typical dynamic response conditions and the rigid body configuration of the panel (i.e., flat plate), are considered as the inputs to the aerothermal module. Using the aerothermal surrogate model, transient heat conduction calculations of the thermal protection panel structure are performed in the aerothermal module. The temperature field snapshot at a certain moment is used as the thermal load input to the aeroelastic solution module. Using the aerodynamic force surrogate model, thermal aeroelastic calculations are carried out for typical deformation snapshots:

[0224] Table 4 Thermal Aeroelastic Calculation Conditions of Thermal Protection Panels

[0225]

[0226]

[0227] In this example, thermal aeroelastic calculations are carried out for four groups of conditions shown in Table 4, where the oncoming flow conditions are all the parameters Ma1 in Zone 1. To simulate the scenario of a sudden increase in dynamic pressure when a high-speed aircraft performs maneuvers such as penetration maneuvers, the oncoming flow Mach numbers for the heat conduction calculations in Conditions 2 and 3 are the same, while the aeroelastic calculations are different.

[0228] The four groups of thermal aeroelastic conditions result in the following typical nonlinear response forms such as buckling, limit cycle motion, period-doubling motion, quasi-periodic motion, and chaotic motion.

[0229] Buckling:

[0230] In Condition 2, under the temperature fields at 60 s, 100 s, and 180 s of heat transfer, the thermal aeroelastic response forms of the panel are all thermal buckling. Figure 15 Shows the time history diagram of the panel response under the temperature field at 60 s.

[0231] Figure 16 Shows the deflection curve of the panel displacement deformation after the thermal buckling response converges to stability under the temperature fields at 60 s, 100 s, and 180 s. Figure 17 Shows the static deformation of the thermal stress under these three temperature fields. It can be seen that as the heat conduction of the thermal protection panel progresses, the higher the temperature inside the temperature field panel, the greater the thermal stress caused by thermal expansion, which also leads to a greater extreme value of the thermal buckling displacement.

[0232] Limit cycle motion:

[0233] In working condition 3, at the temperature field at the 300 s moment of heat transfer, the thermo-aeroelastic response form of the panel is a limit cycle.

[0234] In the time history diagram of the limit cycle motion, the vibration form of the panel shows a certain periodic law. Figure 18 In it, the response curves in some time periods are enlarged for intuitive display. Correspondingly, there is only a single closed phase plane curve in the phase plane diagram.

[0235] Period-doubling motion:

[0236] In working condition 3, at the temperature field at the 100 s moment of heat transfer, the thermo-aeroelastic response form of the panel is period-doubling motion.

[0237] In the time history diagram of the period-doubling motion, the vibration form of the panel shows a certain periodic law. Figure 19 In it, the same enlargement is carried out, but the vibration form is more complex than that of the limit cycle. The phase plane diagram is in the form of a band formed by multiple closed curves. The main frequencies corresponding to the peaks in the frequency spectrum diagram are 10.22 Hz, 20 Hz, 30.22 Hz, and 40 Hz, showing an integer multiple relationship of 2 times, 3 times, and 4 times.

[0238] Quasi-periodic motion:

[0239] In working condition 4, at the temperature field at the 300 s moment of heat transfer, the thermo-aeroelastic response form of the panel is quasi-periodic motion.

[0240] Figure 20 The time history diagram, frequency spectrum diagram, and phase plane diagram of the panel response are shown. The time history diagram and phase plane diagram of the quasi-periodic motion are similar to those of the period-doubling motion. The vibration of the panel shows a certain periodicity (enlarged in the time history diagram), and the phase plane diagram is in the form of a band formed by multiple closed curves. However, there is no integer multiple relationship between the frequencies corresponding to the peaks in the frequency spectrum diagram.

[0241] Chaotic motion:

[0242] In working condition 1, at the temperature field at the 100 s moment of heat transfer, the thermo-aeroelastic response form of the panel is chaotic motion.

[0243] Figure 21 The time history diagram, frequency spectrum diagram, and phase plane diagram of the panel response are shown. In the chaotic response, the vibration of the panel shown in the time history diagram has no obvious law, the distribution of the curves in the phase plane diagram is chaotic, and multiple peaks in the frequency spectrum diagram are continuously distributed within a certain frequency range.

[0244] The present invention takes into account the non-linear characteristics caused by the structural thermal stress and material property degradation induced by the aerodynamic heating effect of the high-speed oncoming flow, and uses the aerodynamic / thermal surrogate model to efficiently and accurately predict and analyze the thermo-aeroelasticity of the thermal protection wall panel. It can provide an important theoretical basis for the accurate prediction of aerodynamic / thermal loads, the design of thermal protection structures, and the flight safety assessment of high-speed aircraft.

[0245] Based on the fast calculation technology of the present invention, the time consumption for calculating the thermo-aeroelasticity of the thermal protection wall panel using the aerodynamic / thermal surrogate model is about 30 minutes. Compared with the time consumption of dozens of hours using high-order methods (such as CFD / CSD), the efficiency of the thermo-aeroelasticity calculation has been greatly improved.

Claims

1. A rapid calculation method for thermal pneumoelasticity of the thermal protection panel of a high-speed aircraft based on an aerodynamic / thermal surrogate model, characterized in that It includes the following steps: Step 1: Set the research object; Take the TPS wall panel on the forebody-inlet surface of a high-speed aircraft as the research object. The TPS wall panel consists of a thermal protection layer and a bottom metal plate. Among them, the thermal protection layer includes a radiation layer and a heat insulation layer, which are used to radiate heat to the environment and block the heat from conducting downward respectively; the two ends of the bottom metal plate are simply supported. When the aircraft is under high-speed oncoming flow, an oblique shock wave is generated at the sharp wedge leading edge, and the airflow behind the wave flows parallel to the sharp wedge surface. The gas parameters at the leading edge of the elastic wall panel are the same as those in the airflow region behind the wave; Step 2: Based on the Kriging surrogate model method, establish an unsteady aerodynamic / thermal surrogate model for the thermal protection wall panel in high-speed flow; Step 2-1: Determine the state space of the system: Parametrize the input and output variables of the system to obtain the input and output parameters expected by the surrogate model, and determine the parameter boundaries; Step 2-2: Sample the parameter state space: Use the Latin hypercube sampling method to select a certain number of state points in the state space as samples and for verification; Step 2-3: Response of sample points: Use a high-order method to calculate the working conditions corresponding to the sample points to generate training data; Step 2-4: Generate a surrogate model: Based on the input and output of the sample state points and response results, establish a surrogate model; Step 2-5: Accuracy evaluation: Evaluate the accuracy of the surrogate model at the evaluation points according to the cross-validation idea; Step 2-6: Apply the surrogate model for prediction: If the requirements are met, the surrogate model can be used; Step 3: Establish a thermo-aeroelastic coupling model; Adopt the thermo-aeroelastic coupling analysis method to decompose the three-field coupling system into two modules, namely two coupled subsystems: the thermal module and the aeroelastic module. Among them, the thermal module is a coupled system of aerothermal and structural transient heat conduction. The heat flux calculated by aerothermal is transmitted to the heat conduction module, and the transient temperature distribution of the structure obtained by heat conduction calculation is then transmitted to the aerothermal calculation to update the heat flux at the next time step; the aeroelastic module calculates the unsteady aerodynamic force by the aerodynamic force model and transmits it to the structural dynamics module, and calculates the elastic deformation by the structural dynamics model and then transmits it to the aerodynamic force calculation to update the aerodynamic force at the next time step.

2. A rapid computational method for thermo-aeroelasticity of a thermal protection panel of a high-speed aircraft based on an aerodynamic / thermal surrogate model, characterized in that, The specific expression of the Kriging model is as follows: Step 4-1: Write the mathematical expression of the Kriging model as a linear superposition of the global estimate and the local difference: f(x) = R(x) + C(x) (1) where x is an n-dimensional vector, whose dimension represents the number of input parameters, that is, the number of design variables; R(x) is a regression function representing the global approximate estimation, and R(x) takes a constant value β; C(x) is a correlation function with a mean of 0, a variance of σ 2 , a non-zero covariance, which characterizes the local differences; therefore, the above equation can be rewritten as: y(x) = β + C(x) (2) Step 4-2: Represent the local difference in the prediction of the unknown quantity x by a stochastic process, use a Gaussian random function as the correlation function to interpolate the sample points, and thereby estimate the trend of the stochastic process; The correlation between C(x i ) and C(x j ) is related to the distance between the two points x i and x j , where i = 1, 2,..., m, j = 1, 2,..., m, and m is the number of sample points used to generate the surrogate model; In the Kriging model, the weighted distance is used: wherein, are the k-th elements in the i / j-th sample input parameter vector, respectively; Combined with the Gaussian random function, x i and x j The correlation function R between them is expressed as: In the formula, n is the number of input parameters, θ, 0 ≤ θ ≤ ∞ is the correlation parameter vector; Step 4-3: Write the predicted output of Kriging as: In the formula, is the estimated value of β; R is an m×m dimensional correlation matrix, and the matrix element R (i,j) = R(x i , x j ); r is an m-dimensional correlation vector, r i (x) = R(x, x i ); y is an m-dimensional column vector composed of the response values of the sampling sample points, as shown in the following formula: y = [y(x1),...,y(x m )] (6) The covariance matrix of C(x) is expressed as: Cov[C(x i ),C(x j )]=σ 2 R[R(x i ,x j )] (7) Solve for the unknown parameters, β and σ 2 The estimates in the least squares sense are obtained from the following two equations: The solution of the correlation parameter θ is a one-dimensional optimization problem, as shown in the following formula: Solve the value of the correlation parameter θ by maximum likelihood estimation.

3. A rapid calculation method for thermo-aeroelasticity of a thermal protection panel of a high-speed aircraft based on an aerodynamic / thermal proxy model according to claim 2, characterized in that The specific content of Step 2-1 is as follows: The inputs of the aerodynamic / thermal surrogate model are the incoming flow Mach number Ma1, the deformation of the wall panel, and the surface temperature of the wall panel; the Mach number is directly used as the model input parameter as a scalar, and the deformation and surface temperature of the wall panel are parameterized through spatial coordinates; Specifically, the deformation w(x) of the panel is represented by the superposition of its first six-order modal functions, as shown in Equation (11), where Φ i is the modal function and h is the panel thickness; the surface temperature distribution T(x) is expressed as an n-order polynomial of the coordinate x in the chord direction of the panel, as shown in Equation (12): T(x) = T1 + T1x + T3x 2 + T4x 3 (12) Finally, the panel deformation is parameterized by the dimensionless modal displacement amplitude and the panel surface temperature is parameterized by the polynomial coefficients T i In addition, the extreme value constraints of the displacement amplitude and the surface temperature are set: the dimensionless deformation amplitude of the panel ≤ 10, 200K ≤ the extreme value of the panel surface temperature ≤ 2000K.

4. A rapid calculation method for thermo-aeroelasticity of a thermal protection panel of a high-speed aircraft based on an aerodynamic / thermal proxy model, characterized in that, The specific content of Step 2-2 is as follows: According to the parameter boundaries and extreme value constraints, the constrained Latin hypercube sampling method is used to sample in the parameter space to select sample points and evaluation points.

5. A rapid calculation method for thermo-aeroelasticity of a thermal protection panel of a high-speed aircraft based on an aerodynamic / thermal proxy model, characterized in that, The specific steps of step 2-3 are as follows: Referring to Equation (1), the aerodynamic / thermal surrogate model based on Kriging is expressed as the sum of the local deviation C(d, X) and the global estimate R(d, X), as shown in the following equation: y(d) = R(d, X) + C(d, X) (13) Among them, y(d) is the prediction of the surrogate model for the physical quantity distribution on the panel surface under a certain state point in the parameter space, that is, under a certain Mach number, a certain deformation, and a certain surface temperature distribution; d is the vector composed of the input parameters corresponding to a certain state point, and X is the input parameter matrix of the training sample points; the elements in d and X are shown in Equations (14) and (15), where n is the number of sample points used to establish the model; 6. A rapid calculation method for thermal pneumoelasticity of a thermal protection panel of a high-speed aircraft based on an aerodynamic / thermal surrogate model, characterized in that, The specific steps of step 2-4 are as follows: Based on the sample point input parameter matrix and the corresponding sample point response data matrix as training data, an aerodynamic / thermal surrogate model is established; the dimensionless parameters aerodynamic pressure coefficient C P and Stanton number St are used as the outputs of the model, and the numerical values of aerodynamic pressure and aerodynamic heat flux are obtained through equations (16) and (17): Among them, is the predicted value of the surface pressure coefficient of the panel, and ρ3 and U3 are the density and velocity in the leading edge region of the elastic panel; Q aero = c p U3ρ3(T s (x) - T 0,3 )y St (d) (17) where y St (d) is the predicted value of the Stanton number on the surface of the wall panel, and c p is the specific heat at constant pressure.

7. A rapid calculation method for thermal pneumoelasticity of a thermal protection panel of a high-speed aircraft based on an aerodynamic / thermal proxy model, characterized in that, The specific steps of step 2-5 are as follows: Based on the idea of cross-validation, 300 evaluation points that do not overlap with the sample points used to construct the surrogate model are selected. By comparing the prediction results of the surrogate model with the CFD calculation results, the accuracy of the model is quantitatively evaluated; Root mean square error and maximum error L ∞ , such as: Among them, i represents the i-th node on the plate, j represents the j-th evaluation point, and p is the total number of nodes on the plate; ROM is the predicted value of the surrogate model, and Full is the value obtained by full-order CFD calculation. Based on the above indicators, four values are calculated as evaluation indicators for quantifying the entire surrogate model, namely: root mean square error mean maximum error maximum maximum error L ∞max and L ∞ proportion of the number of evaluation points less than 10% The calculation formulas are as follows: Among them, n is the number of evaluation points; At the same time, the above four evaluation indicators are tested to verify the convergence of the model as the number of sample points changes.

8. A rapid computational method for thermal pneumoelasticity of a thermal protection panel of a high-speed aircraft based on an aerodynamic / thermal surrogate model, characterized in that, The specific steps of step 2-6 are as follows: Calculate the distribution of the average relative error between the predicted value and the CFD calculation result in the chordwise position of the panel: For the prediction of unsteady aerodynamic forces, based on the prediction of steady-state aerodynamic forces by the surrogate model, a correction term considering unsteady effects based on the third-order piston theory is combined to establish an unsteady aerodynamic force surrogate model for the thermal protection panel; The calculation formula for the unsteady pressure coefficient based on the surrogate model is as follows: Among them is the steady-state pressure coefficient predicted by the surrogate model, C p,vel and is the aerodynamic damping term based on the third-order piston theory, and γ is the specific heat ratio of air; Finally, the calculation formula for the unsteady aerodynamic pressure on the panel surface is:

9. A rapid calculation method for thermal pneumoelasticity of a thermal protection panel of a high-speed aircraft based on an aerodynamic / thermal surrogate model, characterized in that The specific steps of step 3 are as follows: Use finite elements to solve structural dynamics. The basic motion equation in transient dynamics analysis is: where M is the mass matrix, C is the damping matrix, K is the stiffness matrix, ü is the nodal acceleration vector, is the nodal velocity vector, u is the nodal displacement vector, and F is the load matrix varying with time t; In the aerodynamic heat module, due to the presence of the radiation layer, the upper surface is a convective boundary, and there is an aerodynamic heat flux Q aero input and a radiative heat flux Q rad output. Therefore, the net heat flux on the upper surface of the thermal protection wall panel is Q = Q aero -Q rad ; Radiant heat flux is related to the temperature difference between the wall temperature T wall and the ambient temperature T env and is calculated by the formula: where σ is the Stefan-Boltzmann constant, taking σ = 5.669×10 -8 / (W / m 2 / T 4 ), the surface of the structure is non-blackbody, and the radiation emissivity ε is taken as 0.7; Considering heat conduction inside the structure, the transient heat conduction partial differential equation of the two-dimensional panel is as follows: Among them, x is the chordwise direction of the panel, z is the thickness direction of the panel, ρ is the material density, c is the specific heat ratio of the material, T is the temperature, k is the material thermal conductivity, and t is the time; central difference is used for spatial discretization, and forward difference is used for time advancement. The discretized algebraic equation: Among them, i and j are the node numbers in the x and z directions, from left to right and from top to bottom; the front and rear edge end faces and the lower surface are adiabatic boundaries, and their temperatures are equal to the corresponding solid boundary temperatures.

Images