# Software Development for Thickness Optimization of Tile-Type Thermal Protection System

Authors: Zhijin Wang; Tao Cheng; Yu Zhou
Journal: International Journal of Aerospace Engineering
 (2022)

---

## Abstract

During the flight mission of hypersonic aircraft, severe aerodynamic heating will occur on the surface, so thermal protection system (TPS) is required to protect the load-bearing structure of the aircraft. The present paper develops an engineering software for automatic optimization of the thickness of tile-type TPS for reusable aircraft. For requirements on TPS of reusable aircraft in the reentry stage, the method of heat flow-time curve enveloping, automatic material selection, and one-dimensional unsteady heat transfer calculation for multilayer plates under thermal load conditions had been researched, an interactive engineering software had been developed. The software improves the calculation accuracy and calculation efficiency of TPS thickness optimization, and it is suitable for rapid design in the conceptual design stage of the aircraft. Finally, by an example, the function of the software is verified.


---

## Body

## 1. Introduction

With the development of aerospace technology, the flight Mach number of the aircraft increases gradually, and the aerodynamic heating of the surface of the aircraft is becoming more and more serious. To ensure flight safety and protect the hypersonic vehicle’s load-bearing structures and sensitive components within acceptable temperature ranges during the entry or reentry flight, TPS needs to withstand high temperatures and keep the internal instrument intact during flight. The selection of the appropriate TPS material is based on the maximum heat flux experienced on the vehicle so that the selected TPS can always protect the internal load-bearing structure in a safe temperature range. The thickness of the selected TPS material depends on the total thermal load throughout the mission to keep the temperature of critical structures within safe limits. Reference [1] reviewed and discussed the structural design and material design of sandwich structures of different TPS types with various configurations, including corrugated cores, foams, and honeycomb cores. The performance of the TPS sandwich structure in terms of the temperature gradient, deformation limit, and mechanical strength is also discussed, and the fabrication method of the TPS sandwich structure for use in hypersonic vehicles is discussed.The ceramic insulation tile is a very typical form of reusable TPS for hypersonic vehicles, as shown in Figure1. Due to the brittleness of its material, the ceramic insulation tile is not suitable to withstand deformation. In addition, the thermal expansion coefficients of the thermal insulation tile and the skin are different, and the difference is obvious when they are deformed by heat. Usually, a strain isolation pad is arranged between the thermal insulation tile and the skin to prevent the thermal insulation tile from being affected by the deformation of the aircraft skin. However, this thermal insulation tile-type TPS still has shortcomings such as weak damage resistance and high maintenance cost. Therefore, in the thermal protection scheme of the thermal insulation tile TPS, the material selection and thickness design of the TPS are the key points and difficulties in the design process.Figure 1
The configuration of a typical ceramic insulation tile [2].In terms of TPS design research, as early as the 1960s, a software named Charring Material Thermal Response & Ablation Program (CMA) was developed by AeroTherm, in the 1990s, NASA’s Ames Research Center developed The Fully Implicit Ablation & Thermal Response Program (FIAT). In the 21st century, the Ames Center developed the Advanced Engineering Environment (AEE) [3] to connect data from different disciplines, which could automate the multidisciplinary design and analysis of spacecraft. And based on AEE, NASA developed MISSION software [4] for orbit design and CBAero software [5] to generate the aerothermal database for aerothermal analysis, which played a huge role in promoting the design of thermal protection of early aircraft.To improve design efficiency and reduce nonoptimal results caused by assumption errors, Cowart and Olds [6] took the lead in proposing the dynamic TPS size design in 1999. By integrating MINIVER, TPSX, and other software, a numerical solution and optimization toolkit for designing TPS thickness, unit weight, and area percentage was developed. McGuire et al. [7] developed the TPSSIZER software, which relies on AEE, to analyze and design the material, thickness, and weight of TPS while reducing the overall weight of the TPS. Chen et al. [8] developed HYAAT software, which includes the database, aerodynamics module, TPS design module, structural optimization, and other submodules, which can solve the entire functional process from heat flow calculation to TPS material selection design and structural optimization. SpaceWorks [9] conducted relevant research on the concept and preliminary design of nonablative TPS for reusable launch vehicles and developed Sentry software for TPS design. However, due to technical limitations at that time, the software generally has the disadvantage of low computational efficiency. In recent years, many scholars have also carried out a lot of research on the design of thermal protection systems [10–12]. Zhang [13] proposed a smoothing algorithm for the inner shape of the TPS curve, which integrated the automatic design process of the TPS and completed the overall task of the conceptual design stage of the TPS. However, the smoothing algorithm Zhang proposed did not distinguish different material areas, resulting in the thickness increment not correlating with the density of the material in each region. Based on the shape of the aircraft, flight orbit parameters, and material data, Xu et al. [14] developed the automatic design method of hypersonic aircraft TPS by using the method of rapid prediction of aerodynamic heat flow through a one-dimensional heat transfer model. Le and Goo [2] designed and analyzed a robust metallic TPS panel made of 304 stainless steel. On the premise of satisfying the strength requirements, the weight is reduced by 60%, and the critical problem of the design of support brackets in the metallic TPS panel was solved, providing a valuable approach for the design of metallic TPS panels in the early stage of metallic TPS development. However, most of the research is still in the stage of numerical method research with independent functions or using Isight [15] to integrate Abaqus (UIIA) or Patran and other finite element analysis software to do optimization. But this integrated optimization method has low calculation efficiency, lacks interactive operation, and has poor user experience, which cannot reach the level of systematic and engineered tile-type TPS thickness automatic optimization design software.This paper focuses on the automatic optimization of the thickness of tile-type TPS for reentry aircraft. Taking the aircraft geometry and heat flux density data as input and aiming at the selection of heat protection materials and thickness optimization of heat protection materials for the whole aircraft, the paper explored the automatic thermal protection thickness design method and developed the tile-type TPS thickness optimization design software (TTTPSTODS). The software contains the complete process of the geometry file input, heat flow data processing, automatic material selection, and tile-type thickness optimization, and so the computational efficiency is improved while minimizing the mass of the TPS system.


## 2. Brief Introduction of Tile-Type TPS Software

By using C++ language [16], combined with OpenGL (Open Graphics Library) [17] graphics interface for programming under the MFC [18] framework, this paper developed a 3D visual, interactive preprocessing, and postprocessing engineering software that can automatically optimize the thickness of tile-type TPS. The engineering software mainly includes the following modules: (1) preprocessing module including 3D model file of analysis object, original heat flow file, and creation of the material library; (2) heat flow-time curve envelope module; (3) automatic material selection module; (4) input of thermal analysis boundary condition; (5) design and optimization for the thickness of tile-type TPS. After the optimization, the thickness distribution and the mass distribution of the tile-type TPS will be automatically generated and demonstrated.The module composition, interrelation, and calculation flow chart of TTTPSTODS are shown in Figure2.Figure 2
Design process of TTTPSTODS.


## 3. The Main Functional Modules in the TTTPSTODS

### 3.1. Heat Flow-Time Curve Envelope Module

