A design method applied to an ablative thermal protection system

Abstract

The present invention discloses a design method applied to an ablative thermal protection system. Quantify various uncertainty random parameters for the considered uncertainty input parameters and use sampling techniques to generate random input samples; then perform deterministic aerodynamic heat and heat conduction analyses on the random input samples, with the temperature characteristic quantity of the highest back temperature of the thermal protection system during the entire flight process as the output variable; subsequently, based on the generated small number of random samples, construct a surrogate model of the random input sample variables and the output variables, and perform accuracy verification on the surrogate model. Finally, perform rapid calculations on a large number of samples for the surrogate model that meets the accuracy requirements, analyze the probability distribution characteristics of the output features, and perform sensitivity calculation and analysis on the uncertainty parameters. This method can achieve efficient probabilistic design of the ablative thermal protection system during the entire flight process, providing a direction for the design and optimization of ablative thermal protection.

Description

Technical Field

[0001] The present invention belongs to the field related to hypersonic thermal protection design, and particularly relates to a design method applied to an ablative thermal protection system. Background Art

[0002] During the return reentry of high-speed aircraft and their flight in the atmosphere, they face diverse aerodynamic heating environments. As one of the most important subsystems of high-speed aircraft, the thermal protection system plays an irreplaceable role in ensuring the safety of personnel and equipment and the structural integrity of the aircraft. To cope with different thermal environment characteristics, currently, two major types of thermal protection systems have been developed: ablative and non-ablative. The material used in the ablative thermal protection system is a type of pyrolytic carbon material with low density, which undergoes a pyrolysis reaction at high temperatures to absorb heat and generate pyrolysis gases. Its heat protection mechanism mainly includes thermal radiation, heat absorption by heat capacity, heat absorption by pyrolysis reaction, and thermal blocking effect. It can reduce heat conduction by heat absorption through the pyrolysis reaction of the material, release of pyrolysis gases, and surface ablation recession, thereby protecting the internal structure of the aircraft. The structure and forming process of the ablative thermal protection system are relatively simple, and it has the advantages of high efficiency, light weight, simple process, and low cost. Ablative heat protection is used from short-range ballistic missiles to long-range ballistic missiles, spacecraft, recoverable satellites, and low-lift reentry vehicles, and it can be applied to both long-time and tens of seconds, as well as high aerodynamic heating rates. After decades of development, ablative materials have currently become one of the most mature and widely used thermal protection material systems.

[0003] Due to various factors such as ballistic deviation, uncertainty of oncoming flow conditions, dispersion of material properties (especially the thermophysical and thermochemical material properties of ablative thermal protection materials), and deviation of aircraft shape parameters, there are usually varying degrees of uncertainties in the design and use of the thermal protection system, which pose great challenges to the reliability design and evaluation of the system. Due to deficiencies in analysis methods, experimental measurements, and design means, traditional designs usually use safety factors or design margins to consider the influence of uncertainties. However, this method cannot accurately consider the influence of uncertainty distributions on the thermal protection system and usually leads to a relatively conservative design; in some special cases, it cannot fully consider the influence of uncertainties of sensitive parameters, resulting in a reduction in system reliability or even failure to meet engineering requirements. Summary of the Invention

[0004] The present invention aims to provide a refined and practical design method for the ablative thermal protection system of aircraft based on probability techniques in view of the deficiencies of existing design methods.

[0005] The design method applied to the ablative thermal protection system in the present invention includes the following steps:

[0006] S1. Based on the external dimensions and flight mission design parameters of the given aircraft and the selection of thermal protection materials, preliminarily determine the design parameters of the ablative thermal protection system, including the parameters related to the thermal response of the ablative thermal protection system;

[0007] The parameters related to the thermal response include geometric dimension parameters, material property parameters, flight trajectory parameters, and oncoming flow condition parameters.

[0008] S2. Determine the parameters with uncertainty among the parameters related to the thermal response, and quantitatively characterize the uncertainty of these parameters;

[0009] S3. For N parameters with uncertainty, conduct M random samplings according to the probability distribution law corresponding to the quantitative characterization of the uncertainty of each parameter to obtain N groups of sample sequences, with M values in each group;

[0010] S4. Randomly select one sample value from each of these N groups of sample sequences to obtain a random input sample, and then randomly select one sample value from the remaining values of each sample sequence to obtain a new random input sample, and so on, to obtain M random input samples, and the value of any uncertain parameter only appears once in the M random input samples;

[0011] S5. Establish a parametric analysis model for the aerodynamic heat and thermal response of the ablative thermal protection system;

[0012] Perform meshing on the analysis model, and respectively create a nominal surface mesh for aerodynamic heat calculation and a nominal volume mesh for thermal response calculation according to the design dimensions;

[0013] Obtain the surface mesh / volume mesh corresponding to each group of random input samples from the nominal surface mesh / volume mesh according to the mesh deformation algorithm;

