Thermal Protection System Design of a Reusable Launch Vehicle Using Integral Soft Objects

In the present paper, a modelling procedure of the thermal protection system designed for a conceptual Reusable Launch Vehicle is presented. A special parametric model, featuring a scalar field irradiated by a set of bidimensional soft objects, is developed and used to assign an almost arbitrary distribution of insulating materials over the vehicle surface. The model fully exploits the autoblending capability of soft objects and allows a rational distribution of thermal coating materials using a limited number of parameters. Applications to different conceptual vehicle configurations of an assigned thickness map, and material layout show the flexibility of the model. The model is finally integrated in the framework of a multidisciplinary analysis to perform a trajectory-based TPS sizing, subjected to fixed thermal constraints.


Introduction
In the last decade, a growing number of projects have been focused on the development of fully Reusable Launch Vehicles (RLV), designed for a crew-rescue mission.Successful demonstrative flights of private companies, like SpaceX, Virgin Galactic, and Sierra Nevada Corporation, and the activity promoted by the European Space Agency are aimed at improving operability of RLV [1,2].Consequently, a great deal of research effort has been put to design RLV as blended wing-bodies, because of the promising trade-off between aerodynamic efficiency, cross-range, and aeroheating performances during the reentry [3].The EXPERT (European eXPErimental Reentry Testbed) program and the IXV (Intermediate eXperimental Vehicle), which performed an atmospheric lifting reentry from orbital speed, are just examples of such demonstrators developed to predict performances of a full-scale vehicle.Besides, the X-37B, an unmanned lifting body developed by Boeing, has been put into orbit by an Atlas V rocket, performing a successful lifting-guided reentry.Furthermore, the foreseeable opportunity for space tourism represented by experimental flights of Virgin Galactic's SpaceShipTwo and SpaceX has also emerged [4].The requirements currently considered for RLV design are (i) to allow a very low-g (nearly 1.5 g) reentry, with a landing on a conventional runway; (ii) to adopt a light-weight (passive), fully reusable thermal protection system to withstand several flights without any replacement; and (iii) to provide vehicle autonomy to land at predefined locations for rescue issues [2,3].In order to fulfill all those requirements, the duration of reentry flight increases and consequently the integrated heat load absorbed by the structure [2].
The above circumstances may conflict with the adoption of a fully reusable TPS, eventually restricting the choice of insulating materials and penalizing the total mass [5].However, antithetical requirements between room for the payload, weight, and vehicle operability demand a trade-off between vehicle shape and TPS sizing.In preliminary design practice, thousands of design configurations are typically evaluated by an optimization algorithm to find the best fit [5][6][7][8][9][10].As a consequence, a preliminary appraisal of vehicle performances is commonly performed using high-efficiency, low-order fidelity methods that give support to a multidisci-plinary analysis performed with a computational effort which fit the typical timeline of the conceptual design phase [11].
The aerothermal environment is a basic design criterion for either TPS sizing or choice of materials [12,13].TPS sizing is generally performed once the reentry trajectory is assigned, having computed the peak heating flux and the time integrated heat load [14].External rocket propulsion systems allow an RLV a less severe heating due to their different ascent trajectories.
Therefore, an RLV operates an unpowered crew rescue with the TPS sized mainly by its aerothermal reentry corridor.On the other hand, an RLV that performs an ascent phase with air-breathing propulsion systems is sized considering more severe heating conditions [15].Several works dealing with TPS sizing have been published in the literature.Lobbia [8] determined the sizing of a TPS in the framework of a multidisciplinary optimization.Material densities and maximum reuse temperature were computed.TPS mass was estimated assuming the category of materials used for the space shuttle and thickness distribution assigned on a review of HL-20 materials for each component.Trajectorybased TPS sizing has been proposed by Olynick [13] for a winged vehicle concept.The peak heating temperature was determined considering an X-33 trajectory, discretized in a number of fixed waypoints.The resulting aerothermal database was used as an input for a one-dimensional conduction analysis, and several one-dimensional stackups of different materials representative of TPS were consequently dimensioned.Bradford [15] developed an engineering software tool for aeroheating analysis and TPS sizing.The tool is applicable in the conceptual design phase, for reusable nonablative thermal protection systems.The thermal model was based on a one-dimensional analysis, and TPS was modeled considering a stackup of ten different material layers.Mazzaracchio [14] proposed a method to identify the optimal combination between the ablative and the reusable part of TPS.This was dictated by a trade-off between the reentry trajectory and the TPS sizing.In Wurster and Stone [16], a selection of design trajectories, and wind tunnel data were integrated with CFD simulations and engineering prediction to perform TPS sizing for an HL-20 reentry vehicle.In the present work, following an approach proposed in Ref. [17], a soft objectderived representation for TPS thickness and material assignation is introduced.According to the legacy formulation of this technique, originally developed in computer graphics (C.G.) for the rendering of complex organic shapes [18], three-dimensional object surfaces are (implicitly) obtained by defining a set of source points (or even more complex varieties) irradiating a potential field that is subsequently rendered as a chosen isosurface.A major feature offered by this methodology is that the topology of the regions represented by these objects can be arbitrary and easily controlled just by altering the mutual spatial position of the primitives.These very desirable properties have made implicit modelling techniques extremely popular in the field of photorealistic rendering of organic shapes whose topology changes dynamically, and many advances and developments have been proposed over time to improve their capabilities and/or reduce some unwanted side-effects [19,20].In the present paper, following a quite different paradigm developed in [21], the full potential field irradiated by a set of bidimensional soft objects is congruently mapped on a discretized paneled RLV shape.The methodology is able to create arbitrary TPS distributions seamlessly increasing the thickness where critical heat loads are experienced and dropping out elsewhere.A similar, slightly modified procedure can also be applied to create an arbitrary binary map representing membership functions for TPS materials.This binary map can be operated independently on thickness distribution (or locally synchronized with different thickness maps).The present formulation is formalized in the framework of a parametric model which exploits simple variations of parameters to perform the soft object mapping over a discretized surface.Applications of the developed procedure to different vehicle configurations show the flexibility of the method.In addition, the procedure can be easily embedded in a multidisciplinary computation, to perform a trajectory-based TPS sizing on a RLV subjected to fixed thermal constraint.