The hypersonic vehicle may have multiple return flight tracks, the heat flow-time curve corresponding to each design point of each return track is not the same exactly, and the design of the tile-type TPS needs to be considered for the most severe thermal load case at each design point. Therefore, it is necessary to comprehensively consider each design point to obtain the most extreme heat flow-time curve, including the heat flow-time curves of each track. This problem will be solved in the heat flow-time curve envelope module.For this problem, this paper proposes a new heat flow-time curve envelope method. An example is used to demonstrate the method used in this software system. Two heat flow-time curves (corresponding to two return tracks, respectively) of one design point are randomly selected for analysis. Two heat flow-time curve tracks are shown in Figure3.Figure 3
Two heat flow-time curve tracks.
(a)
Track 1(b)
Track 2The envelope method specifically includes the following three steps:(1)
Unification of Heat Flow-Time Curves of Each Track. The heat flow at each time point corresponding to each return track is different, and the corresponding return time is not necessarily the same. Based on the duration of the heat flow-time curve with the longest time, the time of the other heat flow-time curves is extended to this duration (the processing method at time 0 is similar). For example, in Figure 3, the duration of track 1 is 100 s, and the duration of track 2 is 93 s. Therefore, the time of the heat flow-time curve of track 2 is extended to 100 s, and the heat flow of the extended section can be taken as the heat flow of track 1 in the period of 93-100 s. The delayed heat flow-time curve of track 2 is shown in Figure 4.(2)
The Isochronous Approximation of the Heat Flow-Time Curve of Each Track. Because the set of time points of each original heat flow-time curve is not necessarily the same, for each heat flow-time curve, the time span is set as 10 s (can be set according to the calculation accuracy), and the value of heat flow corresponding to each moment is obtained. If the selected time point is not in the time point set of the original data, linear interpolation is performed by using the heat flow values of the two time points closest to the selected time point, and a curve similar to the original heat flow-time curve is obtained. The isochronous approximation curves of heat flow-time for track 1 and track 2 are shown in Figure 5.(3)
Heat Flow-Time Curve Envelope. For all heat flow-time curves after approximation, at each time point, the corresponding maximum heat flow value is selected in all heat flow-time approximate curves as the envelope heat flow value at this time point, and the heat flow-time envelope curve which contains the value of the most severe heat flow load can be obtained. Figure 6 shows the envelope curve obtained from the heat flow-time curves for track 1 and track 2.Figure 4
Delayed heat flow-time curve of track 2.Figure 5
The isochronous approximation curves of heat flow-time.
(a)
Track 1(b)
Track 2Figure 6
Envelope curve of track 1 and track 2.We compared the heat flow-time curves of track 1 and track 2 and their envelope curve in one picture, so the envelope effect can be seen more clearly, as shown in Figure7.Figure 7
Comparison of original heat flow-time curves and heat flow-time envelope curve.It can be seen from Figure7 that the envelope curve obtained by the heat flow-time curve envelope method can envelop the maximum heat flow value at each time point in most time periods. However, due to the chosen time span being 10 s, the envelope curve misses several maximum heat flow points when the heat flow-time curve changes rapidly. If the time span is reduced, the accuracy of the envelope will be improved, but the calculation time of the envelope will increase significantly. Therefore, it is necessary to weigh the advantages and disadvantages of all parties and to choose a reasonable time span. In order to research the influence of the time span selected in the envelope algorithm on the accuracy, the time span of 3 s, 4 s, 5 s, and 6 s are set, respectively, and the corresponding enveloping curves are given in Figure 8.Figure 8
Heat flow-time enveloping curves corresponding to different time spans.
(a)
Time span 3 s(b)
Time span 4 s(c)
Time span 5 s(d)
Time span 6 sAs can be seen from Figure8, when the time span is 5 s, the enveloping curve misses the server point. As the time span gradually increases, the speed of the calculation will increase, but the accuracy will decrease. In the TPS software developed in this paper, users can choose an appropriate time span for calculation according to their own needs.


### 3.2. Automatic Material Selection Module

During the flight mission of the hypersonic vehicle, each area of the body surface will bear different heat flow loads, and the temperature distribution in each area is also different. Therefore, the upper limit temperature of the material is the main factor in the selection of materials in different areas of the fuselage. The automatic material selection module is realized by the principle and algorithm of automatic selection, which is based on the radiation equilibrium temperature of the surface. The prerequisite of automatic material selection is to obtain the maximum temperature of the TPS surface, then the maximum temperature and the upper limit temperature that the material can bear will be compared. For the most efficient use of materials, the upper limit temperature of the material should be closer to the maximum temperature of the TPS surface as much as possible. To obtain the principle of automatic material selection suitable for tile-type TPS, the boundary conditions of the outer surface of a single tile-type are taken as an example to analyze, as shown in Figure9.Figure 9
Boundary conditions on the outer surface of a single tile-type.In Figure9, λ, ε, and h represent the thermal conductivity coefficient, radiation coefficient, and convection coefficient, respectively, corresponding to the three heat transfer ways of the outer surface boundary of a single block of insulation tile, heat conduction, heat radiation, and heat convection. The relationship between the three heat transfer ways is shown in
(1)qconv=qcond+qrad,where qconv, qcond, and qrad represent the heat flux density of heat convection, the heat flux density of heat conduction, and the heat flux density of heat radiation, respectively. qconv is usually given in the form of heat flow, and obtaining qcond requires calculating the heat transfer inside the structure, the heat flow density of thermal radiation can be expressed as
(2)qrad=εσT4−Tf4,where ε is the emissivity of the object, which is related to the type of material and the surface, σ is the Stefan-Boltzmann constant, the value of which is 5.76×10−8W/m2K4.Bring formula (2) into formula (1) and after simplifying, we will get
(3)T=qconv−qcondεσ+Tf44,where T and Tf represent the outer surface temperature of the structure and ambient temperature, respectively. Formula (3) can accurately calculate the outer surface temperature of the structure.Part of the heat conducted by thermal convection is dissipated to the environment or other parts of the body surface, and the other part enters the interior of the structure through thermal conduction. When the outer surface temperature is high, most of the heat will be dissipated into the atmosphere by thermal radiation, only a small part will enter the structure by thermal conduction. Therefore,qcond can be ignored, and formula (3) can be simplified to:
(4)T=qconvεσ+Tf44.Formula (4) is also called the radiation balance formula. When most of the heat is transferred by thermal radiation, the temperature of the outer surface of the structure can also be estimated by formula (4). During the whole time domain, the maximum temperature Tmax of the outer surface of the structure is obtained when qconv takes the maximum value. Before automatic material selection, the thermal properties of various optional thermal protection materials should be input first, including the material upper limit temperature Tlim. Then, the material whose upper limit temperature is closest to the maximum temperature of the outer surface of the structure will be taken, that is, the target material, as shown in
(5)minTmax−Tlim.In the automatic material selection module of the engineering software, the maximum temperature of the outer surface of the structure at each structural point is calculated according to formula (4), and then the material of each design point will be selected by formula (5). The extreme temperature of the material and other material properties are included in the input material library.


### 3.3. One-Dimensional Heat Transfer Algorithm of Multilayer Plate Based on Controlled Volume Method