[0014] S6. Based on the oncoming flow condition parameters in the random input samples obtained from step S4 and other oncoming flow condition parameters with certainty, use a finite element analysis tool to perform spatial discretization and analysis calculation on the aerodynamic heat analysis model of the thermal protection system, and substitute the oncoming flow condition parameters at the current moment to obtain M linear differential equations, and solve to obtain M cold wall heat fluxes at the current moment

[0015] Solve to obtain M cold wall heat fluxes q during the flight mission for t moments of the entire flight mission of the given aircraft c(T＝300K) ;

[0016] S7. Based on the material property parameters and geometric dimension parameters in the random input samples obtained from step S4, use a finite element analysis tool to perform spatial discretization and analytical calculation on the thermal response analysis model of the thermal protection system. The thermal load is the hot wall heat flux converted from the cold wall heat flux calculated in step S6 for the same sample. A total of M nonlinear differential equations with time t as the independent variable are obtained;

[0017] S8. Solve the M nonlinear differential equations obtained in step S7 to obtain M temperature field vectors T(x, y, z, t), that is, the temperature results of all finite nodes in the entire space of the thermal protection system at each calculation moment during the time history;

[0018] S9. Extract the maximum temperature in the back temperature region of the thermal protection system during the entire time history from these M temperature field vectors T(x, y, z, t) as the output temperature characteristic quantity corresponding to these M uncertain input samples;

[0019] S10. From the M uncertain input parameter combination samples obtained in step S9 and their corresponding output temperature characteristic quantities, construct and train a surrogate model for fitting the functional relationship between the random input sample variables and the output temperature characteristic quantities, and evaluate the accuracy of the surrogate model. If the accuracy requirement is not met, return and repeat steps S3 - S10 until the accuracy of the obtained model reaches the accuracy requirement;

[0020] S11. Sample the uncertain parameters n times again, and use the surrogate model obtained in step S10 to obtain n target random output temperature characteristic quantities, and then calculate the reliability of the thermal protection system under this design parameter condition;

[0021] S12. End when the reliability of the thermal protection system is greater than or equal to the required value;

[0022] Otherwise, adjust the design parameters in the thermal protection system, and repeat steps S2 - S11 until the reliability of the thermal protection system is greater than or equal to the required value.

[0023] Further, the uncertain parameters in the material property parameters include the density of the original material, the specific heat capacity, thermal conductivity of the original material and the charred material, porosity, gas permeability, and material heat of dissociation;

[0024] The uncertainty quantification characterization of the geometric dimension parameters is obtained according to the measurement of the thermal protection coating by an optical scanner;

[0025] The uncertain parameters in the oncoming flow condition parameters include the altitude, speed, and angle of attack during the entire flight process of the aircraft.

[0026] Preferably, the surrogate model is a Kriging model.

[0027] Further, the sampling method described in steps S3 and S11 is the Monte Carlo sampling method based on Latin hypercube, specifically including: for a parameter with uncertainty, divide its value range into M intervals with the same probability, randomly select a value from these M intervals with the same probability, and the M values obtained form a set of sample sequences of this parameter. A total of N sets of sample sequences of uncertain parameters are obtained.

[0028] Further, in step S8, the Newton iteration method and the implicit time integration method are used to solve the nonlinear differential equations.

[0029] Further, in step S7, the loading method of aerodynamic heat flux at the boundary condition in the analysis and calculation of the thermal response model is;

[0030]

[0031] where C h is the enthalpy convection exchange coefficient, and the heat flux of the hot wall at this moment is obtained by setting the enthalpy convection exchange coefficient and the equivalent gas recovery enthalpy.

[0032] Further, in step S7, the calculation formula for the hot wall heat flux density is;

[0033]

[0034] where q n is the hot wall heat flux, is the heat flux density at the wall temperature of T 0 i.e., the cold wall heat flux, T 0 generally takes a value of 300K, h w is the wall enthalpy, h r is the gas recovery enthalpy, is the enthalpy value of the gas at T 0 at this time.

[0035] Further, the calculation formula for the reliability of the thermal protection system in step S11 is:

[0036]

[0037] where T lim is the limit temperature that the system can withstand, T max is the highest back temperature, N(T lim >T max ) is the number of samples in which the highest back temperature is lower than the limit temperature.

[0038] Further, in step S11, it also includes the step of performing sensitivity analysis on the parameters with uncertainty.

[0039] Further, the sensitivity analysis is calculated using the Sobol index, including generating a sample matrix of n×2N based on sampling, where the first N columns are matrix A and the last N columns are matrix B;

[0040] Construct a matrix For i = 1, 2,..., N, make the i-th column of equal to the i-th column of B, and the remaining columns come from A. The Sobol first-order and total effect index formulas are:

[0041] Var(Y)=Var(f(A)+f(B))

[0042]

[0043]

[0044]

[0045]