Methodology
2.1.Soft Object Definition.Soft objects constitute a modelling technique introduced in computer graphics to represent three-dimensional objects having complex and organic shapes.Blinn first developed a soft object model to display the appearance of electron density clouds in a covalent bond of a molecule [18].According to the model formulation, curved (closed) surfaces can be modeled defining n ≥ 1 potential fields f i , namely, blobs.Several blobs f i can be connected smoothly by the self-blending property, by performing an algebraic summation of their potential fields [22]: The commonly adopted notation from the field function f i .The strength s i = f i d i accounts for the value of f i in a generic point x, at distance d i from a key point x i .
The potential field F generated by n blobs can be generally taken into account computing a specific isosurface S (subsequently processed) by a raster conversion algorithm: where the threshold T in equation ( 4) selects an isosurface of F. Blinn originally proposed the "blobby molecule," an isotropically decaying Gaussian function modulated by strength, and radius [18].The blobby molecule is a potential 2 International Journal of Aerospace Engineering function with infinite support.This aspect affects the computational effort in a practical implementation, because it has to be evaluated in all points of the space.However, in the literature, several finite support potential functions have been proposed for different modelling purposes [22].The field function f i used in the present work has a finite support and assumes normalized values in the range between 0 and 1 [22]: The parameter p defines the hardness of blob and controls the level of blending between two soft objects.
2.2.Two-Dimensional Integral Soft Object Field.Two-dimensional soft objects preserve the self-blending property.However, it is not always easy to create a rational distribution of a set of independent point source objects for application purposes.Therefore, blobs are conveniently and easily arranged in macroaggregates.Figure 1 shows that a potential field is created superposing n = 6 discrete blobs having radius r and centers mutually placed on a line segment of length l over an equally spaced two-dimensional grid of step δ e = 1/ n blob .If δ e < 2r two or more blobs superpose, the strength of the potential field is obtained summing up the strengths of each blob (yellow colored region).This procedure relies on a similar idea to the one developed in [21] to generate self-stiffened structural panels.Specifically, the full (integral) potential field irradiated by a set of discrete two-dimensional point source blobs generates a seamless potential field.This approach is quite different from isocontour tracking commonly adopted to represent soft objects.