Before the calculation and analysis of TPS temperature, a suitable and simplified heat transfer model should be established first to achieve the purpose of ensuring both accuracy and calculation efficiency. Generally, there are three types of thermal conductivity models: one-dimensional model, two-dimensional model, and three-dimensional model. The advantages and disadvantages of the three heat transfer models are shown in Table1.Table 1
Comparison of advantages and disadvantages of three heat transfer models.
One-dimensional heat transfer modelTwo-dimensional heat transfer modelThree-dimensional heat transfer modelAdvantages(1) The model is simple(2) The calculation is efficient and reliable(1) It can simulate temperature changes and heat transfer in both thickness and lateral directions simultaneously(2) The calculation accuracy of the two-dimensional model is higherIt can simulate the heat conduction process in three directions.Disadvantages(1) The lateral temperature change cannot be simulated(2) It affects the accuracy under some complex boundary conditions(1) The calculation speed is slow(2) It is not suitable for the rapid design of TPS in the conceptual design stage(1) The model establishment is complicated(2) The calculation amount is huge(3) The efficiency is extremely lowConsidering the advantages and disadvantages of the three heat transfer models, thermal insulation tiles are focused on the study of heat conduction and temperature distribution along the thickness direction. Additionally, the temperature difference between adjacent tiles is small, radiation effects from the sides of the tiles can be ignored, so a one-dimensional heat transfer model is sufficient.The thermal protection layer may contain multiple layers of different materials, and a multilayer one-dimensional heat transfer analysis calculation should be performed during the thickness optimization of the thermal protection layer. In the heat transfer analysis module, this software adopts the control volume method which is currently used in mainstream commercial software. The method focuses on the integral balance of the control volume, and it used nodes as the representative of the control volume [19]. Firstly, in the control volume method, the variation law of the variable which needs to be determined in the control volume, and the distribution function in the control volume are set. Then, the control quantity is integrated according to the characteristics of the distribution function, and the algebraic formula between the node variables can be obtained. Finally, the formula will be solved and the calculation results of each node will be got. Since there is no gap between the objects, energy must be conserved, so the calculation results obtained by the control body method must satisfy the principle of energy conservation. In other words, the control volume method can be solved for both uniform and nonuniform meshes. Therefore, the control volume size has a wide range and high flexibility without losing physical authenticity.When controlling the balance of the volume integral, the selected variable distribution usually has two forms: stepped distribution and piecewise linear distribution. As shown in Figure10, it is assumed that the x-axis is the one-dimensional transmission coordinate axis of temperature, and the T-axis is the temperature that changes with time τ, that is, the variable. Considering the gradient value of the calculated variable at the control volume interface, the stepped distribution is generally only used for the distribution of source terms, physical parameters, and variables in the time domain. Assuming three adjacent nodes W, P, and E, the control volume interface is represented by a dashed line, the location of which is represented by w and e, and is usually at the midpoint of the two nodes.Figure 10
Two distribution assumptions in the control volume method.For no internal heat source, the one-dimensional unsteady heat conduction control formula, in which the temperature is only conducted along the thickness direction of the multilayer plate, is shown in the following formula(6)ρc∂T∂τ=λ∂2T∂x2,where ∂T/∂τ,∂T/∂x are changing rate of temperature with time and the changing rate of temperature along the one-dimensional thickness direction, and λ and c are the thermal conductivity and specific heat capacity of the material, respectively.As shown in Figure11, for a uniform grid divided by a multilayer plate with no internal heat source and the temperature is only conducted along the thickness direction, the grid length is Δx, and three adjacent nodes n−1,n,n+1 are selected for the research. ΔX represents the control volume P, e, and w are the boundary points of the control volume P. For a uniform grid, e is the midpoint between nodes n and n+1, and w is the midpoint between nodesn−1 and n, so
(7)ΔX=Δx.Figure 11
Mesh nodes and control volumes.Integrating the control volumeP from time τ to τ+Δτ, we can get
(8)∫we∫ττ+Δτρc∂T∂τdτdx=∫ττ+Δτ∫we∂∂xλ∂T∂xdxdτ.When the temperature of the control volumeP is represented by the node temperature, the distribution of the temperature in the time and space coordinates is a step shape, and the integral on the left side of formula (8) can be expressed as
(9)∫we∫ττ+Δτρc∂T∂τdτdx=ρcPTnτ+Δτ−TnτΔx.When the distribution of temperatureT between nodes is piecewise linear, the integral on the right side of formula (8) can be expressed as
(10)∫ττ+Δτ∫we∂∂xλ∂T∂xdxdτ=∫ττ+ΔτλeTn+1−TnΔx−λwTn−Tn−1Δxdτ.Combining formula (9) and formula (10), and we can get the simplified governing formula
(11)ρcPTnτ+Δτ−TnτΔx=∫ττ+ΔτλeTn+1−TnΔx−λwTn−Tn−1Δxdτ.In formula (11), ρcP is the average density and specific heat capacity in the control volume P, and the specific heat capacity varies with the temperature, λe and λw are the equivalent thermal conductivity at points e and w, respectively, and the process of solution is as follows.In Figure12, e represents the upper boundary of the control volume P, Δxe+ is the distance from e to n+1, and Δxe− is the distance from e to n. The heat flow through the interface at point e can be expressed as
(12)q=Te−Tn+1Δxe+/λn+1=Tn−TeΔxe−/λn=Tn−Tn+1Δxe+/λn+1+Δxe−/λn.Figure 12
Control volume interface.In formula (12), q represents the heat flow through the interface at point e, Tn+1, Tn, and Te represent the temperatures at points n+1, n, and e, respectively, λn+1, λn, and λe represent the thermal conductivity at points n+1, n, and e. The heat flux density between nodes n+1 and n can also be expressed by the interfacial thermal conductivity
(13)q=Tn−Tn+1Δx/λe.According to the compatibility principle, the heat flow of formula (12) and formula (13) should be equal, so
(14)Δxλe=Δxe+λn+1+Δxe−λn.In a uniform grid, wheree is the midpoint of n and n+1, we can get
(15)λe=2λn+1λnλn+1+λn.The thermal conductivity of pointe can be obtained by formula (15), and the thermal conductivity of point w can be obtained in the same way. However, although the thermal conductivity of two points is obtained, the integral with respect to time on the right side of formula (11) cannot be obtained, and further approximation is required.Figure13 is a schematic diagram of the weighting factor, and f is the weighting factor. By introducing the weighting factor, formula (11) can be simplified as
(16)ρcPTnτ+Δτ−TnτΔx=fλeTn+1τ+Δτ−Tnτ+ΔτΔx−λwTnτ+Δτ−Tn−1τ+ΔτΔx+1−fλeTn+1τ−TnτΔx−λwTnτ−Tn−1τΔx.Figure 13
Weighting factor of time steps.As shown in Figure13, f is between 0 and 1. f=0 means that the temperature maintains the value of time τ in the interval of Δτ, and does not jump to Tnτ+Δτ until the time is τ+Δτ, when f=1, the temperature jumps to Tnτ+Δτ at the time of τ and remains unchanged throughout the time of Δτ, when f=0.5, the temperature T will be distributed piecewise linearly over time.Whenf=0, the time τ is expressed as time i, and the time τ+Δτ is expressed as time i+1, the formula (16) can be simplified as
(17)Tni+1=ΔτρcPΔx2λeTn+1i−Tni−λwTni−Tn−1i+Tni.In formula (17), the temperature of each node at the next moment i+1 is obtained from the temperature of the three adjacent nodes including the node at the moment i, and there is no need to solve the algebraic formulas simultaneously. Formula (17) is the discretization formula that controls the display format of the volume method.Formula (17) needs to satisfy the principle of positive coefficients, that is, the coefficients of Tni+1 and Tni are the same as positive. Otherwise, physically unreal solutions will appear, resulting in nonconvergence of computational results, so the coefficients Tni in formula (17) must be positive, that is
(18)1−λw+λeΔτρcΔx2>0.Formula (18) can be further simplified and we can get
(19)Δτ<ρcΔx2λw+λe.It can be seen from formula (19) that when Δx is reduced to 0.5 times of the original, the maximum value of Δτ will be reduced to 0.25 times of the original, that is, reduced by the quadratic, which will greatly reduce the calculation efficiency, so it is unrealistic to reduce the value of the grid length to improve the calculation accuracy.Whenf=1, the time τ is expressed as time i, and the time τ+Δτ is expressed as time i+1, formula (16) can be simplified as
(20)λwTn−1i+1−λw+λe+ρcPΔx2ΔτTni+1+λeTn+1i+1=−ρcPΔx2ΔτTni.The characteristic of formula (20) is that the temperature of each node at the next moment i+1 is obtained from the temperature of the moment i+1 including the three adjacent nodes of the node, so each algebraic equation has three unknowns, which cannot be solved directly. It is necessary to list the algebraic equations of all nodes, and solve the algebraic equation system simultaneously, which is often solved by tridiagonal matrix algorithm (TDMA). Formula (20) is called the discretization equation of the fully implicit format of the control volume method. Carefully observe the formula (20), and we can see that the coefficients of Tni+1 and Tni must be both negative, so Δx and Δτ have no constraint relationship, it must have a real solution in the physical sense, that is, the calculation results must converge.TDMA is used to calculate the case where the coefficient matrix is arranged in a tridiagonal arrangement, including two calculation procedures: the elimination of variables to find the coefficient and back-substitution to find the temperature of the node. Firstly, eliminate one variable, so that only two variables remain in the equation, then these two variables are obtained through boundary conditions, all two-variable results are calculated through iteration, and finally, the temperature of each node is calculated by back-substitution again.Suppose the discretization equation is(21)aiTi−1+biTi+ciTi+1=di,i=1,2,⋯,N,(22)b1c1a2b2c2⋯⋯⋯aN−1bN−1cN−1aNbNT1T2⋯TN−1TN=d1d2⋯dN−1dN.The matrix expression of formula (21) is shown in formula (22). The entire matrix only has coefficients around the diagonal and on both sides of the diagonal, forming a tridiagonal matrix.Wheni=1, we have
(23)T1=−c1b1T2+d1b1=P1T2+Q1.Wheni=2, we have
(24)T2=−c2a2P1+b2T3+d2−a2Q1a2P1+b2=P2T3+Q2.The discrete equation for anyi is
(25)Ti=−ciaiPi−1+biTi+1+di−aiQi−1aiPi−1+bi=PiTi+1+Qi.TheNth equation is
(26)aNTN−1+bNTN=dN.SubstituteTN−1=PN−1TN+QN into formula (26) and simplify, we get
(27)TN=dN−aNQN−1aNPN−1+bN=QN.Therefore, the general form of the two-variable equation is(28)Ti=PiTi+1+Qi,i=1,2,⋯,N.In formula (28):
(29)Pi=−ciaiPi−1+bi,Qi=di−aiQi−1aiPi−1+bi.The coefficients of the boundary nodes are as follows, becausea1=0, so
(30)P1=−c1b1,Q1=d1b1.BecausecN=0, so
(31)PN=0,QN=dN−aNQN−1aNPN−1+bN=TN.Therefore, in the positive process, only all the coefficientsPi and Qi need to be calculated, and finally the QN value is obtained. Since tN=QN, all node temperatures can be calculated by applying the general form of the two-variable equation Ti=PiTi+1+Qi in the inverse process back-substitution.It can be seen from the above operation process that in order to obtain allPi and Qi, the conditions of bi≠0 and aiPi−1+bi≠0 must be satisfied. The system of equations has a unique solution only when the conditions are met.Whenf=0.5, the time τ is expressed as time i, and the time τ+Δτ is expressed as time i+1, the formula (16) can be simplified as
(32)λwTn−1i+1−λw+λe+2ρcPΔx2ΔτTni+1+λeTn+1i+1=λwTni−Tn−1i+λeTni−Tn+1i−2ρcPΔx2ΔτTni.Formula (32) is very similar to the combination of formula (17) and formula (20). It needs to solve the algebraic equation system and also needs to satisfy the principle of positive coefficients. Formula (32) is called the discretization equation of the Crank-Nicolson format (C-N format) of the control volume method. The coefficients of Tni+1 and Tni have the same sign (both positive or negative), so
(33)Δτ<2ρcΔx2λw+λe.Comparing formula (19) and formula (33), when the value of Δx is determined, and the discrete control equation in C-N format is used, the value of Δτ is twice the value of the display format. Carefully explore the variation law of temperature change with time, we can see that at small time steps, the C-N format is closer to a linear variation, and the accuracy will be higher at this time. Undoubtedly, taking into account the computational accuracy and physical authenticity, it is ideal to use the exponential change law, but this will increase the difficulty of dealing with specific problems.Using above three discrete methods, we can obtain the governing formulas describing the heat flow inside the structure, but the boundary conditions are also required to provide initial values for the governing formulas, so as to solve the temperature values of all nodes of the structure at different times. Taking the multilayer plate boundary condition in Figure14 as an example, we can use the law of conservation of energy to discretize the boundary conditions.Figure 14
Boundary conditions on the surface of the structure.The boundary conditions on the surface of the structure is shown in Figure14. There are two heat transfer ways: heat convection and heat radiation. The boundary node N represents an element with a width of Δx/2 (the shaded part in Figure 14), and the temperature presents a stepped distribution with time. For the shaded part, the energy is conserved, so we can get
(34)λTN−1i−TNiΔx+hTf−TNi+εσTf4−TN4i=ρcΔx2TNi+1−TNiΔτ.In formula (34), h is the convection coefficient that changes with time, ε is the emissivity of the material on the boundary, which is a fixed value, and Tf is the ambient temperature.Simplifying formula (34) further, we can get
(35)TNi+1=−2εσΔτρcΔxTN4i+1−2λΔτρcΔx2−2hΔτρcΔxTNi+2λΔτρcΔx2TN−1i+2εσΔτρcΔxTf4+2hΔτρcΔxTf.It can be seen from formula (35) that the temperature of the boundary node at the next moment i+1 is obtained from the temperature of three adjacent nodes including the node at the current moment i, and the temperature value of the node at all moments can be calculated by iterative calculation.