[0046] The beneficial effects of the present invention compared with the prior art are as follows:

[0047] This design method

[0048] (1). The present invention proposes a probability design and reliability evaluation method for ablative thermal protection systems, which can effectively reduce the design margin of ablative thermal protection systems, thereby reducing the weight of the thermal protection system.

[0049] (2). The present invention takes into account various uncertainty parameters, including geometric deviations, ballistic deviations, oncoming flow deviations, and uncertainties in the properties of thermal protection materials, and comprehensively and comprehensively designs and analyzes the ablative thermal protection system.

[0050] (3). The present invention uses engineering algorithms to calculate transient aerodynamic heat loads, which can quickly and approximately solve aerodynamic heat, thereby quickly obtaining the heat flux load of the thermal protection system.

[0051] Therefore, the present invention accurately quantifies and characterizes the uncertainty of parameters, reasonably determines the thickness of ablative materials based on the probability design method, maximally reduces the weight of the thermal protection system on the premise of ensuring reliability, improves the efficiency of the thermal protection system, and meets the requirements of lightweight and efficient heat protection of aircraft. Description of the Drawings

[0052] Figure 1 is the analysis flow chart of the probability design method for the ablative thermal protection system of the present invention;

[0053] Figure 2 is the schematic diagram of the finite element mesh for aerodynamic heat calculation in the embodiment of the present invention;

[0054] Figure 3 Schematic diagram of a thermal response two-dimensional axisymmetric finite element mesh in an embodiment of the present invention;

[0055] Figure 4 Material-related property curves in an embodiment of the present invention;

[0056] Figure 5 Material-related property curves in an embodiment of the present invention;

[0057] Figure 6 Material-related property curves in an embodiment of the present invention;

[0058] Figure 7 Material-related property curves in an embodiment of the present invention;

[0059] Figure 8 Material-related property curves in an embodiment of the present invention;

[0060] Figure 9 Comparison chart of the calculation results of the surrogate model and the deterministic calculation results in an embodiment of the present invention;

[0061] Figure 10 Probability characteristic distribution and cumulative probability distribution curve graph in an embodiment of the present invention. Detailed implementation manners

[0062] The present invention will be described in detail below in conjunction with the accompanying drawings and specific implementation manners.

[0063] The thermal protection system of a hypersonic vehicle is located on the outermost side. The materials of the ablative thermal protection system are mainly divided into silicon-based, carbon-based, and pyrolytic carbonization types. In recent years, silicone rubber and other pyrolytic carbonization materials have been widely used during the reentry of a spacecraft return capsule and a warhead into the atmosphere. Its main heat dissipation mechanisms are thermal radiation, heat absorption by heat capacity, endothermic pyrolysis reaction, and thermal choking effect, providing thermal insulation protection for its internal structure.

[0064] As Figure 1 shown, the design method applied to the ablative thermal protection system in this embodiment includes the following steps:

[0065] S1. According to the shape and size of the vehicle, the design parameters of the flight process, and the selection of thermal protection materials, preliminarily determine the design parameters of the ablative thermal protection system, including the thermal response-related parameters of the ablative thermal protection system; including geometric dimension output, flight trajectory parameters, oncoming flow condition parameters, material property parameters, etc.

[0066] S2. Based on the relevant parameters of the above thermal protection system, determine the parameters with uncertainty therein and perform uncertainty quantification and characterization on the corresponding parameters. Among the material property parameters, the parameters with uncertainty mainly include the density of the original material, the specific heat capacity (enthalpy value), thermal conductivity, porosity, gas permeability, and material heat of decomposition of the original material and the charred material. The uncertainty quantification and characterization can be carried out by experiments or relevant literature.

[0067] The uncertainty characterization of the geometric dimension parameters is obtained from the measurement of the thermal protection coating by an optical scanner and is subjected to uncertainty quantification and characterization. The oncoming flow condition parameters are the altitude, speed, angle of attack, etc. during the entire flight process of the aircraft. In this embodiment, the uncertainty quantification and characterization of these parameters are carried out in combination with relevant literature.

[0068] S3. Perform random sampling according to the probability distribution law of each parameter with uncertainty characterized above. Using Latin hypercube sampling, for N parameters, divide the value range of each parameter into M intervals with the same probability, and randomly select a value from these M intervals with the same probability. The M values obtained form a set of sample sequences corresponding to the uncertainty parameters, and a total of N sets of sample sequences of uncertainty parameters are obtained.

[0069] S4. Randomly select a sample value from each of these N sets of sample sequences to obtain a random input sample, and then randomly select a sample value from the remaining values of each sample sequence to obtain a new random input sample. By analogy, M random input samples can be obtained and the value of any parameter appears only once in the M random input samples.