Modelling of Two-Dimensional Stick Primitives.
A superposition of n point source blobs with key points placed on a geometric segment (straight or curved) is denoted from now on as a "stick."However, as shown in Figure 1(b), equation (1) creates a stick with a support having "bulges."Increasing the number of sticks, the shape of the support becomes more regular.However, the strength is not bounded (see Figure 1(b)).The above drawback is overcome by modifying the definition of the potential field given by equation (1).A bounded potential field, regardless of the number of the blobs used on a stick, is obtained with the relation Equation (6), where F 0 P = 0, expresses the global potential field F j P irradiated by a set of j blobs at a generic point P of space placed at a distance d from the key points, as the max between the previous j − 1 potentials accounted for by the assembly layer F j − 1 P and the current potential G j over the plane disk of radius r: where n blob is the total number of blobs present on the B-grid.Figures 2(a) and 2(b) show the strength field of a twodimensional stick primitive obtained using n blob = 6 and 20, respectively, computed with equation (6).By increasing the number of blobs on a stick, the strength of F is still bounded to a maximum unit value.Figures 2(c) and 2(d) show the same behavior for a tapered primitive having a linear variation of the blob radius along the axis of sticks.

RLV Shape Modelling and Thermal Protection System Sizing Criteria
In the present work, we assume that the generic shape of an RLV is represented by a grid formed by a quadrangular and/or by a degenerated triangular panel grid.Grid points are obtained using a proprietary procedure that authors fully detailed in [23,24].Without going into details of the shape model, we remark that the mesh arrangement over the RLV surface is obtained with no NURBS support surface: a three-dimensional parametric wireframe is created using cubic rational B-splines and used to reconstruct the computational surface grid.The control parameter allows a wide range of shape variations to handle different design objectives (thermal or dynamical) for a reentry mission.Grid topology is equivalent to a spherical surface with no singularities (open poles) and allows a mapping of the points in UV coordinates over an equivalent cylindrical surface.The above considerations ensure a topologically invariant shape.In previous papers proposed by the authors [23,24], a multidisciplinary shape optimization for an RLV comprising a trajectorybased TPS sizing procedure was developed.The TPS was modeled using two insulating materials placed at different locations along the vehicle surface.A different mapping (thickness distribution and longitudinal location) of the two materials with different operational temperature was adopted.The sizing of insulating materials required the computation of aerothermal loads across the reentry trajectory from a LEO orbit up to M ∞ = 2, which is considered the limit below which thermal heating can be neglected.TPS thickness was parametrically sized according to thermal requirements assumed in the optimization; a simple (but very rough), bilinear distribution of the TPS thickness along the longitudinal axis of the vehicle and a linear distribution across each cross section were adopted.The maximum allowable temperature values (depending on the adopted material) for the interior and exterior surfaces of TPS outlined the thermal constraint to be fulfilled by the sizing procedure.International Journal of Aerospace Engineering adjusted according to the normalized dimensions relative to this map.The topological map is emulated introducing a two-dimensional grid (from now on denoted as B-grid) having the same topology tree as the vehicle open grid (number of points, panels, and connectivity) except for the unit size.A geometric mapping between the B-grid and the vehicle grid is established, and elements of the Bgrid are biunivocally mapped onto corresponding elements of the vehicle surface (see Figure 3).Several stick primitives are emulated on the B-grid placing a number of n equally spaced isotropic blobs, with radius r and length l, in normalized units.Stick emulation is performed by overlapping n blobs using the special formulation reported in [21] that ensures a convergent envelope of the finite support and a limited value of the blob strength.An exemplificative spatial distribution of sticks on the B-grid is shown in Figure 3. Position and orientation of each stick is determined by assigning coordinates of centers C i and precession angles θ i , respectively, with respect to a Cartesian frame of reference Oxz oriented as in Figure 3. Therefore, a generic distribution of sticks created on the vehicle grid is equally mapped on the vehicle surface whatever the morphological map considered.In the present case, gray colored regions (1) denote points of the B-grid mapped on the windward side of the RLV shape (see Figure 3), while white regions (2) relate to leeward regions of the vehicle.Regions of the vehicle surface mainly subjected to heating peaks during the reentry maneuver are the (i) nose, (ii) leading edge, and (iii) tail.The global potential field generated by the sticks onto the B-grid is adjusted in a suitable dimensional scale and subsequently mapped on the mesh panels of the vehicle surface grid to obtain an easy and powerful control of the thickness distribution.The proposed methodology is able to create virtually arbitrary TPS distributions and can be easily tuned up to locally increase the thickness where critical heat loads are expected and dropping out elsewhere.A similar, slightly modified procedure is also applied to create an arbitrary binary map distribution of different TPS materials that may be operated independently of the thickness distribution.