## 3.1. Heat Flow-Time Curve Envelope Module

The hypersonic vehicle may have multiple return flight tracks, the heat flow-time curve corresponding to each design point of each return track is not the same exactly, and the design of the tile-type TPS needs to be considered for the most severe thermal load case at each design point. Therefore, it is necessary to comprehensively consider each design point to obtain the most extreme heat flow-time curve, including the heat flow-time curves of each track. This problem will be solved in the heat flow-time curve envelope module.For this problem, this paper proposes a new heat flow-time curve envelope method. An example is used to demonstrate the method used in this software system. Two heat flow-time curves (corresponding to two return tracks, respectively) of one design point are randomly selected for analysis. Two heat flow-time curve tracks are shown in Figure3.Figure 3
Two heat flow-time curve tracks.
(a)
Track 1(b)
Track 2The envelope method specifically includes the following three steps:(1)
Unification of Heat Flow-Time Curves of Each Track. The heat flow at each time point corresponding to each return track is different, and the corresponding return time is not necessarily the same. Based on the duration of the heat flow-time curve with the longest time, the time of the other heat flow-time curves is extended to this duration (the processing method at time 0 is similar). For example, in Figure 3, the duration of track 1 is 100 s, and the duration of track 2 is 93 s. Therefore, the time of the heat flow-time curve of track 2 is extended to 100 s, and the heat flow of the extended section can be taken as the heat flow of track 1 in the period of 93-100 s. The delayed heat flow-time curve of track 2 is shown in Figure 4.(2)
The Isochronous Approximation of the Heat Flow-Time Curve of Each Track. Because the set of time points of each original heat flow-time curve is not necessarily the same, for each heat flow-time curve, the time span is set as 10 s (can be set according to the calculation accuracy), and the value of heat flow corresponding to each moment is obtained. If the selected time point is not in the time point set of the original data, linear interpolation is performed by using the heat flow values of the two time points closest to the selected time point, and a curve similar to the original heat flow-time curve is obtained. The isochronous approximation curves of heat flow-time for track 1 and track 2 are shown in Figure 5.(3)
Heat Flow-Time Curve Envelope. For all heat flow-time curves after approximation, at each time point, the corresponding maximum heat flow value is selected in all heat flow-time approximate curves as the envelope heat flow value at this time point, and the heat flow-time envelope curve which contains the value of the most severe heat flow load can be obtained. Figure 6 shows the envelope curve obtained from the heat flow-time curves for track 1 and track 2.Figure 4
Delayed heat flow-time curve of track 2.Figure 5
The isochronous approximation curves of heat flow-time.
(a)
Track 1(b)
Track 2Figure 6
Envelope curve of track 1 and track 2.We compared the heat flow-time curves of track 1 and track 2 and their envelope curve in one picture, so the envelope effect can be seen more clearly, as shown in Figure7.Figure 7
Comparison of original heat flow-time curves and heat flow-time envelope curve.It can be seen from Figure7 that the envelope curve obtained by the heat flow-time curve envelope method can envelop the maximum heat flow value at each time point in most time periods. However, due to the chosen time span being 10 s, the envelope curve misses several maximum heat flow points when the heat flow-time curve changes rapidly. If the time span is reduced, the accuracy of the envelope will be improved, but the calculation time of the envelope will increase significantly. Therefore, it is necessary to weigh the advantages and disadvantages of all parties and to choose a reasonable time span. In order to research the influence of the time span selected in the envelope algorithm on the accuracy, the time span of 3 s, 4 s, 5 s, and 6 s are set, respectively, and the corresponding enveloping curves are given in Figure 8.Figure 8
Heat flow-time enveloping curves corresponding to different time spans.
(a)
Time span 3 s(b)
Time span 4 s(c)
Time span 5 s(d)
Time span 6 sAs can be seen from Figure8, when the time span is 5 s, the enveloping curve misses the server point. As the time span gradually increases, the speed of the calculation will increase, but the accuracy will decrease. In the TPS software developed in this paper, users can choose an appropriate time span for calculation according to their own needs.


## 3.2. Automatic Material Selection Module

During the flight mission of the hypersonic vehicle, each area of the body surface will bear different heat flow loads, and the temperature distribution in each area is also different. Therefore, the upper limit temperature of the material is the main factor in the selection of materials in different areas of the fuselage. The automatic material selection module is realized by the principle and algorithm of automatic selection, which is based on the radiation equilibrium temperature of the surface. The prerequisite of automatic material selection is to obtain the maximum temperature of the TPS surface, then the maximum temperature and the upper limit temperature that the material can bear will be compared. For the most efficient use of materials, the upper limit temperature of the material should be closer to the maximum temperature of the TPS surface as much as possible. To obtain the principle of automatic material selection suitable for tile-type TPS, the boundary conditions of the outer surface of a single tile-type are taken as an example to analyze, as shown in Figure9.Figure 9
Boundary conditions on the outer surface of a single tile-type.In Figure9, λ, ε, and h represent the thermal conductivity coefficient, radiation coefficient, and convection coefficient, respectively, corresponding to the three heat transfer ways of the outer surface boundary of a single block of insulation tile, heat conduction, heat radiation, and heat convection. The relationship between the three heat transfer ways is shown in
(1)qconv=qcond+qrad,where qconv, qcond, and qrad represent the heat flux density of heat convection, the heat flux density of heat conduction, and the heat flux density of heat radiation, respectively. qconv is usually given in the form of heat flow, and obtaining qcond requires calculating the heat transfer inside the structure, the heat flow density of thermal radiation can be expressed as
(2)qrad=εσT4−Tf4,where ε is the emissivity of the object, which is related to the type of material and the surface, σ is the Stefan-Boltzmann constant, the value of which is 5.76×10−8W/m2K4.Bring formula (2) into formula (1) and after simplifying, we will get
(3)T=qconv−qcondεσ+Tf44,where T and Tf represent the outer surface temperature of the structure and ambient temperature, respectively. Formula (3) can accurately calculate the outer surface temperature of the structure.Part of the heat conducted by thermal convection is dissipated to the environment or other parts of the body surface, and the other part enters the interior of the structure through thermal conduction. When the outer surface temperature is high, most of the heat will be dissipated into the atmosphere by thermal radiation, only a small part will enter the structure by thermal conduction. Therefore,qcond can be ignored, and formula (3) can be simplified to:
(4)T=qconvεσ+Tf44.Formula (4) is also called the radiation balance formula. When most of the heat is transferred by thermal radiation, the temperature of the outer surface of the structure can also be estimated by formula (4). During the whole time domain, the maximum temperature Tmax of the outer surface of the structure is obtained when qconv takes the maximum value. Before automatic material selection, the thermal properties of various optional thermal protection materials should be input first, including the material upper limit temperature Tlim. Then, the material whose upper limit temperature is closest to the maximum temperature of the outer surface of the structure will be taken, that is, the target material, as shown in
(5)minTmax−Tlim.In the automatic material selection module of the engineering software, the maximum temperature of the outer surface of the structure at each structural point is calculated according to formula (4), and then the material of each design point will be selected by formula (5). The extreme temperature of the material and other material properties are included in the input material library.