[0070] S5. Establish a parametric model of the aerodynamic heat and thermal response of the ablative thermal protection system based on the random input samples obtained in step S4 and other parameters with certainty. First, perform meshing on the aerodynamic heat and thermal response analysis models of the thermal protection system. According to the design dimensions, surface meshes for aerodynamic heat calculation and volume meshes for thermal response calculation as shown in Figure 2 and Figure 3 are respectively made. Then, according to the geometric dimension parameters in each random input sample and the mesh deformation algorithm, the finite element meshes of each set of random input samples are obtained.

[0071] S6. According to the flight trajectory, oncoming flow conditions, and the surface mesh of the aircraft, use a finite element analysis tool to perform aerodynamic heat analysis of the thermal protection system.

[0072] First, use Newton's correction theory to calculate the surface pressure distribution of a hypersonic vehicle. The surface of the vehicle model is meshed using unstructured grids. In each triangular element, the modified Newton theory is used. It is assumed that the momentum component of the oncoming flow in the normal direction of the element is completely lost. According to Newton's second law of motion, the aerodynamic pressure within the element is solved, and finally, the aerodynamic pressure distribution on the surface of the vehicle model is obtained by superposition. Based on the Newton steepest descent method, the streamline trajectory on the surface of the hypersonic vehicle is calculated. Using the previously calculated pressure distribution, the reference enthalpy method, and the fitting formula of gas parameters considering the high-temperature gas effect, the gas parameters under the reference enthalpy conditions of the surface at each data point on the streamline trajectory are calculated. The Reynolds number at the target point is obtained by integrating along the streamline, and then the relevant parameters are substituted into the calculation formula to obtain the heat flux density at the target point.

[0073] In this embodiment, the specific process is as follows:

[0074] S6.1 Mesh the model under study using unstructured grids, output the coordinates of all nodes, and output the node numbers of the three nodes of each triangular element in a certain order.

[0075] S6.2 Given input parameters such as the oncoming flow conditions to meet the input requirements of the aerodynamic calculation program, then read in the subroutine for calculating aerodynamic forces and perform calculations using the modified Newton theory. The calculation of the aerodynamic calculation program is based on elements, and the calculation output result is the aerodynamic pressure at each node.

[0076] S6.3 Use the node and element information output in S6.1 and the aerodynamic pressure calculated in S6.2 as input items and substitute them into the heat flux density calculation program.

[0077] The calculation process of the heat flux density calculation program is as follows:

[0078] S6.3.1 Calculate the streamline trajectory passing through the target point on the surface of the vehicle;

[0079] S6.3.2 Calculate the enthalpy behind the shock wave at the stagnation point using the normal shock relation. Determine the stagnation point position and the pressure at the stagnation point according to the streamline trajectory calculated in the previous step, and determine the entropy value at the stagnation point;

[0080] S6.3.3 Using the isentropic assumption, and combining with the high-temperature gas effect, calculate the enthalpy value distribution at each discrete point on the streamline trajectory using the previously obtained pressure distribution on the surface of the vehicle;

[0081] S6.3.4 Calculate the gas flow velocity, wall enthalpy, adiabatic wall enthalpy, and reference enthalpy at each discrete point on the streamline trajectory according to the known conditions and the results obtained;

[0082] S6.3.5 Calculate the gas density, temperature, sound speed, and viscosity coefficient at each discrete point on the streamline trajectory using the reference enthalpy and pressure distribution;

[0083] S6.3.6 Calculate the Reynolds number at the target point according to the previous calculation results by the following formula:

[0084]

[0085] where x 0 represents the stagnation point, and x f represents the target point, that is, integrate along the calculated streamline from the stagnation point to the target point;

[0086] S6.3.7 Substitute the Reynolds number into the calculation formula to obtain the heat flux density at the target point.

[0087] It should be noted that the heat flux calculated in this way is the cold-wall heat flux. In the part of the aerodynamic heat load in the subsequent finite element calculation, it needs to be converted into the hot-wall heat flux, and the influence of radiation is considered to conduct the thermal response analysis.

[0088] The solution of the aerodynamic heat load uses the steady state at each moment to approximately calculate the transient aerodynamic heat. The specific analysis process is as follows:

[0089] ① Input the ballistic data;

[0090] ② Read the altitude, velocity, and angle of attack at the current time step of the trajectory into the aerodynamic heat calculation program, calculate the overall aerodynamic heat of the model in the steady state, and output the calculation results of the current time step;

[0091] ③ Advance to the next time step and repeat process ②;

[0092] ④ After the calculation of the entire trajectory time step is completed, interpolate the heat flux data in the time dimension to obtain the overall cold-wall heat flux density matrix q c(T＝300K) .

[0093] S7. From the random input samples obtained in the above steps, based on the material property parameters and geometric dimension parameters among them, use the finite element analysis software (SAMCEF-Thermal / Amaryllis) to conduct a thermal response calculation and analysis on the thermal protection system. The required thermal load is approximated by converting the cold-wall heat flux under the same sample calculated in step S6 into the hot-wall heat flux, and a total of M nonlinear partial differential equations with the independent variable of time t are obtained.

