Computational Aerodynamic Optimisation of a 2D Space-Shuttle Geometry

The need for cost minimisation in astronautical engineering, to help make human kind an interplanetary species, is not just limited to the flight and manufacturing of spacecraft. Hypersonic aerodynamic optimisation has been an expensive and time-consuming necessity since the beginning of the Apollo missions. A Boltzmann-BGK solver developed at Swansea University is utilised in this research. Firstly, a study was undertaken to determine how the results of an arbitrary simulation would vary depending on the number of nodes in the velocity space mesh, with the higher nodal number producing a more realistic temperature distribution. The best nodal configuration found was then used to find the best nose geometry for a space shuttle during hypersonic re-entry and ascent. An extension on the geometric variables (an increase in the space-shuttle length) was then simulated to evaluate how this changed drag results compared to an identical nose geometry with a shorter length; the results were found to be identical. Additionally, a study was undertaken to evaluate the effect of extending the outer boundary of the simulation from a half boundary to a full boundary encompassing the entire geometry, focusing on the temperature results, with a focus on the rear of the space shuttle. The full boundary extension was found to only change the temperature distribution at the rear of the geometry, the front portion of the geometry had an identical temperature distribution to the previous half-boundary simulation. Finally, the rear of the fully extended boundary simulation was altered to have a rounded rear, this was found to reduce the temperature distribution protruding from the rear.

portion of the aircraft. Due to the reliance on wind tunnels in the sector, there are many difficulties in the hypersonic aerodynamic optimisation process. The necessity to create models to test and then recreate and retest can be very time consuming and expensive, the structural damage caused by overheating is a constant risk and the limited availability of hypersonic wind tunnels means smaller companies will have little chance to test [4,5]. There is a clear opportunity for the development of computational fluid dynamics to remove virtually all the negatives of hypersonic aerodynamic optimisation.
Computational aerodynamic optimisation practitioners have already developed effective programming software based on either Navier-Stokes [6,7] or Euler equations [8]. Due to the molecular variables of a high speed and altitude fluid differing from continuum molecular variables, the use of a more fundamental physical description must be adopted. The Boltzmann equation has been utilised in several research studies [9][10][11][12][13][14][15][16][17], all highlighting the potential of a novel Boltzmann solver due to its accurate use in both rarefied atmospheric conditions at hypersonic speed and more typical, low-speed atmospheric conditions. However, to the authors knowledge, no other Boltzmann-BGK based hypersonic CFD solver has been utilised for space shuttle geometry optimisation except that of Evans, et al. [9,10].
More complete studies on space shuttle aerodynamics during hypersonic flight have been completed successfully in multiple research papers. Studies such as [18] and a similar, more recent study [19], both compare hypersonic CFD data to flight data, whilst these studies are specified for comparison, they do not utilise the Boltzmann equation, hence the novelty of this paper.
Using the Finite Element Boltzmann-BGK solver developed by Evans, et al. [9], an optimised 2D space shuttle geometry will be found by simulating many different space shuttle geometries, all a 2D 'double ellipse' shape. A geometry which minimises drag the most on ascent will be deduced, as will a geometry that minimises temperature on descent. A compromise in the geometry will be found to find the most efficient overall geometry. The program utilises the discontinuous Taylor-Galerkin finite element method of discretisation combined with a novel optimisation approach.

Continuum equation issues
The degree to which a gas can be considered rarefied can be expressed through the value of its Knudsen number, Kn, This is a non-dimensional quantity defined as Where L is any characteristic dimension in the flow and λ is the mean free path of the molecules, i.e. the average distance travelled between collisions.
The typical assumption that the Navier-Stokes equations are no longer valid if the Knudsen number is greater than 0.1 [20] is misleading if any overall dimension of the flow L has been chosen. The limit of the continuum equations should therefore be specified using a local Knudsen number where L is defined as Where L is now the scale length of the macroscopic gradients [20] and ρ is the density of the fluid. A diagram showing the appropriateness of a range of equations at various Knudsen number can be found in Figure 1.
Where the only difference between equations (3) and (4) are the β and π powers.

Boltzmann-BGK equation
The simplification made by Bhatnager, Gross and Krook (BGK) in the Boltzmann equation assumes that the departures of the flow from thermodynamic equilibrium are small. A further assumption deduced from their conclusion is that the effect of the molecular collisions in a non-equilibrium fluid, forces the non-equilibrium molecular velocity distribution back toward a state of equilibrium, at a rate proportional to the molecular collision frequency [10]. The Boltzmann-BGK equation may be written as Where n is the molecular number density, F takes into account any force fields that are present, f = f(r,c,t) is the molecular distribution function of the fluid, (r,t) is a mathematical term proportional to the collision frequency of the molecules, r represents the position of the molecules and f 0 is the local Maxwellian equilibrium distribution function discussed above.
A direct comparison between a Boltzmann-BGK scheme and a continuum scheme can be found in [27] where the Boltzmann-BGK solver is found to be more accurate at modelling shock tube problems.

Boltzmann-BGK Solver
The reader will now be given a brief explanation of some of the key aspects of the Boltzmann-BGK solver developed by Evans, for a more detailed explanation, see [9,10].

Physical space discretisation
The Boltzmann-BGK solver is restricted to two physical space dimension problems. The physical space is referred to as p-space domain, Ω r . The p-space domain is discretised into a linear, discontinuous, triangle element, unstructured mesh ( Figure 2). There are multiple advantages to using a discretisation technique of this nature [9,10] but it is specifically utilised for rarefied hypersonic flow because of the shock capturing properties. For further discretisation techniques, the reader is directed to [28,29].
A diagram showing the appropriateness of a range of equations at various Knudsen number can be found in Figure 1.
Viewing Figure 1, it is clear to the reader that a computation utilising the continuum model will not be appropriate for a solver in rarefied atmospheric conditions. Continuum theory always assumes that a gas is in thermodynamic equilibrium, i.e. in thermal, chemical and mechanical equilibrium, however in this paper, the gas analysed will not have the appropriate time needed to return to equilibrium and so a separate set of equations must be utilised. The inaccuracies of using nonequilibrium Navier-Stokes equations can be found in [21].

Boltzmann equation
A detailed description of the Boltzmann equation theory and application can be found in [25,26]. The reader is again directed to Figure 1, where the accuracy of using a Boltzmann type solver is highlighted graphically. The Boltzmann equation accurately describes a fluid that is not in a state of equilibrium [23].

Maxwell distribution function
At this time, it is appropriate to introduce the Maxwellian distribution function to the reader. This function will be used in the Boltzmann-BGK solver as it allows the equation to still be a non-linear integro-differential equation [9] whilst applying the BGK approximation to simplify the computation.
The Maxwellian distribution form applied to the Boltzmann-BGK solver is Where c 0 is the bulk velocity of the flow, c is the velocity of the flow and T equals the temperature of gas (measured in Kelvin), m is the molecular mass, R is the molar gas constant and k is the Boltzmann constant [9]. Further detail about this function can be found in [9,26].

Velocity space discretisation
The velocity space is referred to as the v-space domain, Ω c , and theoretically is unlimited in extent, however to create a working solution to the solver, a limit on the maximum velocity of molecules is identified; any molecule travelling faster than this is assumed to have negligible impact on the solver accuracy. To apply this principle to the Boltzmann-BGK solution, an estimate of the mean thermal molecular velocity must be made by using moment calculations (see [9,10] for full mathematical description). The relationship between the thermal molecular velocity and other flow variables is found to be c′ is the mean thermal molecular velocity. The finite limit of the radial extent of the velocity domain, r v , can be viewed in Figure 3. Step 1) The piecewise-constant increment is calculated at every physical space element, using the equation Where Δ(nf) re,c is the piecewise-constant increment, ∆t is the timestep (global), k is the summation that extends over three nodes of the element re, N k is the finite element shape function at node k (piecewise-linear) and Element fluxes during a timestep of a half are calculated using the piecewise linear discontinuous equation Step 2 The last step involves a piecewise-linear approximation for ∆(nf) at every physical space element. All element node values of the solution increment over the complete timestep [10] are calculated using the equation A spectral v-space domain is a single, high order element [9] and so gives for more efficient calculations over the full domain. To apply this domain, it must first be mapped from Cartesian or polar coordinates (in real space), to a quadrilateral element in the (η, ζ) plane.
View Figure 4 to see the mapping of the quadrilateral element. The mathematics of this approach can be found in [9,10] but a detailed understanding of this method is not necessary to the reader. This method creates a symmetric distribution of points, with no preference on a radial direction.

Taylor-Galerkin discretisation method
A summarised explanation of the Taylor-Galerkin discretisation will now be highlighted, for a more in-depth description, see [9,10,30].
An approximation is used in this method Where f is the distribution function, n is the

Boundary conditions
Three key boundary conditions have been applied to the Boltzmann-BGK solver: Inflow, outflow and wall. An explanation of each boundary condition and their relevant equations can be found in [9,10,30]. These boundary conditions are applied accurately to the solver by modelling of 1 2 , m n c F + , a part of equation (11). A brief explanation of each boundary condition will now be given.
Inflow: All gas properties necessary to calculate the Maxwellian distribution function in two dimensions, see equation (4), are present. It is now possible to calculate the inter-element flux at the inflow using equation For molecular velocity directed into the domain. For molecular velocity directed out of the domain, equation (9) is utilised.
See Figure 5 for a visual representation of the inflow boundary condition, where all of the gas properties needed for calculations are represented.

Outflow:
The inter-element calculations at the outflow boundary domain is more simplistic then at the inflow, as it is not dependant on the direction of molecular velocity.
The inter-element flux equation is ( ) ( )  Figure 5 shows a visual representation of the outflow boundary condition.

Wall:
The wall condition means that the value for mass flux across the boundary is zero [9,10,30]. This is represented in the equation Where Γ is the physical space domain boundary [9,10,30] and F n,c = (c.n) f(c,r,t).

Parallelisation
For the Boltzmann-BGK solver to be applied to complex problems and work in a reasonable timeframe, the use of parallelisation in the code is necessary. The method chosen is physical space domain decomposition.
The information communicated between each processor is that of the inter-element flux. It is not-ed that the flux information can only be communi-cated in one direction and never both directions at the same time, for example, processor 1 will pass its inter-element flux information to processor 2 (left to right) but then processor 2 can only pass on its inter-element flux information to processor 3 (left to right, not right to left). The direction that the information travels is dependent on the molec-ular velocity vector.
For a more complete description of the parallelisation process used, see [30].

Software and Hardware
The solver developed by Evans is coded in the FORTRAN 90 language, this language was

V-Space Node Value Study
A study was undertaken to determine how the results of an arbitrary simulation will vary depending on the number of nodes in the v-space mesh. Note that the geometry in this study is geometry 1 chosen due to the allocable memory and speed of computation [33]. In this authors research, the solver was compiled and run using a Unix cluster account [34] and utilised the power of the Swansea University High Performance Computing (HPC) facility [35]. The alterations to the space-shuttle geometry were all completed using Matlab software due to its simplicity at plotting two dimensional shapes [36]. Figure 6 shows a diagram of the space shuttle geometry variables that were altered. Values of a, b, c and d were changed during each new geometry simulation.    The v-space mesh was then changed to contain a 120 × 120 mesh, i.e. 14400 nodes in total. Figure 8 shows the temperature distribution for the Lobatto 120 simulation.

Results
The lower v-space mesh simulation needed 1000 more timesteps than the higher v-space mesh simulation to gain results with low variation of lift and drag (view Figure 9 for the drag variation in Table 1. See Appendix for a diagram of points in the velocity space discretisation.

Lobatto 80 Quadrature
This v-space mesh contained a number of nodes equal to 80 × 80, i.e. 6400 nodes. Figure 7 shows the temperature distribution for the Lobatto 80 simulation.  Table 1 that the simulation produces, and whether the simulation has converged to produce accurate results (produces a flat line, see Figure 9 and Figure   10 for examples of this).

Lobatto 120 Quadrature
The data for geometries 1 to 4 are presented in Figure 11 and Figure 12. These simulations have been produced using the conditions stated in Table   2. Figure 11 and Figure 12 show that all geometries produce results for lift ( Figure 12) and drag ( Figure   11) to steady state (minimum 3 significant figures), this is why all lines are flat by 5000 timesteps. This suggests that all results are reliable.

Aims of hypersonic re-entry simulation
This simulation was utilised to evaluate the 4 geometries specified in Table 1. The aim of this and Figure 10 for the lift distribution). However, it is clear from these figures that the variations in both drag and lift results are small (less than 0.01 difference) and so the necessary increase in time needed to use a Lobatto 120 × 120 v-space mesh in a simulation may not be ideal when considering the relative accuracy of the Lobatto 80 × 80 results (which takes far less time to run, around 24 hours less).
Comparing Figure 7 and Figure 8, it is clear that the temperature variations have similar values (the reader is directed to the bottom-right corner of both figures). The Lobatto 120 × 120 simulation has a more condensed temperature distribution then the Lobatto 80 × 80 simulation, this is due to the higher density of v-space meshing which has produced a more realistic temperature distribution. Table 2 shows the variables chosen for all simulations in this research. These values were chosen as these conditions are roughly the most extreme conditions a space shuttle will experience during its product life.

Data
Data is used to evaluate two things: the drag  Table 1) with varying V-Space mesh distributions. Note, Lift (m 2 ) = Lift (N) ÷ dynamic pressure (N/m 2 ). that the geometry experiences needs to be minimised, particularly during hypersonic re-entry due to the large values of temperature produced at simulation was to find the best geometry of the 4 during hypersonic re-entry.  Table 1) has the most success in minimising temperature under identical conditions. All figures have the same temperature scale (45,000 Kelvin is red, 0 Kelvin is dark blue) and Figure 7 has by far the least 45,000 Kelvin temperature regions.

Temperature minimisation: The temperature
Geometry 1 also is the best geometry evaluated in terms of lift minimisation. Viewing Figure 12, geometry 1 flattens out at a steady state solution for lift at a lower value than all other geometries. Note that geometry 2 is the next best geometry for lift minimisation. Figure 16 shows the pressure distribution of the this speed when reentering the Earth's atmosphere. The large temperature can cause structural damage and even failure in some severe cases.
Lift minimisation: Lift needs to be minimised during re-entry. More lift will make it more difficult for a spacecraft to successfully re-enter the atmosphere and come to a controlled landing.

Results
Note that Figure 7 shows the temperature distribution for geometry 1 ( Table 1). Figure   7, Figure 13, Figure 14, and Figure 15, it is clear  Drag minimisation: It is necessary to minimise drag during the ascent phase of a space mission, as this will minimise the amount of thrust needed to reach space. Less drag also means a smaller net force acting upon the structure, decreasing the likeliness of failure of the structure. Figure 11 shows that geometry 1 ( Table 1) is once again the best geometry. It creates the lowest value of drag, 0.036 m 2 . Note, geometry 4 is the second best at minimising drag. best geometry for hypersonic re-entry (geometry 1).

Results
Note, the striations in the temperature and pressure figures are caused by a numerical anomaly produced when using high order spectral discretisation of the velocity space. This is an area of ongoing research in the field of high speed computational aerodynamics.

Aims of ascent simulation
This simulation was again evaluating the 4 geometries specified in Table 1. The aim of this    Figure 17.

drag. A visual representation of this can be seen in
The value of e was 6 (corresponding to Figure  18 and Figure 19). The reader should note that all previous simulations were limited to an e of 1.
The drag results for this body extension simulation are then compared to Figure 7 (same nose geometry, a back extension of 1).
Only the drag results are considered due to issues with the lift not converging in the allotted timesteps. However, the drag converges quickly It is clear from Table 3 that geometry 1 is the ideal geometry for use at any stage in a spaceshuttles flight; therefore this geometry will now be evaluated further by expanding geometric parameters.

Simulation 3 -Increased Shuttle Geometry Parameters Body extension
Geometry 1 variables were then used to form the front nose geometry, however, the body of the geometry was extended to evaluate its effect on   Figure 7. This is the only difference between the full boundary and half boundary temperature distributions, both have an identical distribution when considering only the temperature effect on the front portion of the geometry.
Viewing Figure 21, this geometry has a rounded rear to its design. Comparing this to Figure 20, the rounded rear of Figure 21 causes the temperature at the rear of the geometry to significantly decrease in both magnitude and area.
The rounded rear of Figure 21 also influences the temperature distribution at the front of the geometry. Much like when comparing v-space nodal values, the distribution values are virtually identical, but Figure 21 has a more condensed distribution then that of Figure 20.

Conclusion
This research has outlined the potential of hypersonic aerodynamic analysis using a novel Boltzmann-BGK solver. Issues with the solver have been highlighted, these include: being computationally expensive to run and a numerical anomaly produced when applying high order spectral discretisation of the velocity space.
Issues aside, the power of this solver has been presented by simulating a number of scenarios that are applicable to a space shuttle mission. The success of this solver under different scenarios paves the way for further development in hypersonic Boltzmann based solvers with the hope to a steady state solution in 5000 timesteps and so can be compared to the geometry 1 drag results ( Figure 18).
Results-body extension: Viewing Figure 18, it is clear that both standard geometry 1 and geometry 1 extended back converge to an identical drag value of 0.036 m 2 . Note that geometry 1 extended back converges to a steady state value less efficiently than standard geometry 1, this is expected as having a larger geometry in the simulation involves a larger number of calculations and therefore processing time.

Boundary extension
A final simulation was then completed with the boundary of the simulation extended to surround the entire geometry, so the temperature distribution can be analysed ( Figure 19). The temperature results from this simulation are then compared to the original geometry 1 temperature results, whilst also being compared to an identical geometry with a rounded rear. Note, the effect the rear of the geometry (now included in the full boundary simulation) has on the temperature distribution is a particular focus. Figure 20 shows an identical space shuttle geometry to that of Figure 7, however Figure 20 has a rounded boundary. This boundary extension has the effect of producing a temperature distribution protruding from the rear of the geometry, creating a more realistic temperature distribution than that of the of complete reliance on computational data for aerodynamic optimisation in the future.

Results-boundary extension:
The solver used must first be compared to accurate data from a source outside of computational aerodynamics before further development into a potential 3D aerodynamic package can be produced.
Once further development is completed on this unique solver, a more large scale research project, similar to [18][19] can be completed.