Parametric Model of a Thermal Protection System
5.1.Thickness Modelling.As a demonstrative example, a parametric representation of a thermal protection system is obtained using a limited set of stick primitives (n stick = 5), oriented as shown in Figure 4. Skin sticks characterized by a large radius and limited strength are spread over the skin surface in the longitudinal direction in order to provide a thickness graded baseline.A constant minimum thickness is superposed in all remaining points of the B-grid, ensuring a nonzero value in any point of the grid.Furthermore, additional parametric sticks, specifically positioned and oriented to affect thickness in critical regions, such as the nose, leading edge, and trailing edge, complete the support for the TPS and create a rational distribution of insulating material suitable with a reentry mission.The parametric position of sticks and the axis of orientation are defined by assigning centroid coordinates x c and z c and angle θ th , measured with respect to the system of reference reported in Figure 4.  x c, q=1,2,3,4,5 = 0 0, 0 0, 0 0, 1 0, 1 0 , Skin (q = 1, 2) and nose (q = 3) sticks have a tapered support obtained by imposing a linear variation of the point source blob radius.Conversely, a constant radius is adopted for the leading edge (q = 4) and trailing edge (q = 5) sticks.

Material Modelling.
A similar but completely independent stick-based parameterization has been also defined to model a dynamic distribution map of different insulating materials, denoted here generically as material "0" and material "1."We assume that material "1" outperforms material "0."Therefore, material "1" is a natural candidate to insulate the nose, leading edge, and trailing edge.TPS materials are assigned according to a discrete distribution.Differently than sticks used for thickness distribution, this additional set of

6
International Journal of Aerospace Engineering primitives returns just binary values used to define specific materials.In this case, the field function m th (see relations in ( 9)) assumes a constant value equal to one inside the finite support of a stick and zero elsewhere.Centroid coordinates (mx c,q and mz c,q ) and length l q of these additional sticks are given by parametric relations in (9) with normalized parameters reported in Table 1: mx c, q=1,2,3,4,5 = 0 0, 0 0, 0 0, 1 0, 1 0 , mz c, q=1,⋯,5 = d q min + mt q ⋅ d q max − d q min , ml q=1,⋯,5 = mlt q ⋅ d q max , mth q=1,⋯,5 = 1 9