[0094] Taking the thermal-structural transient thermal analysis model of a thermal protection material such as silicone rubber that undergoes pyrolytic carbonization as an example, it mainly involves classical heat conduction, porous medium fluid, material pyrolysis reaction, and surface ablation recession. The corresponding mechanism equations are as follows.

[0095] 1). Classical heat conduction: heat balance equation, Fourier's law

[0096] When the heat balance equation is satisfied, the difference in conductive heat flowing into and out of the control volume must be equal to the source term plus the capacitive heat caused by heating or cooling of the material. The instantaneous temperature distribution T(x, t) is given by the following equation when the solid occupies volume V and satisfies the heat balance equation:

[0097]

[0098] In this equation, ρh is the enthalpy of the material and can also be expressed by the heat capacity, Q is the volume source term, and q is the heat flux vector related to the local temperature gradient by Fourier's law.

[0099]

[0100] where q and is a vector of three components, is a 3×3 matrix.

[0101] For the material parameters in the pyrolysis state, such as the thermal conductivity, it can be considered to vary linearly from the original state to the carbonized state as shown in the following equation:

[0102] λ = λ v - α(λ v - λ c ) (4)

[0103] where λ v , λ c are the thermal conductivities in the original state and the carbonized state respectively, and α is the generalized density defined as follows:

[0104]

[0105] 2) Porous medium fluid: Darcy's law

[0106] The relationship between the mass flow rate and pressure of the pyrolysis gas is given by the following equation:

[0107]

[0108] The influence of the gas flow term on the heat equation:

[0109]

[0110] where K P is the gas diffusion coefficient, is the gas mass flow rate, and h g is the gas enthalpy. The gas diffusion coefficient can be directly defined or calculated by the following equation:

[0111]

[0112] M gis the molecular weight of the gas, β is the gas permeability, μ is the gas dynamic viscosity, and R is the gas constant.

[0113] 3), Material pyrolysis reaction: Arrhenius theorem

[0114] The pyrolytic carbonized material will undergo thermochemical decomposition at high temperature to absorb heat and generate pyrolysis gas. The pyrolysis kinetic equation generally adopts the Arrhenius equation. The pyrolysis reaction rate of a single Arrhenius law is as follows:

[0115]

[0116] Where ρ is the density in the current state, and the subscripts c and v represent the "carbonized" and "original" states respectively. A is the reaction frequency (pre-exponential factor), E is the activation energy, N is the reaction order, and R is the gas constant. The influence of material pyrolysis absorption on the heat balance equation is as follows:

[0117]

[0118] In the formula is the heat of pyrolysis of the material.

[0119] The release of pyrolysis gas will produce a heat blocking effect, further reducing the heat flux loading. The reduced heat flux is shown in the following formula:

[0120]

[0121] In the formula, η 2 is the heat blocking coefficient, H a is the recovery enthalpy of the heat flux loading gas, and H w is the gas enthalpy value at the wall temperature.

[0122] 4), Surface ablation recession

[0123] Ablation is the process of removing material from the surface of a structure due to the thermal and / or mechanical loads it bears. Since the ablation phenomenon will cause the removal of material, this will require re-meshing of the structural grid. It is mainly divided into three different types of surface ablation, namely mechanical ablation, chemical ablation, and phase change ablation.

[0124] ① Mechanical ablation

[0125] In this case, the material is removed due to mechanical friction with the surrounding "air" and does not absorb and carry away heat. The surface ablation rate generally depends on temperature, pressure, and shear stress as shown in the following formula:

[0126]

[0127] Where a, b, and T E are constants, and τ, P, and T are shear stress, pressure, and temperature.

[0128] ② Chemical ablation

[0129] It refers to the phenomena such as oxidation, sublimation, and the reaction between carbon and nitrogen on the surface of the material under high temperature and high pressure conditions, resulting in mass loss of the surface material, thereby taking away a large amount of heat. At the same time, chemical ablation will further produce a blocking effect to reduce the heat loading. The ablation rate generally depends on the temperature T, pressure P, and the mass flow rate of the pyrolysis gas released at the boundary.

[0130] ③ Phase change ablation

[0131] When the surface material reaches the phase change temperature T ph , the material changes from the solid phase to the liquid phase, absorbing a large amount of heat and being removed from the surface at the same time. If the temperature is lower than the phase change temperature, the ablation rate is equal to zero. If the temperature reaches the phase change temperature T ph , the ablation rate that satisfies the boundary heat balance can be calculated.

[0132] The aerodynamic heat load of the thermal response analysis model adopts the following cold and hot wall conversion formula:

[0133]

[0134] In the formula, q n is the heat flux of the hot wall, is the heat flux density when the wall temperature is T 0 , that is, the cold wall heat flux. T 0 generally takes a value of 300K. h w is the wall enthalpy, and h r is the gas recovery enthalpy. is the enthalpy value of the gas at T 0 .