## 3.3. One-Dimensional Heat Transfer Algorithm of Multilayer Plate Based on Controlled Volume Method

Before the calculation and analysis of TPS temperature, a suitable and simplified heat transfer model should be established first to achieve the purpose of ensuring both accuracy and calculation efficiency. Generally, there are three types of thermal conductivity models: one-dimensional model, two-dimensional model, and three-dimensional model. The advantages and disadvantages of the three heat transfer models are shown in Table1.Table 1
Comparison of advantages and disadvantages of three heat transfer models.
One-dimensional heat transfer modelTwo-dimensional heat transfer modelThree-dimensional heat transfer modelAdvantages(1) The model is simple(2) The calculation is efficient and reliable(1) It can simulate temperature changes and heat transfer in both thickness and lateral directions simultaneously(2) The calculation accuracy of the two-dimensional model is higherIt can simulate the heat conduction process in three directions.Disadvantages(1) The lateral temperature change cannot be simulated(2) It affects the accuracy under some complex boundary conditions(1) The calculation speed is slow(2) It is not suitable for the rapid design of TPS in the conceptual design stage(1) The model establishment is complicated(2) The calculation amount is huge(3) The efficiency is extremely lowConsidering the advantages and disadvantages of the three heat transfer models, thermal insulation tiles are focused on the study of heat conduction and temperature distribution along the thickness direction. Additionally, the temperature difference between adjacent tiles is small, radiation effects from the sides of the tiles can be ignored, so a one-dimensional heat transfer model is sufficient.The thermal protection layer may contain multiple layers of different materials, and a multilayer one-dimensional heat transfer analysis calculation should be performed during the thickness optimization of the thermal protection layer. In the heat transfer analysis module, this software adopts the control volume method which is currently used in mainstream commercial software. The method focuses on the integral balance of the control volume, and it used nodes as the representative of the control volume [19]. Firstly, in the control volume method, the variation law of the variable which needs to be determined in the control volume, and the distribution function in the control volume are set. Then, the control quantity is integrated according to the characteristics of the distribution function, and the algebraic formula between the node variables can be obtained. Finally, the formula will be solved and the calculation results of each node will be got. Since there is no gap between the objects, energy must be conserved, so the calculation results obtained by the control body method must satisfy the principle of energy conservation. In other words, the control volume method can be solved for both uniform and nonuniform meshes. Therefore, the control volume size has a wide range and high flexibility without losing physical authenticity.When controlling the balance of the volume integral, the selected variable distribution usually has two forms: stepped distribution and piecewise linear distribution. As shown in Figure10, it is assumed that the x-axis is the one-dimensional transmission coordinate axis of temperature, and the T-axis is the temperature that changes with time τ, that is, the variable. Considering the gradient value of the calculated variable at the control volume interface, the stepped distribution is generally only used for the distribution of source terms, physical parameters, and variables in the time domain. Assuming three adjacent nodes W, P, and E, the control volume interface is represented by a dashed line, the location of which is represented by w and e, and is usually at the midpoint of the two nodes.Figure 10
Two distribution assumptions in the control volume method.For no internal heat source, the one-dimensional unsteady heat conduction control formula, in which the temperature is only conducted along the thickness direction of the multilayer plate, is shown in the following formula(6)ρc∂T∂τ=λ∂2T∂x2,where ∂T/∂τ,∂T/∂x are changing rate of temperature with time and the changing rate of temperature along the one-dimensional thickness direction, and λ and c are the thermal conductivity and specific heat capacity of the material, respectively.As shown in Figure11, for a uniform grid divided by a multilayer plate with no internal heat source and the temperature is only conducted along the thickness direction, the grid length is Δx, and three adjacent nodes n−1,n,n+1 are selected for the research. ΔX represents the control volume P, e, and w are the boundary points of the control volume P. For a uniform grid, e is the midpoint between nodes n and n+1, and w is the midpoint between nodesn−1 and n, so
(7)ΔX=Δx.Figure 11
Mesh nodes and control volumes.Integrating the control volumeP from time τ to τ+Δτ, we can get
(8)∫we∫ττ+Δτρc∂T∂τdτdx=∫ττ+Δτ∫we∂∂xλ∂T∂xdxdτ.When the temperature of the control volumeP is represented by the node temperature, the distribution of the temperature in the time and space coordinates is a step shape, and the integral on the left side of formula (8) can be expressed as
(9)∫we∫ττ+Δτρc∂T∂τdτdx=ρcPTnτ+Δτ−TnτΔx.When the distribution of temperatureT between nodes is piecewise linear, the integral on the right side of formula (8) can be expressed as
(10)∫ττ+Δτ∫we∂∂xλ∂T∂xdxdτ=∫ττ+ΔτλeTn+1−TnΔx−λwTn−Tn−1Δxdτ.Combining formula (9) and formula (10), and we can get the simplified governing formula
(11)ρcPTnτ+Δτ−TnτΔx=∫ττ+ΔτλeTn+1−TnΔx−λwTn−Tn−1Δxdτ.In formula (11), ρcP is the average density and specific heat capacity in the control volume P, and the specific heat capacity varies with the temperature, λe and λw are the equivalent thermal conductivity at points e and w, respectively, and the process of solution is as follows.In Figure12, e represents the upper boundary of the control volume P, Δxe+ is the distance from e to n+1, and Δxe− is the distance from e to n. The heat flow through the interface at point e can be expressed as
(12)q=Te−Tn+1Δxe+/λn+1=Tn−TeΔxe−/λn=Tn−Tn+1Δxe+/λn+1+Δxe−/λn.Figure 12
Control volume interface.In formula (12), q represents the heat flow through the interface at point e, Tn+1, Tn, and Te represent the temperatures at points n+1, n, and e, respectively, λn+1, λn, and λe represent the thermal conductivity at points n+1, n, and e. The heat flux density between nodes n+1 and n can also be expressed by the interfacial thermal conductivity
(13)q=Tn−Tn+1Δx/λe.According to the compatibility principle, the heat flow of formula (12) and formula (13) should be equal, so
(14)Δxλe=Δxe+λn+1+Δxe−λn.In a uniform grid, wheree is the midpoint of n and n+1, we can get
(15)λe=2λn+1λnλn+1+λn.The thermal conductivity of pointe can be obtained by formula (15), and the thermal conductivity of point w can be obtained in the same way. However, although the thermal conductivity of two points is obtained, the integral with respect to time on the right side of formula (11) cannot be obtained, and further approximation is required.Figure13 is a schematic diagram of the weighting factor, and f is the weighting factor. By introducing the weighting factor, formula (11) can be simplified as
(16)ρcPTnτ+Δτ−TnτΔx=fλeTn+1τ+Δτ−Tnτ+ΔτΔx−λwTnτ+Δτ−Tn−1τ+ΔτΔx+1−fλeTn+1τ−TnτΔx−λwTnτ−Tn−1τΔx.Figure 13
Weighting factor of time steps.As shown in Figure13, f is between 0 and 1. f=0 means that the temperature maintains the value of time τ in the interval of Δτ, and does not jump to Tnτ+Δτ until the time is τ+Δτ, when f=1, the temperature jumps to Tnτ+Δτ at the time of τ and remains unchanged throughout the time of Δτ, when f=0.5, the temperature T will be distributed piecewise linearly over time.Whenf=0, the time τ is expressed as time i, and the time τ+Δτ is expressed as time i+1, the formula (16) can be simplified as
(17)Tni+1=ΔτρcPΔx2λeTn+1i−Tni−λwTni−Tn−1i+Tni.In formula (17), the temperature of each node at the next moment i+1 is obtained from the temperature of the three adjacent nodes including the node at the moment i, and there is no need to solve the algebraic formulas simultaneously. Formula (17) is the discretization formula that controls the display format of the volume method.Formula (17) needs to satisfy the principle of positive coefficients, that is, the coefficients of Tni+1 and Tni are the same as positive. Otherwise, physically unreal solutions will appear, resulting in nonconvergence of computational results, so the coefficients Tni in formula (17) must be positive, that is
(18)1−λw+λeΔτρcΔx2>0.Formula (18) can be further simplified and we can get
(19)Δτ<ρcΔx2λw+λe.It can be seen from formula (19) that when Δx is reduced to 0.5 times of the original, the maximum value of Δτ will be reduced to 0.25 times of the original, that is, reduced by the quadratic, which will greatly reduce the calculation efficiency, so it is unrealistic to reduce the value of the grid length to improve the calculation accuracy.Whenf=1, the time τ is expressed as time i, and the time τ+Δτ is expressed as time i+1, formula (16) can be simplified as
(20)λwTn−1i+1−λw+λe+ρcPΔx2ΔτTni+1+λeTn+1i+1=−ρcPΔx2ΔτTni.The characteristic of formula (20) is that the temperature of each node at the next moment i+1 is obtained from the temperature of the moment i+1 including the three adjacent nodes of the node, so each algebraic equation has three unknowns, which cannot be solved directly. It is necessary to list the algebraic equations of all nodes, and solve the algebraic equation system simultaneously, which is often solved by tridiagonal matrix algorithm (TDMA). Formula (20) is called the discretization equation of the fully implicit format of the control volume method. Carefully observe the formula (20), and we can see that the coefficients of Tni+1 and Tni must be both negative, so Δx and Δτ have no constraint relationship, it must have a real solution in the physical sense, that is, the calculation results must converge.TDMA is used to calculate the case where the coefficient matrix is arranged in a tridiagonal arrangement, including two calculation procedures: the elimination of variables to find the coefficient and back-substitution to find the temperature of the node. Firstly, eliminate one variable, so that only two variables remain in the equation, then these two variables are obtained through boundary conditions, all two-variable results are calculated through iteration, and finally, the temperature of each node is calculated by back-substitution again.Suppose the discretization equation is(21)aiTi−1+biTi+ciTi+1=di,i=1,2,⋯,N,(22)b1c1a2b2c2⋯⋯⋯aN−1bN−1cN−1aNbNT1T2⋯TN−1TN=d1d2⋯dN−1dN.The matrix expression of formula (21) is shown in formula (22). The entire matrix only has coefficients around the diagonal and on both sides of the diagonal, forming a tridiagonal matrix.Wheni=1, we have
(23)T1=−c1b1T2+d1b1=P1T2+Q1.Wheni=2, we have
(24)T2=−c2a2P1+b2T3+d2−a2Q1a2P1+b2=P2T3+Q2.The discrete equation for anyi is
(25)Ti=−ciaiPi−1+biTi+1+di−aiQi−1aiPi−1+bi=PiTi+1+Qi.TheNth equation is
(26)aNTN−1+bNTN=dN.SubstituteTN−1=PN−1TN+QN into formula (26) and simplify, we get
(27)TN=dN−aNQN−1aNPN−1+bN=QN.Therefore, the general form of the two-variable equation is(28)Ti=PiTi+1+Qi,i=1,2,⋯,N.In formula (28):
(29)Pi=−ciaiPi−1+bi,Qi=di−aiQi−1aiPi−1+bi.The coefficients of the boundary nodes are as follows, becausea1=0, so
(30)P1=−c1b1,Q1=d1b1.BecausecN=0, so
(31)PN=0,QN=dN−aNQN−1aNPN−1+bN=TN.Therefore, in the positive process, only all the coefficientsPi and Qi need to be calculated, and finally the QN value is obtained. Since tN=QN, all node temperatures can be calculated by applying the general form of the two-variable equation Ti=PiTi+1+Qi in the inverse process back-substitution.It can be seen from the above operation process that in order to obtain allPi and Qi, the conditions of bi≠0 and aiPi−1+bi≠0 must be satisfied. The system of equations has a unique solution only when the conditions are met.Whenf=0.5, the time τ is expressed as time i, and the time τ+Δτ is expressed as time i+1, the formula (16) can be simplified as
(32)λwTn−1i+1−λw+λe+2ρcPΔx2ΔτTni+1+λeTn+1i+1=λwTni−Tn−1i+λeTni−Tn+1i−2ρcPΔx2ΔτTni.Formula (32) is very similar to the combination of formula (17) and formula (20). It needs to solve the algebraic equation system and also needs to satisfy the principle of positive coefficients. Formula (32) is called the discretization equation of the Crank-Nicolson format (C-N format) of the control volume method. The coefficients of Tni+1 and Tni have the same sign (both positive or negative), so
(33)Δτ<2ρcΔx2λw+λe.Comparing formula (19) and formula (33), when the value of Δx is determined, and the discrete control equation in C-N format is used, the value of Δτ is twice the value of the display format. Carefully explore the variation law of temperature change with time, we can see that at small time steps, the C-N format is closer to a linear variation, and the accuracy will be higher at this time. Undoubtedly, taking into account the computational accuracy and physical authenticity, it is ideal to use the exponential change law, but this will increase the difficulty of dealing with specific problems.Using above three discrete methods, we can obtain the governing formulas describing the heat flow inside the structure, but the boundary conditions are also required to provide initial values for the governing formulas, so as to solve the temperature values of all nodes of the structure at different times. Taking the multilayer plate boundary condition in Figure14 as an example, we can use the law of conservation of energy to discretize the boundary conditions.Figure 14
Boundary conditions on the surface of the structure.The boundary conditions on the surface of the structure is shown in Figure14. There are two heat transfer ways: heat convection and heat radiation. The boundary node N represents an element with a width of Δx/2 (the shaded part in Figure 14), and the temperature presents a stepped distribution with time. For the shaded part, the energy is conserved, so we can get
(34)λTN−1i−TNiΔx+hTf−TNi+εσTf4−TN4i=ρcΔx2TNi+1−TNiΔτ.In formula (34), h is the convection coefficient that changes with time, ε is the emissivity of the material on the boundary, which is a fixed value, and Tf is the ambient temperature.Simplifying formula (34) further, we can get
(35)TNi+1=−2εσΔτρcΔxTN4i+1−2λΔτρcΔx2−2hΔτρcΔxTNi+2λΔτρcΔx2TN−1i+2εσΔτρcΔxTf4+2hΔτρcΔxTf.It can be seen from formula (35) that the temperature of the boundary node at the next moment i+1 is obtained from the temperature of three adjacent nodes including the node at the current moment i, and the temperature value of the node at all moments can be calculated by iterative calculation.