Additional Considerations about Integral Soft Objects for TPS Modelling
In order to better clarify the rationale underlying the proposed methodology, some additional considerations are provided next.As a general premise, it should be emphasized that implicit modelling techniques canonically used in C.G. differ in many respects from those reformulated and used in the present context.About this, a brief description of the approaches followed in C.G. is preliminarily given, to highlight both common points and main differences with respect to the illustrated methods.In C.G., implicit modelling is commonly used for the rendering of complex organic shapes.In these methods, some objects (usually referred to as blobs or primitives) of appropriate dimensionality (2D, 3D), typically represented by their own morphological skeletons, are conceived as "emitters" of suitable finite support field functions, expressed as distance laws in an appropriate norm (usually Euclidean).These primitives are allowed to mutually interact with each other by simply overlapping their finite supports, cumulating in that way the field intensities where superposition takes place.According to an implicit representation, a specific instance (isosurface) of the global field associated to an assigned isovalue can be finally visualized.Depending on the chosen rendering method (e.g., ray-tracing), the isosurface can preserve its implicit formulation or be evaluated through suitable progressive sampling algorithms (i.e., octree) to be translated into discrete polygonal elements, typically triangular meshes.In general, the mathematical structure of the isosurfaces will be characterized by smooth curvatures, highlighting an intrinsic capability to generate automatic fillets in those spatial regions where primitives overlap.Fillets can be controlled locally by introducing a hardness parameter in the field functions, which makes the primitives "harder" or "softer" in blending.As previously stated, in the present work, a methodology just roughly inspired by the aforementioned techniques has been introduced to arbitrarily control thickness and material assignment using a limited number of control parameters.The proposed approach uses the same finite support field functions employed in implicit modelling; therefore, the desirable capability to generate arbitrary topologies by simply controlling the mutual spatial position of primitives is , where the discrete resolution of the tracked isosurface is often adaptive and conditioned by the expected rendering quality and the amount of local curvatures, these primitives are now just projected onto a two-dimensional flat grid with a fixed resolution that acts as an extended finite support; the global field resulting from primitive interaction is then used integrally, not just represented with single contours.Therefore, no implicit formulation is used at all.These circumstances explain the definition of "integral soft objects" (opposed to "implicit modelling") used in the present context to describe the "agents" that operate to assemble the maps of thickness and material.In this new framework, the legacy control parameters for field functions adopted in C.G. are still maintained, although with a different semantic connotation.Specifically, the strength parameter has been reinterpreted either as the maximum thickness value of TPS locally transferred onto the topological map or the maximum integer value of the set of pointers; instead, the hardness parameter becomes the thickness gradient with respect to the distance in the finite support.A major aspect that affects the implementation of the method is given by the low spatial resolution of the topological map where TPS thickness and materials are transferred.Due to the specific methodologies (panel methods described next) used in the resolution of the thermal and fluid dynamic fields on the body at different speed regimes, this grid is extremely coarse if compared to the typical resolutions used in C.G. (only a few hundred quadrangular elements are used).In addition, the resulting thickness and material distributions are sampled only at the centroids of the panels and assumed constant for each panel.These intrinsic limitations have determined the above described implementation choices; those rationales are briefly given here: (i) a stick primitive has been preferred and systematically adopted in this context, because effective TPS thickness distributions likely happen along the UV parametric directions of the topological map.Effectiveness and flexibility of these distributions have been increased by also allowing tapered and steerable sticks.(ii) A surrogate representation of a skeleton-based stick is simply obtained by distributing a finite number of one-point skeleton primitives (circular blobs), mixing each other with the MAX blending operator outlined by equation (6).As mentioned earlier, the rationale of this choice is that the field function of the assembled stick primitive still remains limited, and at the same time, only a small number of circular blobs (slightly more than a dozen) are required to obtain a satisfactory envelope of the finite support without bulges because of the low intrinsic resolution of the grid combined with the actual ranges assigned to the finite supports of primitives.Moreover, this simplified approach allows in perspective the definition of some other primitives based on different skeletons (homeomorphic with the line segment, for example, based on splines) that might potentially ensure more sensitivity in the TPS definition problem with no major modification required.(iii) Although it is virtually possible to use any suitable blending function to merge the fields of different preassembled primitives (for example, the sum of the local magnitudes), the MAX blending function is also still preferred for this purpose, because the intrinsic capability offered by this operator to generate "ziggurat-style" step fields when the involved primitives exhibit different strength values has been considered desirable for an effective TPS thickness distribution.