[0135] The aerodynamic heat flux loading method at the boundary conditions in SAMCEF-Thermal / Amaryllis software is obtained through the following conversion;

[0136]

[0137] In the formula, C h is the enthalpy convection exchange coefficient, that is, the heat flux density is equivalently set to the recovery enthalpy of the aircraft surface and the enthalpy convection exchange coefficient.

[0138] S8. Use the Newton iteration method and the time implicit integration method to solve the M nonlinear partial differential equations obtained in step S7 to obtain M temperature field vectors T(x, y, z, t), that is, the temperature results of all finite nodes in the entire thermal protection system space at each calculation moment during the time history. The specific solution method is as follows:

[0139] The discretized equation to be solved is as follows:

[0140]

[0141] where C(q) is the heat capacity matrix of the discretized system, K(q) is the heat conduction matrix, and q is the unknown temperature / pressure vector; g e×t (q) is the nodal load vector.

[0142] The heat transfer equation is solved using the first-order implicit time integration method, in which the rate of change of the nodal variables is considered constant over a time step.

[0143] The solution scheme proceeds as follows;

[0144] The time domain is divided into consecutive time intervals between n and n+1, where the individual time:

[0145] t γ =(1 - γ)t n +γt n+1 (16)

[0146] The solution vector can be expressed in the same way, resulting in the following expression:

[0147] q γ =(1 - γ)q n +γq n+1 (17)

[0148] The rate of change of the nodal variables is:

[0149]

[0150] Substituting the above expression into the heat transfer equation gives the following expression:

[0151]

[0152]

[0153] The Newton-Raphson method is used to solve it, and the calculation is carried out as follows:

[0154] i+1 q γ = i q γ + i Δq γ (21)

[0155] The calculation method for the correction is as follows:

[0156]

[0157] The residual value is defined as:

[0158]

[0159] The iterative calculation reaches convergence, and the following expression is used to calculate the solution at time n+1.

[0160]

[0161] S9. Extract the highest temperature in the back temperature region of the thermal protection system during the entire time history from these M temperature field vectors T(x, y, z, t) as the output temperature characteristic quantity corresponding to these M random input samples of uncertainty.

[0162] S10. Based on the M random input samples obtained in step S9 and their corresponding output temperature characteristic quantities, construct and train a surrogate model for fitting the functional relationship between the random input samples and the output temperature characteristic quantities. In this example, the Kriging model is used but not limited to it.

[0163] The Kriging model in this embodiment is composed of a global model and a local deviation superimposed, and can be expressed as

[0164] Y(x) = f(x) + Z(x) (25)

[0165] In the formula, Y(x) is an unknown approximate model, f(x) is a known polynomial function, and Z(x) is a random process with a mean of zero and a non-zero covariance of σ 2 , and the covariance is non-zero. f(x) provides a global approximate model of the design space. Generally, it can be taken as a constant β, while Z(x) creates a local deviation on the basis of the global model. The covariance matrix of Z(x) can be expressed as

[0166] Cov[Z(x i ), Z(x j )] = σ 2 R[R(x i , x j )] (26)

[0167] In the formula, R is the correlation matrix, and R(x i , x j ) represents the correlation function of the sample points x i , x j . There are different forms of the correlation function. Here, the Gaussian correlation function is selected, that is

[0168]

[0169] In the formula, n dv is the number of design variables, and θ k is an unknown correlation parameter. Once the correlation function is determined, the response estimate of any test point x is

[0170]

[0171] Here, y is a column vector of length n s (number of sampling points), that is, the response value of the sample data. When f(x) is a constant, f is a unit column vector of length n s . r T (x) is the correlation vector of length n between the test point x and the sampling points i.e., s

[0172]

[0173] In the formula is estimated by the following formula

[0174]

[0175] The variance estimate of the global model is

[0176]

[0177] The relevant parameter θ is determined by maximum likelihood estimation k , that is, to solve the following non-linear unconstrained optimization problem

[0178]

[0179] When θ k is obtained, the correlation vector r between the unknown point x and the known sample data is obtained from Equation (29) T (x), and its response value is obtained through Equation (28) to complete the construction of the Kriging approximation model.

[0180] The accuracy evaluation formula of the surrogate model is as follows:[[]]

[0181]

[0182] In the formula ntest is the number of test samples, y i is the test sample value, is the response value of the surrogate model, is the mean value. When R 2 is closer to 1, it indicates that the fitting accuracy of the surrogate model is higher. In some other embodiments, other models can also be used and trained, but not limited to, to simulate the relationship between random input samples and output temperature characteristic quantities, such as: polynomial regression model, artificial neural network, and Gaussian process regression, etc.