## 4. Tile-Type Thickness Optimization Algorithm

In the tile-type thickness optimization module, the optimization parameter is the thickness of tile-type, and within the duration of the entire heat flow-time envelope curve, the maximum limiting temperature of the inner surface of the structure is set as the constraint condition. This is an optimization problem of a one-variable function, the one-dimensional search method [20] is a numerical method that uses an iterative method to find the minimum value of a one-variable function.
(36)minφα.Whenφα is differentiable, theoretically the optimal solution to this problem can be obtained by φ′α=0. However, this formula tends to be highly nonlinear and it is difficult to find an analytical solution. Moreover, in many practical problems, it is not differentiable, or its derivative expression cannot be written, so iterative numerical method is generally used to solve the above-mentioned minimum value problem, which is the so-called one-dimensional search. One-dimensional search method includes one-point method, two-point method, and 0.618 method. The two-point method requires the objective function value and its first derivative value of two points.Assuming that an interval [α1, α2] has been determined, and
(37)φ′α1≤0,φ′α2≥0.The core of the two-point method is to gradually reduce the interval [α1, α2] until the minimum value of the objective function is obtained with sufficient accuracy. The way to narrow the interval is to use formula (38) to obtain a new estimate of the minimum value at each iteration:
(38)α3=α2−ρα2−α1.In the formula (38), ρ∈0,1, the methods for determining ρ usually include the chord position method and the dichotomy method.The chord method usesφ′α1, φ′α2 and linearly interpolates between α1 and α2:
(39)φ′α=φ′α1+φ′α2+φ′α1α2−α1α−α1.The pointα of φ′α=0 in formula (39) is approximated as α3:
(40)ρ=φ′α2φ′α2−φ′α1The dichotomy rule can simply takeρ as 0.5, that is
(41)ρ=α1+α22The new interval should be determined based on the sign ofφ′α3. Since φ′α1<0 has been assumed, if φ′α3>0, the new interval should be taken as [α1, α3], that is, the optimal point should be between α1 and α3, if φ′α3<0, the new interval should be [α3, α2]. The process of narrowing the interval will be carried out until the accuracy requirements of the actual problem are met or the length of the interval is less than the specified error limit ε. The specific calculation process is shown in Figure 15.Figure 15
1D search iteration flow chart.Because the optimization of thickness is a monotonic problem, the temperature can be directly compared to determine the thickness range. The iterative method is similar to the one-dimensional search, which is called quasi-one-dimensional search. The specific operation method is as follows:(1)
It is needed to give the initial thickness firstly, and using the bisection method to find the second thickness to ensure that the optimal thickness value is between these two thicknesses(2)
To narrow the thickness interval, it needs to use chord method to find a new thickness through linear interpolation, and form a new thickness interval(3)
If the error is too large, it is necessary to repeat the previous step until the thickness interval reaches enough accuracy, or just to find the thickness value whose calculation result is equal to the optimal temperature


## 5. Verifying Example