An Example of TPS Modelling Capabilities
The previously introduced modelling procedure has been applied on a conceptual RLV shape created with the model described in Section 4 and detailed in [23,24].The applicability of the procedure is shown for the arbitrarily chosen distribution of stick primitives that creates a morphologically adaptive TPS on two RLV shapes with different dimensions: RLV-1 with length l tot = 9 8 m, wingspan w s = 5 6 m, and cabin height h = 1 6 m and RLV-2 with length l tot = 15 m, wingspan w s = 9 2 m, and cabin height h = 2 m.
The parameters characterizing the distribution of thickness and the materials are reported in Table 1.Figures 5(a) and 5(b) show the application of TPS modelling over the first configuration (RLV-1), on leeward (a) and windward (b) surfaces.Different colors denote different values of thickness and are represented in a dimensional scale.It can be observed that the thickness map can be easily tuned up for best covering of regions where maximum heat loads occurs (i.e., the nose and leading edge).Figure 5 shows the capability to create arbitrary seamless thickness distributions up to the value of the baseline thickness which has been arbitrarily set equal to th min = 0 05 m (denoted in blue color).This corresponds to a region of the leeward surface not covered by the skin stick.Figures 5(c) and 5(d) show the map of two different insulating materials created with equation (9).Red colors indicate material 1, which is placed on regions of the vehicle subjected to higher heat loads.Comparisons between Figures 6(a) and 6(b) and Figures 6(c) and 6(d) also exhibit the capability of the model to handle independently both the thickness and material distribution.Finally, Figures 6(a) and 6(b) and Figures 6(c) and 6(d) show the same blob distribution adopted either for thickness or for material modelling applied on a different RLV configuration (RLV-2).The procedure creates, as expected, the same TPS distribution both for thickness or materials on two different shapes and is completely independent by their morphology.