[0183] ​The M groups of randomly input samples can be divided into training samples and test samples. If the accuracy of the surrogate model fitting is measured to not meet the accuracy requirements (lower than the set threshold value), then it is necessary to return and repeat step S3-10 until the accuracy meets the requirements;

[0184] S11. Perform n times of Monte Carlo simulations based on Latin hypercube sampling using the Kriging surrogate model constructed in step S10 to obtain n target random output temperature characteristic quantities, namely the highest back temperature, for calculating the reliability of the thermal protection system under this design parameter condition; the calculation formula for the reliability of the thermal protection system is:

[0185]

[0186] In the formula, T lim is the limit temperature that the system can withstand, T max is the highest back temperature, and N(T lim >T max ) is the number of cases where the temperature in the highest back result obtained from the random input samples is lower than the limit temperature.

[0187] Through n times of Monte Carlo simulations based on Latin hypercube sampling using the Kriging surrogate model, the probability distribution characteristics of the highest back temperature of the thermal protection system can also be calculated, and its mean and standard deviation can be estimated as:

[0188]

[0189]

[0190] Perform uncertainty parameter sensitivity analysis through the surrogate model, and calculate using the Sobol index. Generate an n×2N sample matrix based on the Monte Carlo sampling method of Latin hypercube. The first N columns are matrix A, and the last N columns are matrix B. Further construct the matrix For i = 1, 2,..., N, make the i-th column of equal to the i-th column of B, and the remaining columns come from A. The Sobol first-order and total effect index formulas are:

[0191] Var(Y) = Var(f(A) + f(B)) (37)

[0192]

[0193]

[0194]

[0195]

[0196] S12. Determine whether the reliability of the thermal protection system is greater than or equal to the required value. If it is, end the process; otherwise, adjust the design parameters in the thermal protection system and repeat steps (2) to (11) until the reliability of the thermal protection system is greater than or equal to the set minimum requirement.

[0197] In this embodiment, taking the design of an ablative thermal protection system for a spherical component as an example, the method provided by the present invention is used for probability analysis and parameter design. The thermal protection structure model adopts a two-dimensional axisymmetric grid as Figure 3 shown, considering the uncertainty parameters of the incoming flow conditions, material properties, and geometric dimensions.

[0198] Among the incoming flow condition parameters, the uncertain parameters are the incoming flow static temperature, density, and velocity. By referring to relevant literature, the dispersion of the incoming flow condition parameters can reach 20%, so the standard deviation of the incoming flow condition parameters is set to 0.05. The uncertainty quantification of the material property parameters through experiments and relevant literature is shown in Table 1. The geometric deviations of the spherical head model mainly consider the matrix outer radius and the thickness of the coating material as shown in Table 2.

[0199] Table 1 Material property parameters and their uncertainty quantification

[0200]

[0201] Table 2 Geometric dimension parameters and their uncertainty quantification

[0202]

[0203] Use the method in the embodiment for random sampling simulation. Conduct 100 random samplings to obtain 100 groups of random uncertainty parameter and maximum back temperature sample data. Use 50 groups as the training samples to construct a surrogate model, and the remaining 50 groups as the test samples. Evaluate the accuracy of the surrogate model, R 2 = 0.9982, which is higher than the set threshold value of 0.99. Then, use the surrogate model to conduct 10 6 random sampling calculations. Conduct probability characteristic analysis on the results. The mean value is 141.1581 °C and the standard deviation is 25.6475 °C. When calculated according to the above method and T lim = 230 °C, the reliability index of the system is 99.93%. Use the above sensitivity calculation method to sample and generate an n×(2×N) sample matrix, where n is 10 6 and N = 12 is the total number of uncertain parameters. Calculate the Sobol first-order and total effect indices as shown below.

[0204] Table 3 Sobol indices of the uncertain parameters in the spherical head model

[0205]

Claims

1. A design method for an ablative thermal protection system, characterized in that, it includes the following steps: S1. According to the external dimension and flight profile design parameters of the given aircraft and the selection of thermal protection materials, preliminarily determine the design parameters of the ablative thermal protection system, including the thermal response related parameters of the ablative thermal protection system; The thermal response related parameters include geometric dimension parameters, material property parameters, flight trajectory parameters and oncoming flow condition parameters; S2. Determine the parameters with uncertainty among the thermal response related parameters and quantify and characterize these uncertainties; S3. For N parameters with uncertainty, conduct M random samplings according to the probability distribution law corresponding to the uncertainty quantification characterization of each parameter to obtain N groups of sample sequences, each group having M values; S4. Randomly select one sample value from each of these N groups of sample sequences to obtain a random input sample, and then randomly select one sample value from the remaining values of each sample sequence to obtain a new random input sample, and so on, M random input samples can be obtained and the value of any uncertainty parameter only appears once in the M random input samples; S5. Establish a parametric analysis model for the aerodynamic heat and thermal response of the ablative thermal protection system; Perform grid processing on the analysis model, and respectively create a nominal surface grid for aerodynamic heat calculation and a nominal volume grid for thermal response calculation according to the design dimensions; Obtain the surface grid / volume grid corresponding to each group of random input samples from the nominal surface grid / volume grid according to the grid deformation algorithm; S6. Based on the oncoming flow condition parameters in the random input samples obtained in step S4 and other deterministic oncoming flow condition parameters, use a finite element analysis tool to perform spatial discretization and analytical calculations on the aerodynamic heat analysis model of the thermal protection system, substitute the oncoming flow condition parameters at the current moment to obtain M linear differential equations, and solve to obtain the cold wall heat fluxes at the current moment The cold-wall heat flux q at M flight histories is solved for t moments during the entire flight history of the given aircraft c(T＝300K) ; S7. Based on the material property parameters and geometric dimension parameters in the random input samples obtained from step S4, use a finite element analysis tool to perform spatial discretization and analysis calculation on the thermal response analysis model of the thermal protection system. The thermal load is converted into the hot wall heat flux according to the cold wall heat flux calculated for the same sample in step S6, and a total of M nonlinear differential equations with time t as the independent variable are obtained; S8. Solve the M nonlinear differential equations obtained in step S7 to obtain M temperature field vectors T(x, y, z, t), that is, the temperature results of all finite nodes in the entire space of the thermal protection system at each calculation moment during the time history; S9. Extract the highest temperature in the back temperature region of the thermal protection system during the entire time history from these M temperature field vectors T(x, y, z, t) as the output temperature characteristic quantity corresponding to these M uncertain input samples; S10. From the M groups of uncertain input parameter combination samples and their corresponding output temperature characteristic quantities obtained in step S9, construct and train a surrogate model for fitting the functional relationship between the random input sample variables and the output temperature characteristic quantities, and evaluate the accuracy of the surrogate model. If the accuracy requirement is not met, return and repeat steps S3 - S10 until the accuracy of the obtained model reaches the accuracy requirement; S11. Conduct n samplings on each parameter with uncertainty again, and use the surrogate model obtained in step S10 to obtain n target random output temperature characteristic quantities, and then calculate the reliability of the thermal protection system under this design parameter condition; S12. End when the reliability of the thermal protection system is greater than or equal to the required value; Otherwise, adjust the design parameters in the thermal protection system and repeat steps S2 - S11 until the reliability of the thermal protection system is greater than or equal to the required value.

2. The method according to claim 1, wherein, the surrogate model is a Kriging model.

3. The method according to claim 1, wherein, the parameters with uncertainty in the material property parameters include the density of the raw material, the specific heat capacity and thermal conductivity of the raw material and the carbonized material, the porosity, the gas permeability, and the heat of pyrolysis of the material; the quantification and characterization of the uncertainty of the geometric dimension parameters are obtained from the measurement of the thermal protection coating by an optical scanner; the parameters with uncertainty in the oncoming flow condition parameters include the altitude, speed, and angle of attack during the entire flight of the aircraft.

4. The method according to claim 1, wherein, the sampling method in steps S3 and S11 is the Monte Carlo sampling method based on Latin hypercube, specifically including: for a parameter with uncertainty, divide its value range into M intervals with the same probability, randomly select a value from these M intervals with the same probability, and the M values obtained form a set of sample sequences of this parameter. A total of N sets of sample sequences of the uncertain parameters are obtained.

5. The method according to claim 1, wherein, in step S8, the Newton - Raphson iteration method and the implicit time - integration method are used to solve the non - linear differential equations.

6. The method according to claim 1, wherein, in step S7, the loading method of the aerodynamic heat flux at the boundary condition in the analysis and calculation of the thermal response model is; Where C h is the enthalpy convection exchange coefficient. By setting the enthalpy convection exchange coefficient and the gas recovery enthalpy equivalently, the heat flux of the hot wall at this moment, q n is the heat flux of the hot wall, is the heat flux density when the wall temperature is T 0 , that is, the cold wall heat flux, T 0 takes the value of 300 K, h w is the wall enthalpy, h r is the gas recovery enthalpy, is the enthalpy value of the gas at T 0 .

7. The method according to claim 1, wherein, in step S7, the calculation formula for the heat flux density of the hot wall is; where q n is the hot-wall heat flux, is the heat flux density at the wall temperature of T 0 , i.e., the cold-wall heat flux, T 0 is taken as 300 K, h w is the wall enthalpy, h r is the gas recovery enthalpy, is the enthalpy value of the gas at T 0 .

8. The method according to claim 1, wherein, the calculation formula for the reliability of the thermal protection system in step S11 is: where T lim is the maximum temperature that the system can withstand, T max is the highest back temperature, N(T lim > T max ) is the number of samples in which the highest back temperature is lower than the maximum temperature.

9. The method according to claim 1, wherein, in step S11, it also includes a step of performing a sensitivity analysis on the parameters with uncertainty.

Images