In order to verify the correctness and engineering practicability of the tile-type TPS thickness automatic optimization software developed in present paper, we selected a half-wing model, and used UIIA and TTTPSTODS (developed in this paper), respectively, to optimize the structure. Finally, the results from two methods are compared and demonstrated.The half-wing model schematic is shown in Figure16. The upper and lower surfaces of the wing are subjected to time-varying heat flow loads. The model has a total of 92 structural nodes and 78 structural meshes. Material selection and thickness optimization are required for each structural point.Figure 16
The half-wing model schematic.For the optimal design model of the half-wing, the mathematical model of each design point can be expressed as(42)minTlimit−T,s.t.T−Tlimit≤0,xmin≤x≤xmax,where x is the variable, which represents the material thickness at each point, T represents the temperature of the inner surface of the skin at each node. Tlimit represents the preset limit temperature of the inner surface of the skin, xmin and xmax represent the minimum and maximum thickness of the insulation tile at each design point, respectively.


### 5.1. Analysis Process

When using UIIA method, we use the NLPQL [15] (Non-Linear Programming by Quadratic Lagrangian) algorithm, NLPQL is a sequential quadratic programming (SQP) method for solving problems with smooth continuously differentiable objective functions and constraints. It expands the objective function with a second-order Taylor series, linearizes the constraints, and obtains the next design point by solving the quadratic programming. Its general mathematical model can be expressed as
(43)minFxObjectivefunction,s.t.gjx=0j=1,2,⋯,mEqualityconstraints,gjx≥0j=1,2,⋯,mInequalityconstraints,xmin≤xi≤xmaxi=1,2,⋯,nBoundaryconstraints.In the TTTPSTODS software, according to the heat flow-time curve envelope method mentioned in Section3.1, we get the heat flow-time enveloping curve of the wing (Figure 17). And we can get the heat flow distribution on the upper and lower surfaces of the half-wing (Figure 18), then, we can set the initial temperature, ambient temperature, emissivity, and target optimization temperature of the wing structure on the boundary condition interface. And to obtain a load case, it is necessary to select and bind the heat flow data and boundary conditions.Figure 17
The heat flow-time enveloping curve.Figure 18
The distribution of heat flow on the wing surface.
(a)
Upper wing surface(b)
Lower wing surfaceThen, the software divides the outer surface of the wing into low temperature, medium temperature, and high temperature areas according to the heat flow distribution, as shown in Figure19. The user needs to input materials for automatic material selection in the material library, the material properties of “density,” “limit temperature,” “ specific heat capacity ,“ and so on. The material properties used in this example are listed in Table 2.Figure 19
The division of material areas.
(a)
Low temperature area(b)
Middle temperature area(c)
High temperature areaTable 2
Material properties used in this paper.
MaterialDensity (kg/m3)Upper limit temperature (K)Temperature (K)Specific heat capacity (J/kg/K)Thermal conductivity (W/m/K)RCG1820.51533.33116.66627.600.737172.22711.270.778255.55794.940.842394.43899.550.951533.301004.141.045672.211087.831.131811.121192.431.218950.011255.211.2971088.881317.951.3761200.011359.821.4341227.791380.711.4481338.881422.571.5071366.661443.471.5271422.211464.411.551450.021476.941.5641533.321506.231.6141644.431569.011.6891811.101631.751.8681922.211631.752.042LI-22003521755.56116.67292.880.055255.56627.60.074394.44878.640.093533.331054.370.113672.221150.60.132811.111150.60.1529501238.460.181088.891255.20.2061227.781267.750.2411366.671267.750.2771533.331267.750.3271644.441267.750.3671811.111267.750.4331922.221267.750.485HRSI144.181260116.662920.04172.22439.320.045255.55627.60.05394.43878.640.06672.211150.60.09811.121204.990.111088.881255.20.1631227.791267.750.19LRSI48.06773227.78751.20.025311.11950.10.034394.441052.30.042588.891129.680.05RCC1656.471811.11144.44502.082.30255.56711.283.76366.67878.644.99477.781004.165.99533.331087.846.30811.111297.047.291088.891422.567.651366.671548.087.721644.441673.607.651811.111715.447.65Innermost layer1656.47588.8988.8988.890.22200.00552.660.43311.11870.850.62422.221159.742.48Skin26105002739040.53739040.64739040.7The software then uses the radiation balance formula to automatically select materials based on the maximum heat flow value at each structural point. Finally, according to the calculation needs, the user can choose the display calculation format, the implicit calculation format and the C-N calculation format. The process is shown in Figure20.Figure 20
The calculation process of the half-wing model.


### 5.2. Results from Calculation

In this example, the constraint condition is set as the maximum temperature of the inner surface of the skin layer, the value of which does not exceed 449 K (can be set). Figures21–26 show the thickness distribution of the main protective layer, the maximum temperature distribution on the inner surface of the skin, and the maximum temperature distribution on the outer surface of the coating after optimization, by UIIA and TTTPSTODS developed in present paper, respectively.Figure 21
The thickness distribution of the main protective layer on the wing surface from UIIA.
(a)
Upper wing surface(b)
Lower wing surfaceFigure 22
The thickness distribution of the main protective layer on the wing surface from TTTPSTODS.
(a)
Upper wing surface(b)
Lower wing surfaceFigure 23
Temperature distribution on the inner surface of the wing skin from UIIA.
(a)
Upper wing surface(b)
Lower wing surfaceFigure 24
The temperature distribution of the inner surface of the wing skin from TTTPSTODS.
(a)
Upper wing surface(b)
Lower wing surfaceFigure 25
Temperature distribution on the outer surface of the wing surface coating from UIIA.
(a)
Upper wing surface(b)
Lower wing surfaceFigure 26
The temperature distribution of the outer surface of the wing surface coating from TTTPSTODS.
(a)
Upper wing surface(b)
Lower wing surface


### 5.3. Comparison and Analysis

From the above obtained results, we can see that in the high temperature region (leading edge) where the heat flow value is large, the thickness of the main protective layer is the largest, and the closer to the leading edge of the wing, the thickness of the main protective layer increases gradually, and the thickness of the main protective layer on the lower surface of the wing is much larger than that on the upper wing surface. Meanwhile, the temperature values of the inner and outer surfaces of the wing reach the maximum value at the leading edge, and the temperature of the lower wing surface is generally higher than that of the upper wing surface. And let us compare the optimization results obtained by UIIA and TTTPSTODS, we can find that the results obtained from the two methods are basically consistent with the numerical values, which indicated that the methods and algorithms used in each module of the TTTPSTODS are correct.Since the result graphs obtained by the two methods are very similar, it is difficult to find their differences directly from the graphs. In order to compare the data results obtained from the two methods more clearly and intuitively, the data of 7 structural grid points (the highest temperature area) on the leading edge of the wing were selected for comparison, as shown in Table3.Table 3
Thickness of tile-type TPS at each structural point of the leading edge of the wing.
Results from TTTPSTODSResults from UIIAGrid pointThickness/mmExternal surface temperature/KInner surface temperature/KGrid pointThickness/mmExternal surface temperature/KInner surface temperature/K1149.821111.124491151.311153.45430.622154.231192.374492155.121222.91445.633157.151231.154493158.331264.58448.604162.261258.544494163.541303.45435.245165.331293.814495166.871339.76446.176167.241332.864496169.201372.57447.787168.871348.064497171.251389.12445.92It can be seen from Table3 that the maximum relative error of the thickness optimization results obtained by the two methods is only about 1.39%, and the maximum relative error of the maximum temperature of the outer surface of the coating is only about 3.75%, which indicated that both methods have good convergence and accuracy.Under the condition that the accuracy requirements are met, the computational efficiency of the two methods is compared. When using UIIA to optimize the thickness of the main protective layer on the wing surface, the calculation time is about 10.5 h, while the calculation time using TTTPSTODS is about 0.5 h, if it is needed to use Abaqus to verify, it takes about 1.2 h, so the entire calculation process takes 1.7 h in total. It can be seen that the calculation speed of TTTPSTODS is much faster than that of Abaqus. For a half-wing model with only 92 nodes, it will take 10.5 h to optimize by using UIIA. If the whole aircraft is optimized for the tile-type TPS, the number of structural points will increase greatly. The time required for an optimization design using UIIA is immeasurable. Therefore, it is difficult to carry out engineering promotion for optimization using UIIA. However, if we only perform local verification, it only takes less than 10 h by using the one-dimensional heat transfer algorithm and thickness optimization algorithm in TTTPSTODS to calculate the model of thousands of structural points, and the calculation results are guaranteed. In the case of accuracy, the computational efficiency is greatly improved.In addition to the advantages of calculation efficiency, comparing with the optimization method of UIIA, TTTPSTODS also has functions such as interactive input and output, visualization, and the operation interface is clear, the operation method is flexible and convenient, and it is a very friendly interactive visual automatic engineering software.