Application of the TPS Sizing Procedure for a Conceptual RLV Configuration
The TPS modelling procedure developed in the previous sections has been implemented in the ANSYS ® Parametric Design Language [25], to perform a trajectory-based sizing of a TPS for an RLV designed for a LEO reentry mission.A multidisciplinary analysis comprising aerodynamics, heating analysis, trajectory estimation, and mass estimation is implemented to determine the aerothermal loads on the vehicle.The entire flowchart of the procedure has been discussed and detailed in [23,24].For the sake of brevity, here, it will be addressed with reference to the specific application, specifying the assumptions adopted.TPS is designed to withstand heating for an unpowered reentry maneuver performed from   9 International Journal of Aerospace Engineering methodology for a conceptual design configuration, loworder fidelity methods [3,5] are used for each subdiscipline to reduce the computational time [11].
8.1.Aerodynamics.Aerodynamic coefficients for the hypersonic phase of the reentry are computed with a public domain panel flow solver using Newtonian flow theory available in Ref. [26], in the waypoints reported in Table 2.According to Newtonian theory, hypersonic flows are modeled as an ensemble of particles impacting the surface of a body approximated to a flat plate at incidence.Therefore, considering a discrete panelization of the body surface, the pressure force δF i acting on the i-th panel having an area δ A i and directed along the outward unit normal vector to the panel ni is then computed as δF i = −p i δA i ni .The nondimensional force coefficients in a body frame of reference are given by computed as the summation of pressure forces over the total number of mesh panels.Newtonian approximation is acceptable in the range of altitudes and Mach numbers up to M ∞ = 2, as viscous-force contribution in the axial direction on hypersonic bodies decreases with an increasing angle of attack [3,8,27].Validation of the present computations on aerodynamic coefficients has been previously performed in [24], together verifying the mesh independence of aerodynamic coefficients.Different choices of impact methods are adopted in the computation and are indicated in Table 3.
The rationale behind the partitioning of the nose and wingbody region is related to the different flow incidence occurring on the wing-body panels requiring the use of appropriate compression methods adopted along the windward region (see Table 3).A similarity between the flow on the wing-body panels and the flow downstream of an attached shock can be assumed, and the tangent cone method can be used [8].Aerodynamic coefficients of the vehicle at an incompressible Mach number are computed by adopting the integral formulation of the potential generated at a point P by a distribution of singularities (sources and doublets) with strength σ and μ, respectively, being S B and S W being the body and the wake surface, respectively.The freely distributed panel code available in Ref. [28] is adopted for the current computation of this aerodynamic regime.The solver requires a discretization of the geometry by quadrilateral panels with sources and doublets of constant strength.An approximation of equation ( 11) on the computational mesh is The unknown values of doublet strength μ i are obtained by solving the linear system of equations obtained by discretizing equation (11).To make the solution uniquely defined, tangency condition is applied, i.e., V ⋅ n = 0, and the source strength σ i is set equal to σ i = n i ⋅ V ∞ on each control node.Furthermore, a wake surface, still based on quadrilateral panels, has been added to apply the Kutta condition on each panel at the trailing edge.The intensity of the wake panels is set equal to the difference between the doublet strength of the upper and lower doublet panels at the trailing edge.Having computed the doublet strengths, the determination of velocity, pressure, and therefore aerodynamic coefficients is performed.8.2.Mass Estimation.The mass of the thermal protection system m tps has been estimated modelling the TPS with two of the materials adopted for the space shuttle thermal insulation system whose thermal properties are reported in [29,30]: the Reinforced Carbon-Carbon (RCC) composite and High-Temperature Reusable Surface Insulation (HRSI) tiles, made of coated LI-900 silica ceramics.RCC material can be selected for the thermal insulation of the nose, leading edge, or trailing edge, where the peak wall temperature is expected.A mass decoupling model, based on the set of relations reported in [31], is adopted in the current procedure to compute the mass of the vehicle: being m vehicle = m dry + m ecd + m av + m ecl + m par , 14 excluding the contribution of the thermal protection system, where m dry , m ecd , m av , m ecl , and m par are reported in Table 4.
The mass of the thermal protection system m tps , with ρ i , S i , and th i being the density of the insulating material, area of the discrete i-th panel of TPS, and thickness of the i-th TPS panel, respectively, computed with the procedure developed in Section 5.
8.3.Flight Dynamics.The calculation of the trajectory is performed in a nonrotating, inertial, Earth-Centered, Earth-Fixed (ECEF) frame of reference, according to the general system of equation of planetary flight reported in Ref. [32].The vehicle is assumed as a mass point, describing a nonplanar reentry trajectory with a constant bank angle μ a = 45 °, assigned to ensure a cross-range performance.A negative value of initial reentry flight-path angle γ = −1 4 °is given to avoid the skip phenomenon [32].The angle of attack changes International Journal of Aerospace Engineering following an implicitly defined modulation law, reported in Table 2.The initial value of longitude θ, latitude φ, and flight azimuth χ are all set to zero.The vehicle dynamics is described by a four-degree-of-freedom point mass model.Trajectory equations have been integrated by using an implicit Newton-Raphson method to reduce the computational time of the overall procedure, performing a time step sensitivity analysis to ensure convergence of the numerical integration.
8.4.Heating Analysis.The thermal state of the TPS is efficiently determined computing its exterior (w 2 ) and interior wall temperatures (w 1 ) which determines the choice of the material.The thermal analysis is only performed on the windward side as this is the maximum heated region of the TPS; the leeward side of the vehicle is supposed to be the fixed wall temperature region.Heat radiated from the external TPS wall w 2 provides the cooling of the vehicle.Catalytic recombination, low density effects, and thermal radiation from nonconvex surfaces are neglected.The above simplifying assumptions are considered reasonable for a conceptual design phase and determine a conservative overestimation of the TPS mass.A one-dimensional model of the TPS thickness [6] is used to determine TPS interior wall temperatures T w1 exploiting the exterior TPS surface T w2 temperature.
Radiative equilibrium is assumed on wall w 2 [5]: where ε is the surface emissivity and σ is the Stefan-Boltzmann constant.The convective heating at the wall w 2 at the stagnation point q w2 stag is approximated using the following correlation: where C pw2 is the specific heat at the wall, ρ the free-stream density, and V the vehicle velocity.The windward convective heating q w2 win on the vehicle is evaluated using a correlation for spheres, cylinders, and flat plates [33]: where the constant C is specifically characterized for a laminar or turbulent boundary layer.TPS wall temperatures are determined using kinematic trajectory data, in the range 2 ≤ M ∞ ≤ 23 where peak values were expected to occur.Equation ( 18) is solved using a Newton-Raphson method to determine the wall temperature T w2i on each panel of the TPS thickness discretized as shown Figure 7.The interior TPS wall temperature T w1i is computed at each node of the discrete i-th panel of the TPS integrating in time a onedimensional unsteady heat-diffusion model [12,23].Both the initial and boundary conditions assigned as T y w2,i , 0 = T y w1,i , 0 = 285 K, T y w2,i , t = T w,i , ∂T ∂y y w1i = 0, 19 provide a well-posed heat-diffusion problem, numerically solved with a finite difference method.The thermal state of TPS is globally defined assuming, for each computational element i, the following scalar relations: where it has been supposed that T max,w1 = 430 K is the value of the maximum allowable TPS temperature which adheres to structural elements of the vehicle and T max,w2 is the reuse temperature limit depending on the material [30].
T max,w2 RCC = 1920 K, T max,w2 Li-900 = 1644 K 22 Position, axis of orientation, and strength of stick primitives are manually and iteratively assigned, rationally choosing stick parameters, to cover the most heated regions of the TPS surface.A minimum value of the thickness set equal to th min = 0 05 m allows the convergence of the heat conduction problem (see equation ( 19)).A maximum value of thickness is chosen equal to th max = 0 29 m and takes into account the conservative assumptions made in Section 8.5; therefore, it is expected to be oversized but in agreement with the approximations used in the conceptual design phase [5].It is

