CN115982854A - Design method applied to ablation type thermal protection system - Google Patents
Design method applied to ablation type thermal protection system Download PDFInfo
- Publication number
- CN115982854A CN115982854A CN202310002424.0A CN202310002424A CN115982854A CN 115982854 A CN115982854 A CN 115982854A CN 202310002424 A CN202310002424 A CN 202310002424A CN 115982854 A CN115982854 A CN 115982854A
- Authority
- CN
- China
- Prior art keywords
- parameters
- protection system
- heat
- uncertainty
- thermal protection
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
Images
Classifications
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Aerodynamic Tests, Hydrodynamic Tests, Wind Tunnels, And Water Tanks (AREA)
Abstract
The invention discloses a design method applied to an ablation type heat protection system. The uncertainty input parameters considered quantize various uncertainty random parameters and generate random input samples by adopting a sampling technology; then, carrying out deterministic pneumatic heat and heat conduction analysis on a random input sample, and taking the temperature characteristic quantity of the highest back temperature of the heat protection system in the whole flight process as an output variable; and then constructing a proxy model of random input sample variables and output variables on the basis of a small amount of generated random samples, and carrying out precision check on the proxy model. And finally, carrying out rapid calculation on a large number of samples of the proxy model meeting the precision requirement, analyzing the probability distribution characteristic of the output characteristic and carrying out sensitivity calculation analysis on the uncertain parameters. The method can realize the efficient probability design of the ablation type heat protection system in the whole flight process, and provides a direction for the design and optimization of the ablation type heat protection.
Description
Technical Field
The invention belongs to the field related to hypersonic heat protection design, and particularly relates to a design method applied to an ablation type heat protection system.
Background
The high-speed aircraft faces diversified pneumatic heating environments in the process of returning and re-entering and flying in the atmosphere, and the thermal protection system serving as one of the most important subsystems of the high-speed aircraft has irreplaceable effects in the aspects of ensuring personnel equipment safety and aircraft structural integrity. In order to deal with different thermal environment characteristics, two types of ablation-type and non-ablation-type thermal protection systems are developed, wherein the material used in the ablation-type thermal protection system is a pyrolytic carbon material which has low density and can generate pyrolytic reaction to absorb heat and generate pyrolytic gas at high temperature, and the heat protection mechanism mainly comprises heat radiation, heat capacity heat absorption, pyrolytic reaction heat absorption and thermal blockage effect, and the functions of reducing heat conduction can be achieved through material pyrolytic reaction heat absorption, pyrolytic gas release and surface ablation retreat, so that the internal structure of the aircraft is protected. The ablation type thermal protection system has the advantages of simple structure and forming process, high efficiency, light weight, simple process, low cost and the like, ablation heat protection is adopted from short-range ballistic missiles to long-range ballistic missiles, airships, return satellites and low-lift reentrants, and the ablation type thermal protection system is applicable to long time to dozens of seconds and high pneumatic heating rate. After decades of development, ablative materials have become one of the most mature and widespread thermal protection material systems.
Due to various factors such as ballistic deviation, uncertainty of inflow conditions, dispersibility of material properties (especially thermophysical and thermochemical material properties of ablation-type thermal protection materials), deviation of aircraft appearance parameters and the like, uncertainty of different degrees generally exists in design and use of a thermal protection system, and great challenge is brought to reliability design and evaluation of the system. Because of the defects in the aspects of analysis methods, test measurement, design means and the like, the influence of uncertainty is usually considered by adopting a safety coefficient or design margin in the traditional design, but the influence of uncertainty distribution on a thermal protection system cannot be accurately considered by the method, and the conservative design is usually caused; under some special conditions, the influence of uncertainty of sensitive parameters cannot be fully considered, so that the reliability of the system is reduced and even the engineering requirements cannot be met.
Disclosure of Invention
Aiming at the defects of the existing design method, the invention provides a design method of an aircraft ablation type heat protection system which is refined based on probability technology and accords with the actual situation.
The design method applied to the ablation type heat protection system comprises the following steps:
s1, preliminarily determining all design parameters of an ablation type heat protection system according to the overall dimension and flight history design parameters of a given aircraft and the type selection of a heat protection material, wherein the design parameters comprise heat response related parameters of the ablation type heat protection system;
the thermal response-related parameters include geometric parameters, material property parameters, ballistic parameters, and inflow condition parameters.
S2, determining parameters with uncertainty according to the thermal response related parameters, and carrying out uncertainty quantitative characterization on the parameters;
s3, carrying out random sampling on N parameters with uncertainty for M times according to a probability distribution rule corresponding to the uncertainty quantization characterization of each parameter to obtain N groups of sample sequences with M values in each group;
s4, randomly selecting a sample value from the N groups of sample sequences to obtain a random input sample, randomly selecting a sample value from the rest values of the sample sequences to obtain a new random input sample, and so on to obtain M random input samples, wherein the value of any uncertainty parameter only appears once in the M random input samples;
s5, establishing a parameterized analysis model of aerodynamic heat and thermal response of the ablation type thermal protection system;
carrying out gridding processing on the analysis model, and respectively making a nominal surface grid for aerodynamic heat calculation and a nominal body grid for thermal response calculation according to the design size;
obtaining a face grid/body grid corresponding to each group of random input samples from the nominal face grid/body grid according to a grid deformation algorithm;
s6, according to the inflow condition parameters and other inflow condition parameters with determinacy in the random input sample obtained in the step S4, adopting a finite element analysis tool to carry out space dispersion and analysis calculation on the thermal protection system pneumatic thermal analysis model, substituting the inflow condition parameters at the current moment to obtain M linear differential equation sets, and solving to obtain M cold wall heat flow at the current moment
Solving t moments of the given aircraft in the whole flight process to obtain cold wall heat flow q in M flight processes c(T＝300K) ；
S7, according to the material attribute parameters and the geometric dimension parameters in the random input sample obtained in the step S4, performing space dispersion and analysis calculation on a thermal response analysis model of the thermal protection system by using a finite element analysis tool, wherein the heat load is calculated according to the step S6, the hot wall heat flow converted from the cold wall heat flow in the same sample is calculated, and M nonlinear differential equation sets with time t as an independent variable are obtained;
s8, solving the M nonlinear differential equation sets obtained in the step S7 to obtain M temperature field vectors T (x, y, z, T), namely temperature results of all finite nodes in the whole thermal protection system space at each calculation moment in a time history;
s9, extracting the highest temperature of a back temperature area of the thermal protection system in the whole time history from the M temperature field vectors T (x, y, z, T) as output temperature characteristic quantities corresponding to the M uncertain input samples;
s10, constructing and training a proxy model for fitting a functional relation between random input sample variables and output temperature characteristic quantities by using the M uncertain input parameter combination samples obtained in the step S9 and the corresponding output temperature characteristic quantities, evaluating the precision of the proxy model, and returning to and repeating the steps S3-S10 if the precision requirement is not met until the precision of the obtained model meets the precision requirement;
s11, sampling each parameter with uncertainty n times again, obtaining n target random output temperature characteristic quantities by using the proxy model obtained in the step S10, and further calculating the reliability of the thermal protection system under the condition of the design parameter;
s12, ending when the reliability of the thermal protection system is greater than or equal to a required value;
otherwise, adjusting the design parameters in the thermal protection system, and repeating the steps S2-S11 until the reliability of the thermal protection system is greater than or equal to the required value.
Further, the parameters with uncertainty in the material property parameters comprise the density of the original material, the specific heat capacity, the thermal conductivity, the porosity, the gas permeability and the material heat of the original material and the carbonized material;
the uncertainty quantitative representation of the geometric dimension parameters is obtained according to the measurement of the optical scanner on the thermal protection coating;
the parameters with uncertainty in the inflow condition parameters include altitude, speed, angle of attack throughout the flight history of the aircraft.
Preferably, the proxy model is a Kriging model.
Further, the sampling method in step S3 and step S11 is a monte carlo sampling method based on the latin hypercube, and specifically includes dividing a certain parameter with uncertainty into M intervals with the same probability in the value range, randomly selecting a value from the M intervals with the same probability, forming a group of sample sequences of the parameter by the obtained M values, and obtaining N groups of sample sequences of the uncertainty parameter.
Further, in step S8, a newton iteration method and an implicit time integration method are used to solve the nonlinear differential equation set.
Further, in step S7, the pneumatic heat flow loading manner at the boundary condition in the analysis and calculation of the thermal response model is as follows;
in the formula C h And (3) setting the enthalpy convection exchange coefficient and the gas recovery enthalpy equivalent to obtain the hot wall heat flow at the moment.
Further, in step S7, the hot wall heat flux density calculation formula is as follows;
in the formula q n In order to heat the flow of the hot wall,
is a wall surface temperature of T 0 Heat flux density of time, i.e. cold wall heat flux, T 0 The value is generally 300K w Is the wall enthalpy, h r For restoring enthalpy to the gas>
Is qiIs located in T 0 The enthalpy value of the time.
Further, the reliability of the thermal protection system in step S11 is calculated by the following formula:
in the formula T lim Limit temperature, T, to which the system can withstand max At the maximum back temperature, N (T) lim ＞T max ) The amount by which the maximum back temperature in the sample is below the threshold temperature.
Further, step S11 includes a step of performing sensitivity analysis on the parameter with uncertainty.
Further, the sensitivity analysis adopts Sobol index calculation, which comprises the steps of generating an N multiplied by 2N sample matrix based on sampling, wherein the first N columns are a matrix A, and the last N columns are a matrix B;
constructing matrices
For i =1,2,. N, such that £ is @>
Is equal to the ith column of B, and the remaining columns are from a, with the Sobol first order and total effect index equations as:
Var(Y)＝Var(f(A)+f(B))
compared with the prior art, the invention has the beneficial effects that:
the design method
(1) The invention provides a probability design and reliability evaluation method for an ablation type heat protection system, which can effectively reduce the design margin of the ablation type heat protection system, thereby reducing the weight of the heat protection system.
(2) According to the ablation-type thermal protection system, various uncertain parameters including geometric deviation, trajectory deviation, inflow deviation and uncertainty of the properties of the thermal protection material are considered, and the ablation-type thermal protection system is comprehensively designed and analyzed.
(3) The transient aerodynamic heat load is calculated by using an engineering algorithm, and the aerodynamic heat can be quickly and approximately solved, so that the heat flow load of the thermal protection system can be quickly obtained.
Therefore, the method accurately quantifies and characterizes the uncertainty of the parameters, reasonably determines the thickness of the ablation material based on a probability design method, reduces the weight of the thermal protection system as much as possible on the premise of ensuring the reliability, improves the efficiency of the thermal protection system, and meets the requirements of light weight and high-efficiency heat protection of an aircraft.
Drawings
FIG. 1 is a flowchart illustrating an exemplary embodiment of a method for designing probability of an ablation type thermal protection system according to the present invention;
FIG. 2 is a schematic view of an aerodynamic thermal computational finite element mesh in an embodiment of the present invention;
FIG. 3 is a schematic diagram of a thermally responsive two-dimensional axisymmetric finite element mesh in an embodiment of the present invention;
FIG. 4 is a graph of material dependent properties in an embodiment of the present invention;
FIG. 5 is a graph of material dependent properties in an embodiment of the present invention;
FIG. 6 is a material dependent property curve in an embodiment of the present invention;
FIG. 7 is a graph of material dependent properties in an embodiment of the present invention;
FIG. 8 is a graph of material dependent properties in an embodiment of the present invention;
FIG. 9 is a graph comparing the computation results of the proxy model with the deterministic computation results in an embodiment of the present invention;
FIG. 10 is a graph of a probability distribution versus a cumulative probability distribution in an embodiment of the invention.
Detailed Description
The present invention will be described in detail below with reference to the accompanying drawings and specific embodiments.
The thermal protection system of the hypersonic aerocraft is positioned at the outermost side, and ablation type thermal protection system materials are mainly silicon-based, carbon-based and pyrolysis carbonization type. In recent years, in the process of returning a space capsule and a warhead to the atmosphere, a pyrolytic carbonized material such as silicon rubber is widely used, and the main heat dissipation mechanisms of the pyrolytic carbonized material are heat radiation, heat capacity heat absorption, pyrolytic reaction heat absorption and thermal blocking effect, so that the internal structure of the pyrolytic carbonized material is protected from heat insulation.
As shown in fig. 1, the design method applied to the ablation type heat protection system in the embodiment includes the following steps:
s1, preliminarily determining various design parameters of an ablation type heat protection system according to the overall dimension and flight history design parameters of an aircraft and the type selection of a heat protection material, wherein the design parameters comprise heat response related parameters of the ablation type heat protection system; including geometric yield, ballistic parameters, inflow condition parameters, material property parameters, and the like.
S2, according to the related parameters of the thermal protection system, determining parameters with uncertainty and carrying out uncertainty quantitative characterization on the corresponding parameters, wherein the parameters with uncertainty in the material attribute parameters mainly comprise the density of an original material, the specific heat capacity (enthalpy value), the thermal conductivity, the porosity, the gas permeability, the material heat release and the like of the original material and a carbonized material, and the uncertainty quantitative characterization can be carried out by adopting experiments or related documents;
the uncertainty characterization of the geometric dimension parameters is obtained according to the measurement of the thermal protection coating by the optical scanner, and the uncertainty quantitative characterization is carried out. The parameters of the inflow condition are altitude, speed, angle of attack, etc. in the whole flight history of the aircraft, and in the embodiment, uncertainty quantitative characterization is performed on these parameters in combination with relevant documents.
And S3, randomly sampling according to the represented probability distribution rule of each parameter with uncertainty, adopting Latin hypercube sampling, dividing the value range of each parameter into M intervals with the same probability for N parameters, randomly selecting a value from the M intervals with the same probability, forming a group of sample sequences corresponding to the uncertainty parameters by the obtained M values, and obtaining the sample sequences of the N groups of uncertainty parameters.
And S4, randomly selecting a sample value from the N groups of sample sequences to obtain a random input sample, randomly selecting a sample value from the rest values of the sample sequences to obtain a new random input sample, and so on to obtain M random input samples, wherein the value of any parameter only appears once in the M random input samples.
And S5, establishing parameterized models of aerodynamic heat and thermal response of the ablation type thermal protection system by the random input samples obtained in the step S4 and other parameters with certainty, firstly carrying out gridding processing on the analysis models of the aerodynamic heat and the thermal response of the thermal protection system, respectively making a surface grid for aerodynamic heat calculation and a volume grid for thermal response calculation as shown in figures 2 and 3 according to design dimensions, and then obtaining respective finite element grids of each group of random input samples according to geometric dimension parameters and a grid deformation algorithm in each random input sample.
And S6, performing aerodynamic thermal analysis on the thermal protection system by adopting a finite element analysis tool according to the flight trajectory, the incoming flow condition and the aircraft surface mesh.
Firstly, calculating the surface pressure distribution of the hypersonic aircraft by utilizing a Newton correction theory, dividing the surface of an aircraft model by adopting a non-structural grid, utilizing the modified Newton theory in each triangular surface element, supposing that the momentum component of the incoming flow in the normal direction of the surface element is completely lost, solving the pneumatic pressure in the surface element according to a Newton second motion law, and finally, superposing to obtain the pneumatic pressure distribution on the surface of the aircraft model. Calculating the surface flow line track of the hypersonic aircraft based on a Newton steepest descent method, calculating gas parameters under the condition of surface reference enthalpy of each data point on the flow line track by utilizing the previously calculated pressure distribution and reference enthalpy method and a fitting formula of the gas parameters under the condition of considering the high-temperature gas effect, obtaining the Reynolds number at a target point by integrating along the flow line, and then substituting related parameters into the calculation formula to obtain the heat flow density at the target point.
In this embodiment, the specific process is as follows:
s6.1, carrying out mesh division on the researched model by using a non-structural mesh, outputting coordinates of all nodes and outputting node numbers of three nodes of each triangular surface element according to a certain sequence.
S6.2, input parameters such as incoming flow conditions are given to enable the input parameters to meet the input requirements of the aerodynamic force calculation program, then the aerodynamic force calculation subprogram is read in, and calculation is carried out by adopting a modified Newton theory. The calculation of the aerodynamic force calculation program is performed by taking a surface element as an object, and the calculation output result is the aerodynamic pressure on each node.
And S6.3, taking the node and surface element information output by the S6.1 and the pneumatic pressure obtained by calculation of the S6.2 as input items, and carrying the input items into a heat flow density calculation program.
The calculation flow of the heat flow density calculation program is as follows:
s6.3.1, calculating a streamline track of the surface of the aircraft passing through a target point;
s6.3.2, computing enthalpy after a shock wave at the stagnation point by using a normal shock wave relational expression, determining the stagnation point position and the pressure at the stagnation point according to the streamline track computed in the previous step, and determining an entropy value at the stagnation point;
s6.3.3, calculating enthalpy value distribution at each discrete point on a streamline track by using an isentropic assumption and the obtained surface pressure distribution of the aircraft and combining a high-temperature gas effect;
s6.3.4, calculating the air flow speed, the wall enthalpy, the heat insulation wall enthalpy and the reference enthalpy at each discrete point on the streamline locus according to the known conditions and the result;
s6.3.5, calculating the airflow density, the temperature, the sound velocity and the viscosity coefficient at each discrete point on the streamline track by using the reference enthalpy and the pressure distribution;
s6.3.6, calculating the Reynolds number at the target point according to the previous calculation result by the following formula:
wherein x is 0 Representing stagnation point, x f Representing the target point, i.e. the integration from the stagnation point to the target point along the calculated streamline;
and S6.3.7, taking the Reynolds number into a calculation formula to obtain the heat flow density of the target point.
It should be noted that the heat flow thus calculated is a cold wall heat flow, and in the following finite element calculation, the aerodynamic heat load part needs to be converted into a hot wall heat flow, and the thermal response analysis is performed by considering the influence of radiation.
The solution of the aerodynamic heat load adopts the steady state at each moment to approximately calculate the transient aerodynamic heat, and the specific analysis flow is as follows:
(1) inputting the ballistic data;
(2) reading the height, the speed and the attack angle of the trajectory at the current time step to an aerodynamic heat calculation program, calculating the overall aerodynamic heat of the model under the steady state condition, and outputting the calculation result of the current time step;
(3) advancing to the next time step, and repeating the process (2);
(4) the whole trajectory time step is calculated, and the heat flow data is interpolated on the time dimension to obtain an integral cold wall heat flow density matrix q c(T＝300K) 。
And S7, randomly inputting the samples obtained in the steps, and performing Thermal response calculation analysis on the Thermal protection system by using finite element analysis software (SAMCEF-Thermal/Amaryllis) according to material property parameters and geometric dimension parameters, wherein the required Thermal load is approximated by converting cold wall heat flow under the same sample obtained by calculation in the step S6 into hot wall heat flow, and M nonlinear partial differential equation sets with time t as an independent variable are obtained in total.
Taking a transient thermal analysis model of a thermal structure of a pyrolytic charring type thermal protection material such as silicon rubber as an example, the model mainly relates to classical heat conduction, porous medium fluid, material pyrolysis reaction and surface ablation retreat, and a corresponding mechanism equation is shown as follows.
1) Classical heat conduction: equation of thermal equilibrium, fourier's law
When the thermal equilibrium equation is satisfied, the difference in the conduction heat into and out of the control volume must be equal to the source term plus the capacitive heat caused by the heating or cooling of the material. The instantaneous temperature distribution T (x, T) is given by the following equation when the solid occupies the volume V and satisfies the heat balance equation:
in this equation, ρ h is the enthalpy of the material and can also be expressed by the heat capacity, Q is the volumetric source term, and Q is the heat flux vector associated with the local temperature gradient by fourier's law.
Wherein q and
is a vector of 3 components>
Is a 3 x 3 matrix. />
For the material parameters in the pyrolysed state, such as thermal conductivity, one can consider the linear to charred state from the virgin state as shown in the following equation:
λ＝λ v -α(λ v -λ c ) (4)
in the formula of v 、λ c The thermal conductivity in the raw state and the carbonized state, respectively, and α is a generalized density defined as follows:
2) Porous medium fluid: law of darcy
The relationship between the mass flow rate and the pressure of the pyrolysis gas is given by:
influence of the gas flow term on the thermal equation:
wherein K P Is the gas diffusion coefficient of the gas,
is the gas mass flow rate, h g Is the gas enthalpy. The gas diffusion coefficient can be defined directly or calculated by the following formula:
M g is the molecular weight of the gas, β is the gas permeability, μ is the gas dynamic viscosity, and R is the gas constant.
3) And material pyrolysis reaction: arrhenius theorem
The pyrolysis carbonization material can perform thermal chemical decomposition and desorption to absorb heat and generate pyrolysis gas at high temperature, the pyrolysis kinetic equation generally adopts an arrhenius equation, and the pyrolysis reaction rate of a single arrhenius law is as follows:
where the density ρ is the density in the current state, and subscripts c and v represent the "carbonized" and "pristine" states, respectively. A is the frequency of the reaction (referred to as the pre-factor), E is the activation energy, N is the order of the reaction, and R is the gas constant. The effect of material thermal desorption on the thermal equilibrium equation is as follows:
in the formula
Is used for heat clearing of materials.
The release of pyrolysis gas creates a thermal choke effect that further reduces the loading of heat flow as shown in the following equation:
in the formula eta 2 Is the thermal blockage coefficient, H a Restoring enthalpy to the loaded hot-fluid gas, H w Is the gas enthalpy at the wall temperature.
4) Ablation back off of the surface
Ablation is a process of removal from the surface of the structure due to the thermal and/or mechanical loads to which the material is subjected, as ablation phenomena can lead to the removal of material, which would require structural meshing. There are three main types of surface ablation, namely mechanical ablation, chemical ablation, and phase change ablation.
(1) Mechanical ablation
In this case, the material is removed and does not absorb the carry-away heat due to mechanical friction with the surrounding "air". The surface ablation rate is generally dependent on temperature, pressure and shear stress as shown in the following equation:
wherein a, b and T E Is constant, and τ, P and T are shear stress, pressure and temperature.
(2) Chemical ablation
The surface material is oxidized, sublimated, reacted between carbon and nitrogen and the like under the conditions of high temperature and high pressure to cause the surface material to generateMass loss to carry away a large amount of heat, while chemical ablation further produces a blocking effect to reduce heat loading
The ablation rate is generally dependent on the temperature T, the pressure P and the pyrolysis gas release mass flow->
(3) Phase change ablation
When the surface material reaches the phase transition temperature T ph When the material changes phase from solid to liquid, it absorbs a lot of heat while being removed from the surface. The ablation rate is equal to zero if the temperature is below the phase transition temperature, T, if the temperature reaches the phase transition temperature ph Then the ablation rate that satisfies the boundary thermal equilibrium can be calculated.
The cold-hot wall conversion formula adopted by the pneumatic heat load of the thermal response analysis model is as follows:
in the formula q n In order to heat the flow of the hot wall,
is a wall surface temperature of T 0 Heat flux density of time, i.e. cold wall heat flux, T 0 The value is generally 300k, h w Is the wall enthalpy, h r To restore enthalpy to the gas>
Is a gas at T 0 Enthalpy value of the time.
The pneumatic heat flow loading mode at the boundary condition in SAMCEF-Thermal/Amaryllis software is obtained by the following conversion;
in the formula C h The enthalpy-convection exchange coefficient, that is, the heat flow density equivalent is set as the aircraft surface recovery enthalpy and the enthalpy-convection exchange coefficient.
S8, solving the M nonlinear partial differential equations obtained in the step S7 by adopting a Newton iteration method and a time implicit integration method to obtain M temperature field vectors T (x, y, z, T), namely the temperature result of all finite nodes in the whole thermal protection system space at each calculation time in a time history, wherein the specific solving method is as follows:
the discretized equation solved is:
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 node load vector.
The heat transfer equation is solved in a first order implicit time integration method in which the rate of change of the node variable is considered constant in one time step.
The solution is carried out in the following manner;
the time domain is divided into successive time intervals between n and n +1, where a single time:
t γ ＝(1-γ)t n +γt n+1 (16)
the solution vector can be represented in the same way, resulting in the following expression:
q γ ＝(1-γ)q n +γq n+1 (17)
the rate of change of the node variables is:
substituting the above expression into the heat transfer equation, the following expression can be obtained:
using Newton iteration method to solve, adopting the following method to calculate:
i+1 q γ ＝ i q γ + i Δq γ (21)
the calculation method of the correction is as follows:
the residual value is defined as:
the convergence to arrival is iteratively calculated, using the following expression to calculate the solution at time n + 1.
S9, extracting the highest temperature of the back temperature area of the thermal protection system in the whole time history from the M temperature field vectors T (x, y, z, T) to serve as the output temperature characteristic quantity corresponding to the M uncertain random input samples.
And S10, constructing and training a proxy model for fitting the functional relationship between the random input samples and the output temperature characteristic quantities by using the M random input samples obtained in the step S9 and the corresponding output temperature characteristic quantities, wherein the proxy model is not limited to a Kriging model.
The Kriging model in this embodiment is formed by superimposing a global model and a local bias, and can be expressed as
Y(x)＝f(x)+Z(x) (25)
Wherein Y (x) is unknownF (x) is a known polynomial function, Z (x) is the mean zero and the variance σ 2 A random process with covariance not equal to zero. f (x) provides a global approximation model of the design space, which can be generally taken as a constant β, while Z (x) creates local biases based on the global model. The covariance matrix of Z (x) can be expressed as
Cov[Z(x i ),Z(x j )]＝σ 2 R[R(x i ,x j )] (26)
Wherein R is a correlation matrix, R (x) i ,x j ) Represents a sample point x i 、x j The correlation function of (2). The correlation functions being of different forms, i.e. Gaussian correlation functions being selected
In the formula n dv To design the number of variables, θ k Are unknown relevant parameters. Once the correlation function is determined, the response of any trial point x is estimated as
Where y is a length n s The column vector of (number of sample points), i.e., the response value of the sample data. When f (x) is constant, f is a length n s The unit column vector of (2). r is T (x) Is a test point x and a sampling point
Has a length of n s Correlation vector, i.e.
In the formula
Estimated from
The variance estimate of the global model is
Determination of the relevant parameter θ by maximum likelihood estimation k I.e. solving the nonlinear unconstrained optimization problem
When theta is measured k After the calculation, the correlation vector r between the unknown point x and the known sample data is obtained by equation (29) T (x) And obtaining the response value through the formula (28) to complete the construction of the Kriging approximate model.
The precision evaluation formula of the proxy model is as follows:
in the formula ntest Is the number of test samples, y i Is to test the value of the sample,
is a proxy model response value, <' > is>
Is the mean value when R 2 Closer to 1, the higher the accuracy of the proxy model fit. In other embodiments, other models may be used and trained to model the relationship between random input samples and output temperature characteristics, such as: polynomial regression models, artificial neural networks, gaussian process regression, and the like.
The M groups of random input samples can be divided into training samples and testing samples, if the fitting precision of the measured proxy model does not meet the precision requirement (is lower than a set threshold value), the step S3-10 needs to be returned again until the precision meets the precision requirement;
s11, carrying out n times of Monte Carlo simulation based on Latin hypercube sampling by using the Kriging agent model constructed in the step S10 to obtain n target random output temperature characteristic quantities, namely the highest back temperature, and the characteristic quantities are used for calculating the reliability of the thermal protection system under the condition of the design parameters; the reliability of the thermal protection system is calculated by the following formula:
in the formula T lim Limit temperature, T, to which the system can withstand max At the maximum back temperature, N (T) lim ＞T max ) The amount of temperature below the threshold temperature in order to obtain the highest back result from the random input samples.
The Monte Carlo simulation based on Latin hypercube sampling is carried out for n times through a Kriging agent model, the probability distribution characteristic of the highest back temperature of the thermal protection system can be calculated, and the mean value and the standard deviation can be estimated as follows:
and carrying out uncertainty parameter sensitivity analysis through the proxy model, and calculating by adopting a Sobol index. An N multiplied by 2N sample matrix is generated by a Monte Carlo sampling method based on the Latin hypercube, the first N columns are matrix A, and the last N columns are matrix B. Further constructing the matrix
For i =1,2, · a, N, so that->
Is equal to the ith column of B, with the remaining columns from a. The Sobol first order and total effect index formula is as follows:
Var(Y)＝Var(f(A)+f(B)) (37)
s12, judging whether the reliability of the thermal protection system is larger than or equal to a required numerical value or not, and ending when the reliability of the thermal protection system is larger than or equal to the required numerical value; otherwise, adjusting design parameters in the thermal protection system, and repeating the steps (2) to (11) until the reliability of the thermal protection system is greater than or equal to the set minimum requirement.
In this embodiment, taking an ablation-type thermal protection system of a spherical member as an example, probability analysis and parameter design are performed by using the method provided by the present invention. The thermal protection structure model adopts a two-dimensional axisymmetric grid as shown in fig. 3, and uncertainty parameters of inflow conditions, material properties and geometric dimensions are considered.
The uncertain parameters in the inflow condition parameters are inflow static temperature, density and speed, the dispersity of the inflow condition parameters can reach 20% by referring to relevant literature data, and then the standard deviation of the inflow condition parameters is set to be 0.05. The uncertainty quantitative characterization of the material property parameters by experiments and related literature is shown in table 1. The geometric deviation of the ball head model is mainly considered by the external radius of the substrate and the thickness of the coating material, and is shown in table 2.
TABLE 1 Material Property parameters and their quantitative characterization of uncertainty
TABLE 2 geometric parameters and their quantitative characterization of uncertainty
Random sampling simulation is carried out by using the method in the embodiment, and 100 times of random sampling are carried out to obtain 100 groups of random uncertainty parameters and the highest back temperature sample data. 50 groups are used as training samples to construct a proxy model, the rest 50 groups are used as testing samples, and precision evaluation R is carried out on the proxy model 2 =0.9982, above a set threshold value of 0.99. Then, a proxy model is used 10 6 Performing random sampling calculation, performing probability characteristic analysis on the result, wherein the average value is 141.1581 ℃, the standard deviation is 25.6475 ℃, and calculating according to the method to obtain T lim =230 ℃, the reliability index of the system is 99.93%. Sampling by the sensitivity calculation method to generate a sample matrix of N x (2 × N), wherein N is 10 6 N =12 is the total number of parameters with uncertainty, and Sobol first order and total effect index are calculated as shown below.
TABLE 3 Sobol index of uncertainty parameters in ball head model
Claims (10)
1. A design method applied to an ablation type thermal protection system is characterized by comprising the following steps:
s1, preliminarily determining all design parameters of an ablation type heat protection system according to the overall dimension and flight history design parameters of a given aircraft and the type selection of a heat protection material, wherein the design parameters comprise heat response related parameters of the ablation type heat protection system;
the thermal response-related parameters include geometric parameters, material property parameters, ballistic parameters, and inflow condition parameters.
S2, determining parameters with uncertainty according to the thermal response related parameters, and carrying out uncertainty quantitative characterization on the parameters;
s3, carrying out random sampling on N parameters with uncertainty for M times according to a probability distribution rule corresponding to the uncertainty quantization characterization of each parameter to obtain N groups of sample sequences with M values in each group;
s4, randomly selecting a sample value from the N groups of sample sequences to obtain a random input sample, randomly selecting a sample value from the rest values of the sample sequences to obtain a new random input sample, and repeating the steps to obtain M random input samples, wherein the value of any uncertainty parameter only appears once in the M random input samples;
s5, establishing a parameterized analysis model of aerodynamic heat and thermal response of the ablation type thermal protection system;
performing gridding processing on the analysis model, and respectively making a nominal surface grid for aerodynamic heat calculation and a nominal body grid for thermal response calculation according to the design size;
obtaining a face grid/body grid corresponding to each group of random input samples from the nominal face grid/body grid according to a grid deformation algorithm;
s6, according to the inflow condition parameters and other inflow condition parameters with determinacy in the random input sample obtained in the step S4, adopting a finite element analysis tool to carry out space dispersion and analysis calculation on the thermal protection system pneumatic thermal analysis model, substituting the inflow condition parameters at the current moment to obtain M linear differential equation sets, and solving to obtain M cold wall heat flow at the current moment
Solving t moments of the given aircraft in the whole flight process to obtain cold wall heat flow q in M flight processes c(T＝300K) ；
S7, according to the material attribute parameters and the geometric dimension parameters in the random input samples obtained in the step S4, performing space dispersion and analysis calculation on a thermal response analysis model of the thermal protection system by using a finite element analysis tool, wherein the heat load is calculated according to the step S6, the heat wall heat flow is converted from the cold wall heat flow in the same sample, and M nonlinear differential equation sets with time t as an independent variable are obtained;
s8, solving the M nonlinear differential equation sets obtained in the step S7 to obtain M temperature field vectors T (x, y, z, T), namely temperature results of all finite nodes in the whole thermal protection system space at each calculation moment in a time history;
s9, extracting the highest temperature of a back temperature area of the thermal protection system in the whole time history from the M temperature field vectors T (x, y, z, T) as output temperature characteristic quantities corresponding to the M uncertain input samples;
s10, constructing and training a proxy model for fitting a functional relation between random input sample variables and output temperature characteristic quantities by using the M uncertain input parameter combination samples obtained in the step S9 and the corresponding output temperature characteristic quantities, evaluating the precision of the proxy model, and returning and repeating the steps S3-S10 if the precision of the proxy model does not meet the precision requirement until the precision of the obtained model meets the precision requirement;
s11, sampling each parameter with uncertainty n times again, obtaining n target random output temperature characteristic quantities by using the proxy model obtained in the step S10, and further calculating the reliability of the thermal protection system under the condition of the design parameter;
s12, ending when the reliability of the thermal protection system is greater than or equal to a required value;
otherwise, adjusting the design parameters in the thermal protection system, and repeating the steps S2-S11 until the reliability of the thermal protection system is greater than or equal to the required value.
2. The method of claim 1, wherein the surrogate model is a Kriging model.
3. The method of claim 1, wherein the parameters of uncertainty in the material property parameters include density of the starting material, specific heat capacity, thermal conductivity, porosity, gas permeability of the starting material and the charred material, and heat of material heat;
the uncertainty quantitative representation of the geometric dimension parameters is obtained according to the measurement of the optical scanner on the thermal protection coating;
the parameters with uncertainty in the inflow condition parameters include altitude, speed, angle of attack throughout the flight history of the aircraft.
4. The method according to claim 1, wherein the sampling method in step S3 and step S11 is a latin hypercube-based monte carlo sampling method, and specifically includes dividing a certain parameter with uncertainty into M intervals with the same probability in the value range, randomly selecting a value from the M intervals with the same probability, and obtaining M values to form a group of sample sequences of the parameter, thereby obtaining N groups of sample sequences of the uncertainty parameter.
5. The method of claim 1, wherein in step S8, the system of nonlinear differential equations is solved by using a newton iteration method and an implicit time integration method.
6. The method according to claim 1, wherein in step S7, the pneumatic heat flow loading manner at the boundary condition in the analytical calculation of the thermal response model is;
in the formula C h For enthalpy convection exchange coefficient, by arranging an enthalpy convection exchange systemThe number and the gas recovery enthalpy are equivalent to obtain the hot wall heat flow at the moment.
7. The method of claim 1, wherein in step S7, the hot wall heat flux density is calculated as;
in the formula q n In order to heat the flow of the hot wall,
is wall temperature of T 0 Heat flux density of time, i.e. cold wall heat flux, T 0 The value is generally 300k, h w Is the wall enthalpy, h r For restoring enthalpy to the gas>
Is a gas at T 0 Enthalpy value of the time.
8. The method according to claim 1, wherein the reliability of the thermal protection system in step S11 is calculated by the following formula:
in the formula T lim Limit temperature, T, to which the system can withstand max At the maximum back temperature, N (T) lim ＞T max ) The amount by which the maximum back temperature in the sample is below the threshold temperature.
9. The method according to claim 1, wherein step S11 further comprises the step of performing sensitivity analysis on the parameter with uncertainty.
10. The method of claim 9, wherein in step S11, the sensitivity analysis is calculated by using Sobol index, and comprises generating an N × 2N sample matrix based on sampling, where the first N columns are matrix a and the last N columns are matrix B;
constructing matrices
For i =1,2,. N, such that £ is @>
Is equal to the ith column of B, with the remaining columns from a. The Sobol first order and total effect index formula is as follows:
Var(Y)＝Var(f(A)+f(B))
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310002424.0A CN115982854A (en) | 2023-01-03 | 2023-01-03 | Design method applied to ablation type thermal protection system |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310002424.0A CN115982854A (en) | 2023-01-03 | 2023-01-03 | Design method applied to ablation type thermal protection system |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN115982854A true CN115982854A (en) | 2023-04-18 |
Family
ID=85959404
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310002424.0A Pending CN115982854A (en) | 2023-01-03 | 2023-01-03 | Design method applied to ablation type thermal protection system |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115982854A (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115494107A (en) * | 2022-10-13 | 2022-12-20 | 重庆大学 | Fine thermophysical property test method for ablation type silicone rubber thermal protection material |
-
2023
- 2023-01-03 CN CN202310002424.0A patent/CN115982854A/en active Pending
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115494107A (en) * | 2022-10-13 | 2022-12-20 | 重庆大学 | Fine thermophysical property test method for ablation type silicone rubber thermal protection material |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Raj et al. | Sensitivity of ice accretion and aerodynamic performance degradation to critical physical and modeling parameters affecting airfoil icing | |
| West IV et al. | Uncertainty quantification of hypersonic reentry flows with sparse sampling and stochastic expansions | |
| Mikielewicz et al. | Temperature, velocity and mean turbulence structure in strongly heated internal gas flows: comparison of numerical predictions with data | |
| Molavi et al. | Inverse identification of thermal properties of charring ablators | |
| Cross et al. | Two-dimensional modeling of ablation and pyrolysis with application to rocket nozzles | |
| Mazzaracchio et al. | A probabilistic sizing tool and Monte Carlo analysis for entry vehicle ablative thermal protection systems | |
| West IV et al. | Uncertainty analysis of radiative heating predictions for Titan entry | |
| Pederson et al. | The Sedov blast wave as a radial piston verification test | |
| Zibitsker et al. | Fully-coupled simulation of low temperature ablator and hypersonic flow solver | |
| CN115982854A (en) | Design method applied to ablation type thermal protection system | |
| Zibitsker et al. | Validation and analysis of a coupled fluid-ablation framework for modeling low-temperature ablator | |
| Tufts et al. | Parabolized stability equation analysis of crossflow instability on HIFiRE-5b flight test | |
| Dulikravich et al. | Inverse problems in aerodynamics, heat transfer, elasticity and materials design | |
| Mohammadiun et al. | Real-time evaluation of severe heat load over moving interface of decomposing composites | |
| Tahmasbi et al. | Inverse identification of temperature-dependent thermal conductivity coefficients in an orthotropic charring composite | |
| Bensayah et al. | Heat transfer in turbulent boundary layers of conical and bell shaped rocket nozzles with complex wall temperature | |
| Hu et al. | Shock wave standoff distance of near space hypersonic vehicles | |
| Lu et al. | Physics-based compressive sensing approach to monitor turbulent flow | |
| Sciacovelli et al. | A priori tests of RANS models for turbulent channel flows of a dense gas | |
| Chen et al. | Aerothermoelastic analysis of a hypersonic vehicle based on thermal modal reconstruction | |
| Liu et al. | Uncertainty analysis of inflow parameters for mid-lift-to-drag vehicle transition | |
| Zettl et al. | Rapid steady-state pressure prediction for ultra high-speed vehicles | |
| Mifsud et al. | A case study on the aerodynamic heating of a hypersonic vehicle | |
| Irick et al. | In-Situ Thermal ROM-Based Optimization Using Borg MOEA: A Preliminary Study | |
| Dai et al. | Heat Flux Prediction of Radiation Balance Wall by Deep Convolutional Neural Networks |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |