Development and Application of Computational Tool Using Local Surface Inclination Methods for Preliminary Analysis of Hypersonic Vehicles

: This work presents a computational tool for preliminary analysis of hypersonic vehicles, based on local surface inclination methods: the HipeX. This program was developed for reading standard triangulation language (STL) geometry files and calculating pressure coefficient and temperature distributions over vehicle’s surface using the Newtonian, modified Newtonian or tangent-wedge methods. Validations were made with a cylinder and a sphere, where only the Newtonian method was applied, and with experimental data from Apollo capsule at Mach 10, where the Newtonian and the modified Newtonian methods were applied. These validations presented the code capability to read geometries as well as to estimate aerodynamic force coefficients. A preliminary application was to predict the aerodynamic force coefficients of a generic hypersonic vehicle over constant dynamic pressure trajectories of 23,940, 60,000 and 95,760 N/m 2 with zero angle of attack. With a fixed dynamic pressure of 60,000 N/m 2 , this vehicle was tested over several Mach numbers and with angle of attack variation from –10 to 10 deg. Zero angle of attack investigation showed minor changes on the force coefficients with altitude, while the variation of angle of attack produced more pronounced variations on these parameters. Maximum flow temperatures over vehicle’s surface were estimated ranging from 850 to 5,315 K.


INTRODUCTION
molecules, dissociation and ionization. Also, the very high temperature inside the boundary-layer makes the viscous flow analysis a very important issue (Gulli et al. 2012;Pish et al. 2019). Indeed, several Navier-Stokes solvers can be found elsewhere such as LAURA (Cheatwood and Gnoffo 1996), the software VULCAN for scramjet applications (White and Morrison 1999), FUN3D (Anderson and Bonhaus 1994) and the ANSYS Fluent (ANSYS 2019). Despite that, rapid estimates for pressure distribution can be obtained using local surface inclination methods (Anderson 1989). So far, the HipeX tool has the capability of using the following methods: Newtonian, modified Newtonian, and tangent-wedge. Newtonian and modified Newtonian methods are most used for blunted body geometries (where normal shockwaves are established) and tangent-wedge method is most used for slender geometries (where plane oblique shockwaves are established).
In this paper, it is discussed how the geometry is imported, it is given a short explanation of the applied methods, how force coefficients are calculated, the allowed strategies for trajectories and software interfaces. Moreover, results from the code validation with simple geometries as well as with Apollo capsule reentry data are shown. Finally, a generic hypersonic vehicle is tested over prescribed trajectories.

METHODOLOGY
The HipeX computer program consists of three interfaces: the geometry interface, a trajectory interface, if the user decides to use it, and an interface for the methods.

GEOMETRY READING
Geometry files in the standard triangulation language (STL) format were chosen for this work due to its simplicity to represent surfaces and because it is widely used on computer-aided design (CAD) tools. In such format, the solid can be represented in two different ways: binary and ASCII. Regardless the choice, the STL file format divides the solid surface into various triangles, storing the positions of the vertices and the components of the normal vectors. Fig.s 1 and 2 show the typical configurations of STL file format for ASCII and binary, respectively. All geometries were developed using Autodesk Inventor (Autodesk 2018) as well as the STL files. Regardless the choice, the program is able to separate the coordinates of each vertex of the triangle and the components of the normal vector. After this storage, the program uses this data to calculate the volume, the area of each triangle created in the STL tessellation, total body area, vertex and face quantities.
The volume of the body, enclosed by the surface, is calculated using the divergent theorem. So, it is calculated through the relation: where N is the total number of triangles and A, B, and C are vertices of the i-th triangle, with: where n i is the vector normal to the surface of the i-th triangle and the area of this triangle is given by: From this relation it is also possible to obtain the total surface area through:

LOCAL SURFACE INCLINATION METHODS
Local surface inclination methods are in general based on simplifications of the flow field. In such methods, there is no need for a detailed calculation of the flow around the vehicle since the local surface inclination is enough to provide the distribution of the pressure coefficient.
Three methods were utilized in this work: Newtonian, modified Newtonian and tangent-wedge (Anderson 1989). Each of them is able to calculate the pressure coefficient of the body surface. The pressure coefficient is related with the freestream pressure P 0 , velocity V 0 , and density ρ 0 as follows: where p is the local pressure. Basically, these methods depend on the flow deflection angle θ, which can be determined from the vector operation (see also Fig. 3): Negative flow deflection angles (θ <0) correspond to regions inside the shadow of the incident flow, where it is assumed that C p =0. In the Newtonian method, the simplest one, but applicable to two-or three-dimensional bodies in hypersonic regimes, the pressure coefficient is given by: The modified Newtonian method is also applicable for two-or three-dimensional bodies in hypersonic flight, but is preferable for blunt-nosed bodies (Anderson 1989). It corrects Newtonian method based on a better estimate of the maximum value of C p (C pMAX ) found at the stagnation point. So, the local pressure coefficient is expressed by: C pMAX is calculated considering that the maximum pressure is achieved when the freestream flow is processed through a normal shock wave and through an isentropic compression, this process corresponds to the Rayleigh pitot tube theory, so that (Anderson 1989): where Y 0 is the ratio of specific heats, assumed to be constant, and M 0 is the flight Mach number.
The accuracy of the Newtonian and modified Newtonian methods tend to improve as the flight Mach number increases as pointed out by Anderson (1989).
The third method which was adopted is the tangent-wedge. It is applicable to slender bodies with two-dimensional shapes on supersonic or hypersonic flow regimes, and obviously extensible for three-dimensional shapes with planar shockwaves. The method assumes that the pressure at the surface can be obtained from a tangent wedge with an equivalent flow deflection angle, by doing such assumption; it is possible to use the oblique shockwave relationships. So, this method is capable of calculating the pressure coefficient as well as other thermodynamic properties (like flow temperature and density).
In fact, it is possible to find the shockwave angle β iteratively from the relation (Anderson 2001): Furthermore, the density ratio across the shockwave can be expressed by: The pressure ratio is expressed by: And, the temperature ratio is obtained via:

CALCULATION OF TOTAL FORCE DUE TO PRESSURE
Once C P is calculated through one of the local surface inclination methods described earlier; the total force due to pressure can be determined by the sum of the contributions of each surface element: The force coefficients on the directions X, Y and Z can be found by: where S ref is the reference area.
Knowledge of coordinate systems is a base point for proper calculation of forces. In particular, the body and wind coordinate systems were used for these simulations (Fig. 4), where it was adopted the same definitions as Zipfel (2007). The body coordinate system is defined by x-axis positive towards the vehicle nose and z-axis is positive down, y-axis is found with the right-hand rule. In the wind coordinate system, the x-axis is defined as parallel and in the direction of the relative wind velocity of the vehicle (with magnitude V 0 ), the wind z-axis is positive down, the y-axis is obtained from the right-hand rule.
For given angles of attack (α) and sideslip (β), also shown in Fig. 4, the transformation matrix to interchange vectors from the body to the wind coordinate system is: On the other hand, to transform vectors from the wind to the body coordinate system, one can use: where, one can see the following properties: and, Since most of calculations in the code developed is made in the body frame, it is convenient to transform a given flow velocity to the body coordinate system: Moreover, since the aerodynamic forces are directly related to the forces on the wind coordinate system, one can obtain the drag, lateral force and lift coefficients by applying:

TRAJECTORY STRATEGIES
Three trajectories were chosen for a possible aerospace flight at hypersonic regime: constant dynamic pressure; constant Reynolds number and constant Mach number.
Constant dynamic pressure trajectories work as base for understanding the load limits, very important for the structural design, since all aerodynamic coefficients are proportional to this quantity. The dynamic pressure (q 0 ) is defined by the atmospheric density and flight velocity on a given altitude as in the equation: The choice of constant Reynolds number can be justified by the fact that laminar or turbulent flow regimes can be inferred by this parameter. In fact, there are many other variables involved in the laminar to turbulent flow transitions; however, as a preliminary estimate, it is a common practice to consider the critical Reynolds number as the sole parameter involved on the transition prediction. Although a precise value of this quantity is generally obtained from experiments, Heiser et al. (1994) indicates 10 7 as a good estimate for the beginning of the turbulent flow. The Reynolds number based on freestream properties is defined by: where, L is a reference length and μ 0 is the dynamic viscosity of air, calculated using the Sutherland's law (Anderson 1989). Constant Mach number trajectories are also relevant, and its behavior can be linked to the variations of the atmospheric temperature by means of its dependence on the speed of sound (a 0 ): Freestream thermodynamic properties were calculated using COESA 1976 atmospheric model (NASA 1976), considering geometric altitude.

SOFTWARE INTERFACES AND FLOWCHARTS
To read the geometry of a solid, the HipeX uses an interface created as seen in Fig. 5. The program starts with the chosen solid, followed by the selection of units. To improve the file compatibility, rotation angles can be supplied by the user so that the solid fits to the coordinate system desired for the simulations. After the reading process, the program uses the acquired data to calculate the volume, the area of each triangle, the total area of the solid's surface, and the numbers of vertices and faces. To check out the importing process, the solid can also be inspected. The flowchart for the reading geometry program is shown in Fig. 6. After reading the geometry, the interface for the methods could be executed. Figure 7 shows the interface for the methods created where the user can choose which method will be used. Also, it is possible to choose if it is a single or a multiple point evaluation. If it is a single point evaluation, the user must provide the flight velocity, angle of attack, and angle of sideslip, reference length and area, and altitude. For a multiple trajectory points, the area and the length must be followed by a trajectory file in XLSX format. If the user decides, the trajectory file can be obtained from the trajectory interface. The flowchart for the surface inclination methods can be seen in Fig. 8.
The interface for trajectory calculations can be seen in Fig. 9, where initial and final altitudes must be provided by the user; the number of discretization points and reference length are also needed as inputs (see Fig. 10 for the flowchart). As stated before, the software gives three possibilities for trajectories: constant dynamic pressure; constant Reynolds number and constant Mach number. The user can choose the most appropriated trajectory and insert the value for the constant. Then, the file is saved in an XLSX file. To check out the calculation, plots can also be made by choosing the most pertinent properties for the user.

RESULTS AND DISCUSSIONS
Parameters for several geometries were evaluated to ensure that the code can be applied in real scenarios. Firstly, elementary geometric bodies with known analytical solutions were used to validate both geometry reading and Newtonian method application processes. Secondly, analysis of Apollo aerodynamic coefficients was conducted to compare the evaluations of HipeX with experimental data. Finally, as application, an arbitrary aerospace vehicle was tested over a prescribed trajectory at zero angle of attack, and also tested for several angles of attack at different Mach numbers.

EVALUATIONS OF SIMPLE GEOMETRIES
It was possible to compare the total area and volume found by the program with the total area and volume from analytical formulas for a cylinder and a sphere (Fig. 11). The results show that the built functions provide values that are very close to the exact values (Tables 1 and 2). The values of the axial force coefficient C X are compared with the expected values provided by Anderson (1989), only for the Newtonian method. The results are shown in Table 3. The Newtonian method was chosen due to its independence on flight parameters. Again, very good agreement was observed. It is also important to keep in mind that these results depend on the number of vertices and triangles coming from the exporting process used in the CAD. The cylinder was represented by 285 vertices and a total of 566 elements, in binary format, 1.0 m of diameter and 1.0 m high (Fig. 11a). The sphere was composed with 4,514 vertices, with a total of 9,024 elements, in binary format, 1.0 m diameter (Fig. 11b).

APOLLO CAPSULE REENTRY
In this study, the purpose was to observe the behavior of the force coefficients in the X and Z directions while the angle of attack was varied from 0 to -30 deg., at Mach 10. The Apollo capsule geometry data was taken from Moss et al. (2006) (Fig. 12) and it was represented by 4,560 vertices, 9,116 elements with a reference area of 12.02 m 2 , equivalent to the capsule's maximum cross-section area. A typical run time was less than 200 s for the evaluations over the entire studied range, using an i7-2.20 GHz processor with 8 GB of RAM.
The results are seen in Figs. 13 and 14. The results were compared with those presented by Hirschel and Weiland (2009), where most of the aerodynamic data were obtained from wind tunnel tests. Considering the coefficient of axial force (Fig. 13), the modified Newtonian method presented better results in general. Also, it can be noted that for lower angles of attack the estimate of C X improves, in fact, the relative deviation when compared to the results of Hirschel and Weiland (2009) is less than 10% when AOA = 0 deg. and less than 5% for AOA = -30 deg., considering the modified Newtonian method. In the case of the normal force coefficient (Fig. 14), the behavior was different, the modified Newtonian method produced a greater relative deviation for lower values of AOA, for example, at AOA = -30 deg. a relative deviation less than 41% was found while the Newtonian method showed a relative deviation of less than 36%. However, the deviations of these two methods decreased notably when the angle of attack approached the null value.
The use of the tangent-wedge method was discarded due to the bow shock expected for the Apollo capsule. This method is more appropriate for attached shockwaves, present on slender geometries. To better assess the code capability, the distribution of C p according to the Newtonian method for the values of angle of attack 0, -10, -20 and -30 deg. are shown in Fig. 15. For a single point calculation, the run time for this geometry was less than 21 s, including the time for saving results in a spreadsheet file. For the angle of attack 0 deg., the maximum value of C p =2 was reached at the nose center, while in regions located at shadow regions it assumed a null value. When the angle of attack varied, the position of maximum C p transited away from center.

APPLICATION: GENERIC HYPERSONIC VEHICLE
As an application, an aerospace vehicle, shown in Figs. 16 and 17, was designed to test the software capability. The vehicle geometry is formed by a combination of intersecting planes, with sharp edges and faceted surfaces intended to reduce the design and manufacturing costs, with a similar design philosophy presented by Boehrk (2014). The model has 1.0 m of length, 0.894 m of wing span, with a planform area of 0.5 m². The upper ramp has a 10-deg. slope while the lower ramp has a compression angle of 19 deg. The rudder and elevons have symmetric diamond profiles. The STL mesh was composed of 84 elements and 44 vertices. Firstly, flights with constant dynamic pressures of 23,940, 60,000 and 95,760 N/m 2 with zero angle of attack were simulated. The reasons for these choices were based on the suggested hypersonic airbreathing corridor in terms of dynamic pressure (between 500 and 2,000 lbf/ft 2 ) (Heiser et al. 1994), 60,000 N/m 2 (1,253 lbf/ft 2 ) is an intermediate value. Secondly, the same geometry was tested for various Mach numbers, but keeping a constant dynamic pressure of 60,000 N/m 2 , in which the angle of attack ranged from -10 to 10 deg. The altitude was varied from 30 to 50 km, this altitude range would mimic a typical hypersonic test flight; the freestream properties were obtained from the trajectory interface. Mach number and velocity varied as depicted in Fig. 18(a) and (b), respectively. The test matrices are shown in Tables 4 and 5. Due to the fact that the chosen geometry only generates attached planar shocks, the most appropriated method was the tangent-wedge, which provided the pressure distribution and consequently the forces applied on the vehicle surface. Moreover, it was possible to obtain the inviscid flow temperature at regions adjacent to the vehicle.   60,000 9 -10 -+10 60,000 12 -10 -+10 60,000 15 -10 -+10 60,000 25 -10 -+10 For each of the simulated configurations of Tables 4 and 5, the computation time was less than 47 s, with a 1000 points discretization.
The variation of the force coefficient on the X-direction with the altitude can be seen in Fig. 19, the choice of this coefficient was supported by the fact that the drag coefficient is closely related to it. As the altitude varied from 30 to 50 km, the C X varied from -0.0243 to -0.0192, for the case of a dynamic pressure of 23,940 N/m 2 , this was the greatest variation. The dynamic pressure of 60,000 N/m 2 produced a variation from -0.0213 to -0.0190 while in the case of 95,760 N/m 2 one can see that it varied from -0.0204 to -0.0189. From one to another dynamic pressure, it can be seen that the difference was more pronounced at low altitudes, corresponding to lower Mach numbers (Fig. 18a). As the altitude increases, the Mach number increases and the axial force coefficient becomes less sensitive to the Mach number variation. This fact is in accordance with the hypersonic Mach number independence principle, brought by a very thin shock layer, which is a result from the decrease of the shock wave angle with Mach number. To further investigate the aerospace vehicle behavior, plotted in Fig. 20 are the values of the force coefficient on Z-direction, C Z , for the three cases of dynamic pressures, as the altitude increases. In this parameter, as in the axial force coefficient, lower Mach numbers corresponded to a higher dependence on the altitude and also on the dynamic pressure. For the case of 23,940 N/m 2 , this parameter varied from -0.0346 to -0.0321. The intermediate case, 60,000 N/m 2 , produced a variation from -0.0328 to -0.0321. Meanwhile, the highest dynamic pressure of 95,760 N/m 2 corresponded to a variation from -0.0324 to -0.0320. Again, as the altitude increases, the normal force coefficient becomes less sensitive to the altitude and dynamic pressure variations. Next, the variation of C X as the angle of attack changes is shown in Fig. 21 for Mach numbers 9, 12, 15 and 25, all corresponding to a constant dynamic pressure of 60,000 N/m 2 . A large variation of this parameter can be seen as the angle of attack was varied. The largest variation was notable seen for positive angles. As the angle of attack ranged from -10 to +10 deg., the axial force coefficient varied from -0.0163 to -0.0382 in the case of Mach number 9. In a similar fashion, at a higher Mach number of 25, C X varied from -0.0146 to -0.0367, almost doubled. The curve behavior was slight nonlinear for negative values of angle of attack, meanwhile, for positive values, the variation tended to be more linear. This effect could be explained by the larger angle at the lower surface of the vehicle which causes higher pressures, a 19-deg. ramp, while at the upper surface, the ramp has a 10-deg. inclination, producing mild modifications as the flow is turned negatively.   Figure 23 shows the contours of pressure coefficient and static temperature on the fluid adjacent to the vehicle surface at Mach numbers of 9 and 25, both on a constant dynamic pressure trajectory of 60,000 N/m 2 with zero angle of attack. In Fig. 23(a), corresponding to the Mach number 9, the maximum value of the pressure coefficient was 0.2763 at the lower compression ramp; this region also was related with the maximum pressure coefficient at Mach 25, as seen in Fig. 23(b), assuming 0.2584 in this case, a decrease of less than 7%. In the regions of shadow, corresponding to zero or negative deflections of the flow with respect to the freestream, C p is null.  Knowledge of the vehicle's thermal loads is of paramount importance to a preliminary design, so that surface materials, cooling systems, or even reactive heat protection systems would be designed to withstand the high heat transfer rates associated to the hypersonic flight. Temperature estimates were given by Eq. 15. Indeed,Fig. 23(c) shows that flow temperature could be very high, about 850 K, in the case of Mach number 9, while in the case of Mach number 25, the maximum temperature reached on the flow was as high as 5,315 K (Fig. 23d). Moreover, it can be seen that the largest thermal load was in the 19-deg. compression ramp, in both cases. Also, the method predicts that temperatures drop severely when the flow is parallel to freestream flow, corresponding to shadow regions. It is important to notice, however, that it is expected an over prediction of temperatures since no high-temperature effects were taken into account, for instance, excitation of the vibrational mode of the diatomic molecules and dissociation tend to diminish the flow temperature. Also, if boundary layer effects were accounted, considerable changes on the flow temperature around the vehicle could be found.