## 5.1. Analysis Process

When using UIIA method, we use the NLPQL [15] (Non-Linear Programming by Quadratic Lagrangian) algorithm, NLPQL is a sequential quadratic programming (SQP) method for solving problems with smooth continuously differentiable objective functions and constraints. It expands the objective function with a second-order Taylor series, linearizes the constraints, and obtains the next design point by solving the quadratic programming. Its general mathematical model can be expressed as
(43)minFxObjectivefunction,s.t.gjx=0j=1,2,⋯,mEqualityconstraints,gjx≥0j=1,2,⋯,mInequalityconstraints,xmin≤xi≤xmaxi=1,2,⋯,nBoundaryconstraints.In the TTTPSTODS software, according to the heat flow-time curve envelope method mentioned in Section3.1, we get the heat flow-time enveloping curve of the wing (Figure 17). And we can get the heat flow distribution on the upper and lower surfaces of the half-wing (Figure 18), then, we can set the initial temperature, ambient temperature, emissivity, and target optimization temperature of the wing structure on the boundary condition interface. And to obtain a load case, it is necessary to select and bind the heat flow data and boundary conditions.Figure 17
The heat flow-time enveloping curve.Figure 18
The distribution of heat flow on the wing surface.
(a)
Upper wing surface(b)
Lower wing surfaceThen, the software divides the outer surface of the wing into low temperature, medium temperature, and high temperature areas according to the heat flow distribution, as shown in Figure19. The user needs to input materials for automatic material selection in the material library, the material properties of “density,” “limit temperature,” “ specific heat capacity ,“ and so on. The material properties used in this example are listed in Table 2.Figure 19
The division of material areas.
(a)
Low temperature area(b)
Middle temperature area(c)
High temperature areaTable 2
Material properties used in this paper.
MaterialDensity (kg/m3)Upper limit temperature (K)Temperature (K)Specific heat capacity (J/kg/K)Thermal conductivity (W/m/K)RCG1820.51533.33116.66627.600.737172.22711.270.778255.55794.940.842394.43899.550.951533.301004.141.045672.211087.831.131811.121192.431.218950.011255.211.2971088.881317.951.3761200.011359.821.4341227.791380.711.4481338.881422.571.5071366.661443.471.5271422.211464.411.551450.021476.941.5641533.321506.231.6141644.431569.011.6891811.101631.751.8681922.211631.752.042LI-22003521755.56116.67292.880.055255.56627.60.074394.44878.640.093533.331054.370.113672.221150.60.132811.111150.60.1529501238.460.181088.891255.20.2061227.781267.750.2411366.671267.750.2771533.331267.750.3271644.441267.750.3671811.111267.750.4331922.221267.750.485HRSI144.181260116.662920.04172.22439.320.045255.55627.60.05394.43878.640.06672.211150.60.09811.121204.990.111088.881255.20.1631227.791267.750.19LRSI48.06773227.78751.20.025311.11950.10.034394.441052.30.042588.891129.680.05RCC1656.471811.11144.44502.082.30255.56711.283.76366.67878.644.99477.781004.165.99533.331087.846.30811.111297.047.291088.891422.567.651366.671548.087.721644.441673.607.651811.111715.447.65Innermost layer1656.47588.8988.8988.890.22200.00552.660.43311.11870.850.62422.221159.742.48Skin26105002739040.53739040.64739040.7The software then uses the radiation balance formula to automatically select materials based on the maximum heat flow value at each structural point. Finally, according to the calculation needs, the user can choose the display calculation format, the implicit calculation format and the C-N calculation format. The process is shown in Figure20.Figure 20
The calculation process of the half-wing model.


## 5.2. Results from Calculation

In this example, the constraint condition is set as the maximum temperature of the inner surface of the skin layer, the value of which does not exceed 449 K (can be set). Figures21–26 show the thickness distribution of the main protective layer, the maximum temperature distribution on the inner surface of the skin, and the maximum temperature distribution on the outer surface of the coating after optimization, by UIIA and TTTPSTODS developed in present paper, respectively.Figure 21
The thickness distribution of the main protective layer on the wing surface from UIIA.
(a)
Upper wing surface(b)
Lower wing surfaceFigure 22
The thickness distribution of the main protective layer on the wing surface from TTTPSTODS.
(a)
Upper wing surface(b)
Lower wing surfaceFigure 23
Temperature distribution on the inner surface of the wing skin from UIIA.
(a)
Upper wing surface(b)
Lower wing surfaceFigure 24
The temperature distribution of the inner surface of the wing skin from TTTPSTODS.
(a)
Upper wing surface(b)
Lower wing surfaceFigure 25
Temperature distribution on the outer surface of the wing surface coating from UIIA.
(a)
Upper wing surface(b)
Lower wing surfaceFigure 26
The temperature distribution of the outer surface of the wing surface coating from TTTPSTODS.
(a)
Upper wing surface(b)
Lower wing surface


## 5.3. Comparison and Analysis

From the above obtained results, we can see that in the high temperature region (leading edge) where the heat flow value is large, the thickness of the main protective layer is the largest, and the closer to the leading edge of the wing, the thickness of the main protective layer increases gradually, and the thickness of the main protective layer on the lower surface of the wing is much larger than that on the upper wing surface. Meanwhile, the temperature values of the inner and outer surfaces of the wing reach the maximum value at the leading edge, and the temperature of the lower wing surface is generally higher than that of the upper wing surface. And let us compare the optimization results obtained by UIIA and TTTPSTODS, we can find that the results obtained from the two methods are basically consistent with the numerical values, which indicated that the methods and algorithms used in each module of the TTTPSTODS are correct.Since the result graphs obtained by the two methods are very similar, it is difficult to find their differences directly from the graphs. In order to compare the data results obtained from the two methods more clearly and intuitively, the data of 7 structural grid points (the highest temperature area) on the leading edge of the wing were selected for comparison, as shown in Table3.Table 3
Thickness of tile-type TPS at each structural point of the leading edge of the wing.
Results from TTTPSTODSResults from UIIAGrid pointThickness/mmExternal surface temperature/KInner surface temperature/KGrid pointThickness/mmExternal surface temperature/KInner surface temperature/K1149.821111.124491151.311153.45430.622154.231192.374492155.121222.91445.633157.151231.154493158.331264.58448.604162.261258.544494163.541303.45435.245165.331293.814495166.871339.76446.176167.241332.864496169.201372.57447.787168.871348.064497171.251389.12445.92It can be seen from Table3 that the maximum relative error of the thickness optimization results obtained by the two methods is only about 1.39%, and the maximum relative error of the maximum temperature of the outer surface of the coating is only about 3.75%, which indicated that both methods have good convergence and accuracy.Under the condition that the accuracy requirements are met, the computational efficiency of the two methods is compared. When using UIIA to optimize the thickness of the main protective layer on the wing surface, the calculation time is about 10.5 h, while the calculation time using TTTPSTODS is about 0.5 h, if it is needed to use Abaqus to verify, it takes about 1.2 h, so the entire calculation process takes 1.7 h in total. It can be seen that the calculation speed of TTTPSTODS is much faster than that of Abaqus. For a half-wing model with only 92 nodes, it will take 10.5 h to optimize by using UIIA. If the whole aircraft is optimized for the tile-type TPS, the number of structural points will increase greatly. The time required for an optimization design using UIIA is immeasurable. Therefore, it is difficult to carry out engineering promotion for optimization using UIIA. However, if we only perform local verification, it only takes less than 10 h by using the one-dimensional heat transfer algorithm and thickness optimization algorithm in TTTPSTODS to calculate the model of thousands of structural points, and the calculation results are guaranteed. In the case of accuracy, the computational efficiency is greatly improved.In addition to the advantages of calculation efficiency, comparing with the optimization method of UIIA, TTTPSTODS also has functions such as interactive input and output, visualization, and the operation interface is clear, the operation method is flexible and convenient, and it is a very friendly interactive visual automatic engineering software.


## 6. Conclusion

In the present paper, an efficient automatic engineering software is developed for the thickness optimization of tile-type TPS of hypersonic aircraft, which adopts a reliable heat flow-time envelope method and develops numerical calculation method for thermal protection tile thickness optimization. The software has an interactive visual interface, which is flexible and convenient to use. The calculation accuracy can be set according to specific requirements, the calculation time can be controlled, and the design efficiency can be improved.Finally, a half-wing model is designed, using the Isight integrating Abaqus optimization method and the tile-type TPS thickness optimization software developed in this paper, respectively. The design results are compared and analyzed, which verifies the calculation correctness and calculation efficiency of the engineering software developed in this paper. It is proved that the tile-type TPS thickness automatic optimization design software developed in the present paper is a reliable and practical engineering design software.



---

*Source: 3663567-2022-10-14.xml*