Conclusions
In the present paper, a special modelling procedure of the thermal protection system designed for a conceptual Reusable Launch Vehicle has been developed.A set of macroaggregates of point source blobs organized in envelopes of finite supports and with a bounded strength has been successfully created on the topological map associated with the computational grid.Applications of the modelling procedure to different design configurations highlighted the sensitivity and powerful control to radically change the TPS using a limited number of parameters.This feature has been illustrated executing a multidisciplinary analysis regarding a conceptual configuration where a simple iterative and manual sizing of the thermal protection system has been successfully performed, without adopting any optimization procedure.

Data Availability
The data supporting the analysis are from previously reported studies and datasets, which have been cited in Ref. [20,21].Specifically, the multidisciplinary analysis discussed in Section 7 has been performed on an RLV shape formulating a multiobjective optimization problem.The problem set-up, including the range of design variables, has been detailed in [21].Thermal properties of materials used in the TPS sizing procedure are available in [27].

4. 1 .Figure 1 :
Figure 1: Schematic representation of a stick obtained using point source blobs with centers placed on the linear segment of length l (a); strength field generated by self-blending property (b).

Figure 2 :
Figure 2: Stick primitives obtained with n blob = 6 and 20: constant radius (a, b); variable radius (c, d).The stick support becomes more regular increasing n blob ; the strength field remains bounded to the unit value.

Figure 4 :
Figure4: Arbitrary stick distribution with a longitudinal gradient onto the B-grid adopted for thermal protection system modelling.

8. 5 .
Results.The multidisciplinary analysis previously described has been applied to perform sizing and material distribution of a TPS for an RLV with dimensions chosen equal to length l tot = 8 6 m, wingspan w s = 8 5 m, and cabin height h = 1 5 m.Thickness distribution of Figures8(a) and 8(b) is obtained by setting n stick = 3.

Table 1 :
Parameters adopted in the modelling of TPS configurations of Figures5 and 6.

Table 2 :
Waypoints adopted for aerodynamic computation.

Table 3 :
Selected impact methods for aerodynamic computations.the variable δ i vanishes for negative values of the difference T max,iwi − T max,wi where