CONCLUSIONS
The presented work provides a description of the HipeX code, which is an alternative with a relatively low computational cost, based on local surface inclination methods, for aerothermodynamic analysis of hypersonic vehicles. Firstly, it was discussed how the program reads geometry files in STL format, how the surface inclination methods (Newtonian, modified Newtonian and tangent-wedge) are used to obtain the pressure coefficient, how to calculate the total pressure force acting on vehicle's surface, what are the trajectories that could be adopted (constant dynamic pressure, constant Reynolds number or constant Mach number), and the software interfaces.
A validation of the methodology was made on two ways: in the first one, simple geometries with known analytical results for the Newtonian method were considered. Secondly, comparisons were made with experimental data from Apollo capsule and Newtonian and modified Newtonian methods.
As a primary application, a generic hypersonic aircraft was designed and tested over constant pressure trajectories of 23,940, 60,000 and 95,760 N/m 2 with the tangent-wedge method, firstly with zero angle of attack and secondly with a mild variation of this parameter on selected Mach numbers. The results were focused on the values for the force coefficients. Although during the zero angle of attack investigation the force coefficients showed only minor changes, the angle of attack variation produced considerable variations on these parameters. Moreover, it was possible to investigate how the flow temperature was distributed over the vehicle's surface. Temperatures of the inviscid flow tangential to the vehicle surface were also evaluated. However, high fidelity simulations are indispensable to better predict the thermal state of the surface. Indeed, it should be addressed that the surface temperature is dependent on several parameters not discussed in this work.
Finally, the importance of preliminary tools for aerospace vehicle design relies on the possibility of screening configurations that would best fit the proposed mission, knowledge of advantages or disadvantages can help further development of the selected vehicle configurations. Complex high-fidelity models are of paramount importance to understand the aerodynamics of such vehicles; however, it is possible to estimate many aerodynamic performance parameters with local surface inclination methods with a fraction of the computational cost of a more detailed analysis.

ACKNOWLEDGMENTS
Editors and authors are thankful to Fundação Conrado Wessel for providing the financial support for publishing this article.