Sustainable space technologies—Strategies toward a predictive aerothermal design of re-useable space transportation systems Open Access

With the advent of new space, the private space industry is seeing a rapid growth of companies looking to capitalize on using space in non-traditional ways. The next decade will see increasing numbers of satellite constellations, space tourism, humans returning to the moon, and the development of building blocks required for future Martian exploration. Contrary to the need for these companies to maximize profit, it is their social responsibility to create sustainable platforms to safeguard the use of space in the future. Re-useable launch and landing systems are able to satisfy both of these requirements. At present, their main system architectures are rockets with vertical takeoff and vertically landing (VTVL) stages to launch payloads and capsules, as well as lifting bodies or gliders for atmospheric entry missions. The relatively new concept of vertically landing rockets was first demonstrated in the Delta Clipper program during the early to mid-1990s¹ and has recently reached technological maturation with the SpaceX Falcon 9.² This family of launchers offers system flexibility, cost efficiency, and reliability on a scale that has never been seen before. Similar booster configurations are also the subject of current research outside of the United States,^3,4 with some studies including advanced concepts that, for example, aim to combine the VTVL concept with airbreathing engines.⁵ Classic designs for the entry of re-useable space transportation systems are based on aerodynamic shapes such as capsules or gliders. The latter offers a substantial increase in cross-range capabilities and reduced peak loads during entry. However, this comes at the cost of higher system complexity and lower volumetric efficiency. The return of launcher stages in glide configuration requires less fuel but is difficult to aerodynamically control due to the unfavorable aft center of gravity location, and to date, there have been no successful demonstrations of these configurations. Future potential single-stage orbit systems like Skylon⁶ are often associated with new airbreathing propulsion concepts and are either limited to small payloads or are still subject to significant technological obstacles in their development. Despite this variety of concepts and missions, the aerothermodynamic design of such vehicles presents similar challenges. Flight trajectories run through continuum and rarefied flow regimes, with chemical and thermal non-equilibrium effects present at high flight velocities. Radiative heat from the surrounding air or from rocket exhaust plumes can also contribute significantly to the total thermal load. Finally, the transition from laminar to turbulent boundary layers substantially increases surface heat flux and modifies the aerodynamic performance of the vehicles. The current state-of-the-art design tools can only approximately account for these complex physical effects. Numerical simulation techniques rely on physico-chemical models, which are often difficult to validate or attempt to describe phenomena such as shock-turbulence interaction or laminar-turbulent transition, which are not yet fully understood, while ground-based testing is limited to small geometric scales and often short test times. As a result, a closely coupled approach of experimental and numerical investigation appears to be the only choice to reduce the significant uncertainties resulting from these design tool predictions. This paper is organized into an overview of current architectures in Sec. II, an introduction to the aerothermodynamic environment and the physical effects encountered during flight in Sec. III, and a review of current experimental and numerical design methodologies in Sec. IV. A summary of different current research foci and an outlook for future design methodologies are given in Secs. V and VI.

One of the main problems with re-useable space transportation systems is the safe return of the vehicle to Earth. High flight velocities during re-entry result in specific total energies of about 30 MJ/kg for orbital velocities and 2.5 MJ/kg for a typical first launcher stage. Since atmospheric drag is the main braking mechanism, this energy is dissipated into the surrounding air, which results in very high temperatures in the shock layer around the vehicle. In addition, retro-propulsion systems rely mostly on atmospheric drag to reduce the required fuel for the entry maneuver. From the vehicle point of view, heating is not the only concern, as the deceleration is limited by the g-loads the payload and vehicle structure can support. Mainly to control thermal loading and maximize atmospheric drag, classical entry vehicles are extremely blunt, which results in poor aerodynamic performance and limited control authority. This and additional constraints on the entry angle in the outer atmosphere limit possible trajectories to a narrow range. Generally, too shallow entry angles result in trajectory overshoot⁷ where the vehicle is not captured by the atmosphere, and too steep angles cause undershoot where the vehicle impacts the ground with a large residual velocity. These boundaries of the re-entry corridor are schematically depicted in Fig. 1.

The second main parameter to control aerothermodynamic loads on blunt vehicles is the entry angle. Generally, a steep entry results in higher thermal peak loads, but, due to the smaller exposure time, less integral load. However, the use of classical rocket-based launch systems imposes constraints on the design of entry configurations. The vehicle's mass is limited by the launcher's payload capabilities, and the vehicle's diameter is limited by the size of the payload bay. This results in a severely limited design space for entry vehicles, and usually, blunt shapes are chosen to reduce thermal peak loads while retaining maximum volumetric efficiency. The design space for capsules is further constrained by the flyable angle of attack, which is limited by the side wall angle to avoid direct impingement of hot gases downstream of the front part. Recent attempts to extend the flexibility of capsule-like shapes are focused around the development of deployable structures based on foldable designs or ballutes, which aim to decrease β and, therefore, the encountered entry loads.^8–11 

The need to increase mission flexibility by achieving higher cross-range and targeting capabilities and more control over the maximum heating rates and g-loads during flight has led to several more aircraft-like configurations such as lifting bodies, where the aerodynamic lift is primarily generated by the shape of the fuselage. Prominent examples of wingless configurations are IXV¹² and HyFlex.¹³ Aerodynamic performance can be further increased by blended wing-body shapes like the X24¹⁴ or X38^15,16 and ultimately by the design of delta-winged gliders, of which the Space Shuttle is the most famous example. The X38, as a typical example of a winged lifting body, is shown in Fig. 2. It was a technology demonstrator for a future crew return vehicle for the International Space Station. The additional control surfaces provided improved control authority and cross range capability. Nevertheless, the aerodynamic performance was not sufficient for soft landing and, hence, a parachute system for capsules was required.

Recent developments of new materials for thermal protection systems with very large lateral heat conductivity have also enabled new aerodynamically efficient designs based on sharp-edged configurations like the SHEFEX vehicle.¹⁷ Another advantage of the faceted shape used in these concepts is the reduced complexity of the thermal protection system due to the applicability of relatively large straight panels.

The new concept of retro-propulsion effectively reduces both peak and integral loads at the cost of additional fuel and engines, which are required for the retro-burn maneuver. Hence, this concept is specifically attractive for launcher stage recovery, where most of the required technologies are already present.^3,18 There are two main mechanisms for load reduction. First, the retro-maneuver reduces the flight velocity and, hence, the specific energy of the vehicle, which needs to be dissipated by aerodynamic drag and, second, the retro-plumes efficiently shield the vehicle from the incoming high-energy free stream.

Retro systems and gliders dramatically move the aerodynamic boundary in Fig. 1. In a limited way, also lifting bodies, which leads to a considerably larger extent of feasible flight corridors in all cases.

A summary of the most important properties of the different systems is given in Table I. The different categories are Q_MAX: maximum heat loading, Q_INT: integral heat loading during entry, η_vol: volumetric efficiency, and g-load: maximum mechanical load during entry flight. Note that all values are estimations for typical atmospheric entry conditions in which lifting bodies and gliders are operated at large angles of attack during the hypersonic flight phase. An aircraft-like mission with a slender body at low angles of attack would result in a significant increase in the integral heat load, Q_INT, which is due to the increased importance of viscous effects.

A typical flow field around capsules is depicted in Fig. 3. The main features are a thin shock layer on the windward side and a system of expansion and re-compression waves on the leeward side. The blunt shape of the capsule generates a very strong bow shock, which is a very efficient means to dissipate the kinetic energy into the thermal energy of the surrounding gas and produce large amounts of wave drag. In fact, only 1%–5% of the total energy is effectively transferred as a heat load to the vehicle surface. The remaining amount dissipates into the atmosphere. Secondary shock structures form in the wake of the capsule and are visualized by contours of negative velocity divergence. While they have no influence on the aerodynamic performance, regions of hot gas form between these shocks, which may cause undesired aft body heating by the radiation of the hot gas. As an example, this effect caused the aft body heating during the Viking mission to be more than two times larger than the initial prediction.¹⁹ Despite the simple and robust layout of capsules, instabilities, especially at blunted cones, occur due to the rapid transition of the sonic line downstream of the bow shock from the nose to the shoulder region. This can be triggered by subtle changes in the gas properties in the shock layer and was initially observed by Gnoffo et al. for the Mars Pathfinder probe.^20,21 Similar effects were also found during post-flight investigations of the Viking Lander.²² A detailed analysis, including ground-based testing, is also given by Hornung et al.^23,24 and highlights the need to accurately predict the thermo-chemical state in the shock layer, even for simple aerodynamic shapes.

Besides the relatively thin high temperature shock layers that are present at blunt capsule-like geometries, there are other physical effects experienced by slender and winged configurations that become progressively more important at large flight Mach numbers. When the high velocity flow is slowed by viscous effects in the boundary layer, the main part of the lost kinetic energy is dissipated into the thermal energy of the gas. Due to the constant pressure in the wall normal direction, this results in a strong decrease in density, which immediately follows from the ideal gas equation of state ρ = p/RT. This results in a rapid growth of the boundary layers, which can exert a major displacement effect on the outer inviscid flow called viscous interaction. For example, the thickness of an adiabatic compressible laminar boundary layer on a flat plate scales with the square of the free stream Mach number.²⁵ Another important effect is the presence and development of entropy layers, which are a result of the highly curved shock wave at the blunt leading edges of vehicle geometries, such as wings. The entropy of the flow increases across a shock wave, and a streamline passing through a nearly normal portion of the shock will see a larger entropy rise than a neighboring streamline passing through a more oblique portion. Hence, a layer of strong entropy gradients forms downstream of the blunt leading edges. An illustration of the canonical flow fields over a wedge and a cone for a free stream Mach number of 8 is given in Fig. 4. It can be seen that the entropy layer thickness (depicted by the red line) decreases for conical configurations due to the convergence of the stream surfaces close to the body (effect of mass conservation), whereas it stays constant at 2D wedges.^26,27 Close to the leading edge, boundary layers (depicted by the blue contour) are generally thin compared to the entropy layer. The influence of the entropy layer vanishes when it is merged into the boundary layer. As visible in Fig. 4, this happens much earlier in conical flows (e.g., rockets) than in 2D wedge flows (e.g., wings). In the present example, both layers reach the same thickness in the conical flow at an axial distance of about x/R_NOSE = 7 and much further downstream in the wedge flow. Analytical estimations for this behavior are given, e.g., by Stetson.²⁸ 

To further illustrate the consequences of the presence of the entropy layer, Fig. 5 shows the wall normal profiles of flow velocity, total enthalpy, and entropy of the wedge configuration in Fig. 4 at x/R_NOSE = 10. The region of the classical boundary layer that is dominated by viscous effects is visible in the profile of total enthalpy H/H_∞. Between the edges of this boundary layer and the entropy layer exists a large region with boundary-layer-like velocity gradients that are not related to viscous effects. Hence, large portions of the velocity profile are dominated by the entropy layer effect. This affects many classical boundary layer analyses because the conditions at and the location of the outer boundary layer edge are not clearly defined anymore. Furthermore, entropy layers, which are generated by swept leading edges, can cause strong adverse cross-flow effects.²⁹ Now, flow properties close to the vehicle surface are not only governed by viscous effects but are strongly affected by the design of leading edges, which can significantly affect aerodynamic performance through earlier boundary layer separation at adverse pressure gradients, altered laminar to turbulent transitions, and the initiation of undesired cross-flow.

Another critical phenomenon is the shock–shock interaction. If a shock wave (e.g., generated by the fuselage) impinges on the bow shock of a wing, as illustrated in Fig. 6, the heat and pressure loads can be severely amplified. There are six basic gasdynamic interaction patterns, the so-called Edney types.³⁰ The most dangerous configuration is the type IV interaction, where the impinging outer shock hits the local bow shock close to the stagnation streamline, where the undisturbed downstream flow of the bow shock would be subsonic. This interaction is characterized by the formation of a strong and confined supersonic jet downstream of the interaction point, which causes extremely high surface loads. Local amplification of heat fluxes by factors of 5 and up to 19 was observed in wind tunnel experiments at Mach 6 to 20 conditions.^31,32 The first in-flight confirmation of the severity of shock interference heating took place in October 1967 when the NASA X15A2 experimental rocket-driven aircraft suffered severe damage to one of its fins,^30,33,34 as illustrated in Fig. 7. However, most of the current re-entry configurations use highly swept wings, where shock interference is dominated by Type VI interactions³⁴ and the resulting amplification of surface heat loads is less severe. As an example, a factor of 2.7 was observed at the Space Shuttle leading edge in the vicinity of fuselage-wing shock interference³⁵ in a similar situation as illustrated in Fig. 3. Flow patterns caused by shock interference occur on small geometrical scales and are difficult to accurately predict with design tools.

Another gasdynamic phenomenon is the interaction between shock waves and boundary layers. Since both can be found in every supersonic flow, this type of interaction is commonplace and relevant for practical designs. It occurs when an externally generated shock wave, or one that is generated by strong curvature, jumps on the body surface and interacts with a boundary layer. This may have major consequences for surface loads and entire flow fields.³⁶ The influence of the shock imposes a strong local adverse pressure gradient, which leads to the thickening of the boundary layer or even the development of local flow separation zones. For separated cases, an additional re-attachment shock is formed downstream of the interaction zone. Such flow patterns strongly promote the transition from laminar to turbulent boundary layers and cause strong heat flux amplification at the re-attachment point.³⁶ The entire flow pattern is also a strong function of the local wall cooling and the effect of overlayed entropy layers.

Besides capsules and gliders, the alternative option for re-useable entry technology is retro-propulsion. Here, a wide range of additional gasdynamic effects occur. These are related to the interaction of the propulsive jets with the supersonic free stream. As an illustration, the flow structure around a generic first-stage retro-burn maneuver at high altitude³⁷ is depicted in Fig. 8.

The flow field in this figure shows that a detached bow shock forms upstream of the exhaust plume of the three operational engines. The large supersonic exhaust plume efficiently shields the vehicle from the direct impact of the high energy incoming free stream. After passing the bow shock, the free stream air is deflected by the jet plume, and a contact surface between exhaust and air is created in the stagnation zone (visible in the figure through the temperature gradient between the Mach disk and bow shock and the brown exhaust boundary). The plume itself forms a characteristic barrel shock and a Mach disk. Due to the large expansion of the plumes at low atmospheric backpressure, they merge and interact in the vicinity of the launcher base, which results in a pattern of plume interaction shocks. This interaction generates local stagnation zones between the plumes and can give rise to a significant backflow of hot gases to the base plate. It is seen that the entire vehicle is immersed in its exhaust gases. The stagnation zone is characterized by high temperatures. The rocket exhaust is stagnated and, hence, the combustion chamber temperature is almost recovered. Additionally, the incoming air flow is heated by the strong bow shock. Despite the very large temperatures of the gas surrounding the entire vehicle, in the order of 2000 K, the vehicle experiences moderate to low heat loads due to the low density. In fact, the most critical point concerning both thermal and mechanical loads is the shut-down of the retro-engines when the vehicle loses the shielding effect of the retro plumes and is directly exposed to the still-high supersonic free stream.³ It has to be noted that the current knowledge of the aerothermodynamics of slender rocket-like retro-propulsion configurations is still very limited.^3,38–40 Previous studies focused on blunt body technologies,^41,42 which are mainly relevant for entry into low-density atmospheres, and the development of rockets is strongly driven by private companies, with the resulting limited sharing of flight and design data.

A detailed discussion of the physics of turbulence is not the focus of this paper. Nevertheless, it is an important problem for vehicle design as the heat transfer rate to cooled surfaces increases by a factor of 3–5 compared to laminar flow. Furthermore, the friction drag significantly increases for slender vehicles, and boundary layer separation characteristics change, which may strongly influence the performance characteristics and operational margins of control surfaces and air intakes. The state of the boundary layer (laminar or turbulent) has important consequences for vehicle design. Besides the impact on friction drag and surface heat loads, it also defines allowable surface roughness, steps, and gaps to maintain an aerodynamically smooth surface. Laminar-to-turbulent transition mechanisms are very complex⁴³ and not yet fully understood. The main mechanisms are a natural transition with the growth of streamwise unstable two-dimensional waves (Tollmien–Schlichting waves) until an amplitude where non-linear effects cause breakdown to turbulence. Furthermore, there are cross-flow instabilities that can occur on swept surfaces. They are primarily driven by the inviscid effects of surface shape and external pressure gradients, which produce curved streamlines. Inside the boundary layer, the velocity is reduced by viscous effects, but the pressure gradient being imposed by the inviscid outer flow remains unchanged. Hence, there is a generation of secondary flow perpendicular to the inviscid streamline. Because the cross-flow velocity is zero at the wall and the boundary layer edge, an inflection point exists, which gives rise to an additional instability mode.⁴⁴ An additional principal mechanism is the second mode or Mack mode transition.⁴⁵ It occurs for slender bodies at flight Mach numbers above 6 due to trapped acoustic waves inside the boundary layer. Other transition mechanisms include the bypass transition, which involves the direct initiation of unstable disturbances by large free stream turbulence, which mainly occurs in wind tunnels, and the effect of Görtler vortices, which tend to generate unstable vortex structures at concave surfaces. Generally, the prediction of laminar to turbulent transitions for hypersonic flows is extremely difficult and is the subject of current research. Satisfactory engineering correlations do not exist except for canonical problems like flat plate boundary layers. Indeed, popular correlations to Reynolds numbers, momentum thickness Reynolds numbers, and boundary layer edge Mach numbers vary by two orders of magnitude.⁴⁶ Presently, in vehicle design, this means that one has to consider both laminar and turbulent states over a wide range of flow conditions and assume the worst-case scenario for vehicle design, which often leads to excessive safety margins.

Generally, during hypersonic flight and atmospheric entry, the kinetic energy of the atmospheric gas (in the body fixed reference system) is dissipated into thermal energy when it passes the bow shock of a vehicle. This results in a dramatic and sudden increase in thermodynamic energy. Subsequently, collisional interactions between particles give rise to a cascade of several chemical and thermodynamic kinetic processes. A schematic representation of key processes occurring in the shock layer of a typical re-entry vehicle along the stagnation streamline is shown in Fig. 9. For simplicity, the additional effects of ablating shock layers are left out.

First, the kinetic energy of the gas particles is rapidly converted into the translational energy of the molecules when they collide with the dense shock layer gas. Inter-particle collisions then excite first the rotational and then the vibrational and electronic modes of the molecules. Rotational and translational energy modes reach an equilibrium state rapidly, in the order of tens of collisions. Vibrational–vibrational excitation is considerably slower and takes thousands of collisions.⁴⁷ During the excitation process, molecules collect vibrational energy to the point where the inter-molecular bonds are overcome and dissociation occurs. Molecules in higher vibrational states have a higher probability to dissociate due to their lower dissociation potential—this behavior is called preferential dissociation. Generally, oxygen dissociation occurs first at a temperature of 2000–2500 K, whereas nitrogen dissociation effects are negligible below temperatures of 4000 K. Already at lower temperatures, nitric oxide is formed by exchange reactions (e.g., N₂ + O = NO + N). Furthermore, collisional processes between molecules and atoms excite the bound electrons to elevated states, and the gas radiates electromagnetic energy as the electrons spontaneously decay to less energetic states. The coupling between these excitation and de-excitation processes occurring at different rates can lead to non-equilibrium electronic population distributions.⁴⁸ Given sufficiently energetic collisions, bound electrons can eventually be stripped from their parent nuclei. In air, electrons are first produced by the associative ionization of the NO molecule. The resulting NO⁺ ion then transfers its charge to other neutral species through several charge-exchange reactions.⁴⁷ When the resulting density of the electrons reaches a certain threshold value, a cascade of electron production can occur through the highly efficient electron-induced ionization of N and O atoms.

Sufficiently far downstream of the bow shock, the plasma reaches a state of local thermodynamic and chemical equilibrium. The relaxation length to equilibrium strongly depends on the shock layer density and temperature. As the gas enters the boundary layer, heat convects into the vehicle surface, and the gas temperature drops. This results in the recombination of atoms and ions. The potential catalytic effects of the vehicle surface promote this process. The presence of ablative thermal protection systems and the associated injection of ablation products can further significantly influence thermo-chemical processes in the boundary layer.

A velocity–altitude map that relates the occurrence of the thermo-chemical effects from Fig. 9 to the four characteristic flight trajectories of re-useable space transportation systems is shown in Fig. 11. Due to discrepancies between the boundaries of the phenomena and unknown underlying assumptions in previous publications,^49–52 the plot is redone here based on simple shock layer analysis. First, the atmospheric conditions for a given flight altitude are evaluated based on the US 1976 standard atmosphere. Then, the conditions behind a normal shock wave are computed for a given flight velocity and the atmospheric conditions at the given altitude. These post-shock conditions are calculated with the assumption of thermal equilibrium between the translational, rotational, vibrational, and electronic excitation of the molecules and for a frozen chemical composition. In this post-shock state, a chemical non-equilibrium analysis in a constant-pressure control volume is performed using an 11 species—23 reactions chemistry rate mechanism for Earth entry.^53,54 As the pressure downstream of the bow shock is approximately constant along the stagnation streamline, the model of a constant pressure reactor is a useful approximation to analyze chemical relaxation to equilibrium, which greatly simplifies the calculation procedure compared to full computational fluid dynamics (CFD) solutions. From these results, different characteristic boundaries of aerothermodynamic phenomena are extracted and briefly discussed in the following paragraphs.

First, the boundaries for vibrational excitation (yellow dotted lines labeled with e-vib) are directly computed for the post-shock state. The figure includes the thresholds where the energy of the vibrational modes is 0.1 RT and 0.5 RT, where R is the specific gas constant. This vibrational energy refers to the mass-fraction weighted sum of the individual contributions of the N₂ and O₂ molecules. Note that the maximum vibrational energy of a diatomic molecule is RT.

Next, to assess the significance of chemical reactions, the molar fraction of atoms in the equilibrium state at infinite time was taken. The plot contains the limits of 5% and 20% of the atomic molar fraction in the equilibrium shock layer gas (which corresponds to almost complete dissociation of oxygen) depicted by the dotted red lines labeled with X-atoms.

Then, the ionization limits are depicted by the blue dotted lines, which refer to a mole fraction of 1% and 5% of ionized species in the equilibrium state.

The virtual equilibrium value of the vibrational energy, $evibeq$⁠, is computed from the local temperature of the reactor analysis, which was performed under thermal equilibrium assumptions. The initial condition for e_vib is its free stream value. The vibrational relaxation length, r-vib, is then determined from the location where the vibrational temperature deviates less than 5% from its local equilibrium value. The relaxation time, τ, is estimated from the Millikan–White approximation⁵⁵ with Park’s high temperature correction.⁵⁶ The goal of the analysis was to obtain an estimate for the total length scale needed for thermal relaxation; hence, it was performed for the slowest relaxing nitrogen molecule.

To validate the present analysis procedure, cross-check CFD computations of the shock layer around spherical bodies were performed with the DLR TAU code.⁵⁷ The radius of the body was chosen such that the resulting shock layer thickness allows relaxation to reach equilibrium. Example results for a flight velocity of 4000 m/s at an altitude of 60 km are shown in Fig. 10. The solid lines are the CFD results for the translational/rotational temperature of the shock layer gas and the vibrational temperature of N₂. The vertical dashed lines correspond to the estimated relaxation lengths, r-chem and r-vib, from the above analysis. This flight condition is characterized by strong chemical non-equilibrium, and the chemical relaxation takes about 1 m distance downstream of the shock. The moderate vibrational nonequilibrium needs a relaxation distance of ∼0.2 m. The agreement between the reactor analysis and the CFD is strong in this case. A maximum error of about 40% was observed in a large number of similar cross-checks, which is acceptable for the approximate purpose of Fig. 11 since both chemical and thermal relaxation length scales vary by several orders of magnitude.

The boundaries for the ratio of the radiative to the convective heat flux in Fig. 11 are computed using the shock layer properties obtained in the reactor analysis. A constant shock layer thickness of 0.1 m was assumed, and the radiation database PARADE⁵⁸ was used to obtain the total radiative flux in a wavelength range of 100–6000 nm and a spectral resolution of 10 000 points at equal frequency increments. Self-absorption at large densities was accounted for by using an infinite slab model⁵⁹ to compute the total transfer to the surface. The convective transfer was estimated with the Verant correlation,⁶⁰ which is based on numerous CFD analyses of the stagnation point heat flux in spherical configurations over a very comprehensive range of flow conditions. The resulting laminar convective heat transfer rate was gradually amplified up to a factor of three for Reynolds numbers above 10⁶ to account for the onset of turbulent flow.

To indicate the beginning of turbulent flow, the flight Reynolds number based on a reference length of 1 m, Re, is indicated by the brown-shaded area. Below a value of Re = 10⁶, the flow is likely to be laminar. For large configurations like rocket stages, the transition region to turbulent flow will shift slightly toward higher altitudes due to the increased length scale of the vehicle.

Finally, four typical reference trajectories are included in the figure. They range from super-orbital lunar return missions, over high lift (glider), and ballistic (capsule) entry from low earth orbit to the typical return trajectory of a VTVL stage. All trajectory data are taken from the very comprehensive book on space systems by Suresh and Sivan,⁶¹ except for the first stage entry, which corresponds to a generic configuration that was designed and analyzed in the RETALT project of the European Commission Horizon 2020 program.^3,18

From Fig. 11, the main aerothermodynamic phenomena that occur during the different mission types and phases can be identified.

At altitudes above 85 km, the heating is dominated by shock layer radiation. Vibrational and chemical non-equilibrium effects (yellow and red dashed lines) are negligible for most of the trajectory. They become important only for flight velocities between 2200 and 3000 ms/s (molecular vibration) and between 3800 and 2800 m/s (non-equilibrium chemistry). Down to an altitude of about 50 km, the vehicle is surrounded by strongly ionized flow. The most critical parts (concerning heating) of the re-entry occur in the equilibrium regime. This is due to the high temperatures (caused by the large flight velocities) and the high density (caused by the relatively low altitude of the high velocity part) in the shock layer, which both promote rapid equilibration of chemical kinetics and thermal excitation. The boundary layer transition to turbulent flow is likely to occur at altitudes below 30 km and at moderate velocities.

For both LEO entries, ionization is only significant at altitudes above 90 km. Above 80 km, the radiative heat flux dominates over the convective contribution. However, due to the low densities, the absolute level remains low at these flight conditions, and the contribution to the integral heat load is usually negligible. Contrary to the hyperbolic entry, the main part of the flight path is within the chemical non-equilibrium regime. Thermal non-equilibrium occurs at altitudes of around 50 km. Transition to turbulent flow occurs at comparatively low flight velocities where the total heating rate is not critical due to the reduced total energy of the flow.

Ionization, radiation, and chemical dissociation are negligible during the entire trajectory. Note that this analysis is only valid for the aerodynamic flight phases without propulsion. During a retro-burn, complex thermochemical interactions occur due to the interaction of the atmospheric gas with the exhaust plumes. In addition, radiative heat loading might occur due to the presence of soot in the exhaust plumes of kerosene-fueled rockets. These effects are excluded in the present analysis. The map in Fig. 11 also indicates the presence of strong thermal non-equilibrium effects. The actual importance of those is uncertain in this regime. It is known that the applied modeling approach of vibrational relaxation is not well validated in this region, which specifically concerns the applied correlations for the vibrational relaxation times at the present low temperatures for the moderate flight velocities.⁶² For practical application, it appears to be justified to neglect vibrational non-equilibrium effects in the low temperature limit below around 2500 K.

The results in Fig. 11 represent only a rough global estimate of the likelihood of the occurrence of aerothermal effects. This especially applies to the estimates of the non-equilibrium boundaries for chemical and thermal relaxation. The underlying analysis is for blunt body shock layers. Slender bodies such as gliding vehicles tend to be subjected to enhanced non-equilibrium effects. This is because the associated oblique shocks result in decreased post-shock temperatures, which increase both the chemical and thermal relaxation times. The Reynolds-number based indication for the onset of turbulent flow relates to small vehicles with a characteristic length scale of 1 m, and a slight shift toward higher altitudes occurs for large-scale configurations. Hence, for accurate statements, specific and detailed analyses are needed. This also applies for the high velocity and high altitude ranges where an accurate analysis of shock layer radiation is required. A very small region downstream of the shock can contribute significantly to the total heat load due to the strongly non-linear dependence of emitted radiation on the temperature. Even if this is associated with a small relaxation length in the present map, non-equilibrium might still be important to quantify. This is also true for small regions of ionization, which cause the blackout effect during entry flight. Nevertheless, the results still give a general overview of the expected phenomena. Note that the lower right part of the figure with the associated strong influence of radiative heat flux is not relevant for practical applications, as the total heat loads in this regime are far too large for any practical vehicle design.

It has already been established that during flight, spacecraft are subjected to a wide spectrum of freestream conditions, ranging from low Mach number incompressible flows to high Mach number compressible flows with thermo-chemical effects. No single ground-based test facility is able to account for these environments through the entire flight envelope, and as a result, it is critical that engineers have access to a variety of specialist wind tunnels, each of which has been designed to simulate a specific phase of flight. This fact is best highlighted by the Apollo wind tunnel testing program,⁶³ where 25 facilities were used to cover the various vehicle configurations through a Mach number range of 0 to 20. Since the 1960s, there have been significant advances in numerical techniques and computing power. This means that such a large test campaign will likely never be repeated. However, it does not detract from the importance of conducting such experiments, particularly for validation purposes.

Wind tunnel selection is largely governed by the facility's ability to simulate the desired flight conditions. This can be quantitatively evaluated using similarity parameters, which are values used to describe the flow state. For low enthalpy flows (less than 2 MJ/kg), perfect gas assumptions are valid, with only Mach and Reynolds number similarity required.⁶⁴ The Mach number represents the ratio of the flow velocity to the speed of sound in the gas, and it is important for replicating compressible flow effects as well as shocks. The Reynolds number describes the ratio of inertial to viscous forces for a given flow, where low values indicate the fluid behavior is likely laminar and high values suggest turbulence. Attaining Reynolds number similarity between flight and experiment means that the boundary layer is replicated, which is particularly important in understanding phenomena like laminar to turbulent transition. For low enthalpy conditions, both the Mach and Reynolds numbers seen in flight should be replicated experimentally in a wind tunnel; however, this is often not practically possible. Strong limitations are imposed on model size due to the available dimensions of the test section and mounting infrastructure, which in turn limits the Reynolds number of the flow due to its dependence on model characteristic length. Over the years, aerodynamicists have employed strategies to overcome these limitations and alter the Reynolds number of the flow to be more indicative of flight, which include boundary layer trips,⁶⁵ changing the working fluid of the facility,⁶⁶ or creating pressurized wind tunnels, where the influence of Reynolds number can be isolated for a given Mach number due to density variations.⁶⁷ 

As the flow enthalpy is increased above 2 MJ/kg, the freestream Mach number becomes progressively irrelevant due to hypersonic similarity.²⁵ On the other hand, high-enthalpy effects emerge, specifically excitation of vibration, dissociation, and ionization modes as the air travels through strong shocks. This means that the flow velocity becomes a critical factor in place of the Mach number, as well as the Damköhler number, which relates the vehicle body length and the dissociation relaxation distance, which is approximately equal to the inverse of mass or molar density. This means that when correctly replicating the distribution of species present due to chemical dissociation reactions initiated by a strong shock, ρL similarity must be satisfied. In other words, increasing the freestream density can offset the impact of having a sub-scale model when completing tests looking to replicate aerothermochemistry effects. It is worth noting that, provided the flight velocity and gas viscosity are matched in the ground-based facility and ρL scaling is adhered to, the flight Reynolds number is also preserved, making these two similarity parameters tightly coupled. While ρL scaling duplicates the properties of a dissociation dominated flow, it does not preserve chemical equilibrium properties in the shock layer or regions with dominant ternary recombination reactions (e.g., the boundary layer). The reproduction of the thermodynamic state of the test gas is required to match the equilibrium conditions, whereas ternary recombination dominated flows scale with ρ²L. The consideration of similarity parameters can make facility selection difficult, and full-scale testing is the only means to reproduce all of them. In cases where matching all key criteria is not possible, engineering judgments need to be made to acknowledge and understand the limitations of the collated data. In some cases where specific component subsystems are of interest, wind tunnel experiments can be designed such that the correct conditions upstream of an isolated vehicle part are created such that the entire model is not needed. This reduces, or in some cases eliminates, the scaling factor required by the test article for the experimental campaign, allowing, e.g., higher Reynolds numbers to be realized.

Figure 12 gives an overview of different ground based test facilities and their approximate operating ranges in terms of altitude and velocity. Once again, this figure highlights the need for various test facilities to cover the entire range of flight conditions for vehicles returning from LEO. For reference, recall from Fig. 11 that the Space Shuttle reached speeds of around 7.8 km/s, while a first-stage entry is limited to around 2.2 km/s at an altitude of 100 km.

In terms of similarity parameters for low enthalpy flows, Mach number similarity can be easily achieved in blowdown wind tunnels. Upstream of the test section is a plenum, which is fed with test gas from a high-pressure vessel. The flow is then expanded through a nozzle, where the desired test Mach number is determined by the area ratio between the nozzle exit plane and the throat. At the exit of the test section, the mass flow is regulated through the use of a diffuser before flowing into a dump tank. As a blowdown tunnel is not a continuous flow facility, the test times are governed by the size of the pressure reservoir relative to the nozzle throat area, with typical test times in the order of seconds to minutes. Due to the high pressures and long test times in these tunnels, high Reynolds numbers can be achieved, which makes these facilities suitable for performing comprehensive aerodynamic analyses. As the requirement for Mach number increases, so does the pressure ratio between the high-pressure and low-pressure sides of the tunnel.⁶⁸ For crudely designed blowdown tunnels that eject flow to sea level ambient conditions, the high-pressure required to achieve the desired conditions in the test section becomes too large. This can be alleviated by attaching a vacuum storage tank to the exhaust side, allowing higher pressure ratios to be realized. A basic diagram of this tunnel is provided in Fig. 13.

It is important to consider that, due to the principle of energy conservation, increasing the flow of kinetic energy must be offset by a reduction in thermal energy. Therefore, the acceleration of the air through an expansion in a nozzle will result in a decrease in temperature at the exit. Where air is used as a test fluid and large expansion ratio nozzles are used for high Mach number operation, the air must first be passed through a dryer bed and heated to prevent moisture as well as oxygen and nitrogen condensation. For this reason, blowdown facilities are considered “cold,” as the temperature of the air in the test section is often just above the condensation temperature to avoid liquefaction during the expansion process. This results in a reduced speed of sound for the test gas. Hence, in these tunnels, Mach number similarity can be achieved with comparably low flow velocities, resulting in lower enthalpy and stagnation temperatures when compared with flight conditions. This is where the terms hypersonic and hypervelocity need to be clarified. While the definition of hypersonic is somewhat controversial, it is generally considered to be around Mach 5 and above. Hypersonic Mach numbers can be easily achieved in blowdown tunnels due to the temperature loss during the expansion process through the wind tunnel nozzle. Hypervelocity is, therefore, the correct term to distinguish the artificially high Mach number flows with low freestream velocities from the high Mach number, high velocity flows experienced in flight. This leads to one of the main disadvantages of blow-down tunnels, which is their inability to produce shock-layer temperatures consistent with vehicles traveling in the hypervelocity regime.

To remedy this deficiency, there are traditionally two methods used by experimental facilities to achieve high flow enthalpies through the heating of the test gas. The first makes use of a conventional blow-down tunnel architecture with the addition of electrical arc heating and is aptly named the arc heated tunnel. High power is required to run the heater, with some facilities using up to 70 MW⁶⁹ and being able to realize enthalpies up to 45 MJ/kg. One major advantage of these facilities is the long test times, which are afforded by the large blow-down air storage systems. On the other hand, the operating pressure and mass flow are limited due to the heating process.⁷⁰ As a result, the freestream densities are low, which restricts the facility's capabilities in terms of Reynolds number and ρL scaling.⁷¹ In addition, arc heaters use copper electrodes, which experience some electric erosion. This causes flow contamination through the presence of copper particles and flow non-uniformity. For this reason, these facilities are often used for materials testing, particularly thermal protection systems and ablators, rather than extracting accurate load cell measurements for aerodynamic coefficient derivations. An overview of an arc heated facility is given in Fig. 14.

The second method of achieving high enthalpy flows is through the exploitation of shock and expansion waves in shock and expansion tunnels (schematics in Figs. 15–17). Shock tunnels consist of three main parts, the driver section (driver or compression tube), the driven section (shock tube), and the test section. The high pressure driver section is separated from the shock tube by a diaphragm, which is designed to rupture once a certain pressure differential has been reached. At this time, an incident shock propagates down toward the test section and is reflected at the end of the shock tube. During this time, the driver and the test gas are separated by a contact surface. After the incident shock wave is reflected at the end of the shock tube, the test gas is brought to rest, and the shockwave penetrates the contact surface. This leaves a zone of high temperature, stationary test gas at the end of the shock tube, causing the secondary diaphragm to burst, allowing the test gas to be expanded into the test section of the tunnel (see Fig. 15).

Initial designs of shock tunnels selected room temperature helium or hydrogen as the driver gas, facilitating flows with enthalpies of around 5 MJ/kg in the test section. Heating the driver gas enables one to increase this value to ∼12.5 MJ/kg. However, at these enthalpies, nitrogen dissociation does not occur, and as a result, a significant portion of the re-entry trajectory of spaceplanes could not be assessed experimentally.⁷² This led to a design evolution of the shock tunnel to include a free-moving piston in the compression tube (see Fig. 16). High pressure air is stored in a reservoir behind the piston, which initiates the piston acceleration down the compression tube, rupturing the diaphragm between the compression tube and shock tube and initiating the process explained previously for a conventional shock tunnel. The major advantage of this approach is the rapid compression and heating of the driver gas due to the motion of the piston, allowing a significant increase in the flow enthalpy in the test section (above 20 MJ/kg⁵²). Referring back to the previously mentioned ρL scaling factor, higher expansion ratios of the nozzle result in lower flow densities in the test section and, therefore, can influence this similarity parameter. Considering this, it needs to be noted that as long as the flow is greater than approximately Mach 8, the vehicle aerodynamics are only weakly coupled with the Mach number.⁷² In practice, this means that lower Mach numbers can be used to aid the ρL scaling. Typical test times depend on tunnel configuration and range from 0.5⁷⁴ to 150 ms.⁷⁵ Shock tunnels can be modified to be operated in non-reflecting mode, which results in the direct expansion of the central core of the test gas to hypersonic conditions. This facilitates the reproduction of the exact shock conditions experienced in re-entry at true local pressure at the expense of short test times and a poor ability to simulate aerodynamic flows. These facilities are primarily used for studying the immediate post-shock relaxation and radiation processes.

Finally, expansion tube facilities (see Fig. 17), are able to attain higher freestream velocities than shock tunnels through the addition of an acceleration tube. The test is initiated when the primary diaphragm ruptures from a pressure increase in the driver gas. This causes the formation of a primary shock, which propagates through the shock tube section of the tunnel. The secondary diaphragm bursts once it is reached by the primary shock. This results in the simultaneous formation of a secondary shock wave, which travels down the constant area expansion section, as well as an unsteady expansion, which moves back into the test gas. This means that the test gas undergoes two processes: first, compression and heating by the primary shock, before being accelerated and cooled by the unsteady expansion. The basic distinction from other impulse facilities is that the unsteady expansion process adds energy and total pressure to the test gas through an unsteady cascade from the upstream to the downstream slugs of gas. The final thermodynamic state of the test gas is a function of the individual section pressures and gas compositions, giving flexibility to the possible test conditions. In addition, the chemistry of the flow can be better regulated because the test gas is not subjected to high intermediate temperatures. When compared to the standard reflected shock tunnel, higher velocity flows and enthalpies can be realized, in some cases in excess of 100 MJ/kg,⁷⁶ which corresponds to a free stream velocity of more than 14 km/s. As shown in Fig. 17, there is no nozzle present, and the high speed gas is directly discharged into the test section. This results in test times in the order of 50⁷⁷ to 8000 µs.⁷⁸ Often, a diverging nozzle is attached to the end of the expansion tube to increase the size of the core flow region.

The descriptions of the test facilities given above highlight that the main determining factor for instrumentation and data acquisition requirements will be the test duration. This inherently influences the initial decision of facility selection, as while a facility may be able to replicate the desired freestream conditions, the required diagnostic tools may not be deployable or available. For low enthalpy, long-test-duration facilities such as blowdown tunnels, the measurement techniques used for subsonic and supersonic flows are still valid and applicable. Where stagnation temperatures do not exceed 1000 K, even hot wire anemometry for freestream turbulence characterization is possible.⁷⁹ As the facility test times decrease into the millisecond and microsecond ranges, some challenges begin to arise. Traditional load cell-based force measurement techniques do not have sufficient response times, which makes the determination of aerodynamic loads in impulse facilities extremely difficult. Recently, there have been advances in methods including stress-wave balances⁸⁰ and instrumented free-flying models accompanied by visual tracking techniques,⁸¹ which allow the execution of force measurements. Typically, diagnostic tools can be divided into three categories, intrusive, semi-intrusive, and non-intrusive. Intrusive techniques include the addition of measurement equipment to the surface of the model and include thin film surface thermometers, thermopiles, and coaxial thermocouples.⁸² Semi-intrusive methods apply to uniform surface coverings like temperature sensitive paint (TSP) or pressure sensitive paint (PSP), while non-intrusive techniques generally include purely optical methods such as infrared thermography or spectroscopy. While each of the sensor based tools has its own limitations in terms of response time, durability, temperature sensing ranges, and cost,⁸³ additional considerations need to be made for surface coatings. Traditionally, TSP and PSP have suffered from long response times relative to the test duration of impulse tunnels, and it is only in the past 5–10 years that fast response TSP and PSP have been developed.^84,85 One significant advantage of surface coatings is their surface resolution, which is determined by the pixel size in the photographs and represents a departure from the small number of sensors that provide a discrete dataset. However, there are some strong limitations to these surface coatings when they are considered for use under high enthalpy conditions. Both TSP and PSP can encounter issues once test gas self-luminosity sets in and emission bands similar to those of the surface coatings are emitted, making imaging a non-trivial task.^86,87 Furthermore, issues pertaining to the robustness of the coating, including surface abrasion, cracking, or charring, can cause issues in result processing. Finally, because PSP is fundamentally an oxygen concentration measurement, at high flow enthalpies where the concentrations of O₂ are not known or are changing, quantitative PSP is not possible.⁸⁷ In long-term test facilities, heat flux analyses can be easily conducted, but one drawback is the influence of lateral heat conduction on the wall normal heat transfer. The impact of this can be reduced by designing thick walled models with low thermal conductivity.⁸⁸ The lateral conduction problem is alleviated when testing in short duration facilities, particularly shock tunnels, as the test time is long enough for a steady state flow to be established but short enough to limit lateral heat conduction.⁷⁹ In any case, knowledge of the flow field is required to interpret the local sensor results since they can be strongly affected by features such as shock impingement or flow separation. For this reason, Schlieren photography also plays an important role in correlating the measured values with a flowfield, particularly for the determination of shock location, standoff distance, and any interactions between flow features. Spectroscopy can be used to identify the species and their concentrations in the flowfield. Either the emission of the excited test gas is measured or the absorption of probe light by the test gas. While emission spectroscopy is a relatively simple technique, its main limitation is that the gas must be luminous, which is typically only seen in high enthalpy facilities. The results generated from emission spectroscopy are intensity as a function of wavelength. The identification of species present in the gas can be deduced from these data because different atoms and molecules have different wavelengths associated with their light emissions. Using fitting techniques based on known emission spectra, these data can be used to identify the species present, their concentrations, as well as their electronic, vibrational, and rotational states. On the other hand, the absorption technique relies on the amount of light that is absorbed by the test gas. This requires the experimental setup to have a light source on one side of the tunnel and a detector on the other. This tool has the ability to determine pressure, temperature (translational, rotational, and vibrational), velocity, species concentrations, and density. From these data, mass flow rates can also be found. Typically, this method is used to measure these data for a single species; however, investigations using more sophisticated laser and optical setups allow one to measure the properties of multiple species.⁸⁹ 

The selection of a ground based facility and diagnostic tools is guided by the type of vehicle to be analyzed and the trajectory points that can be covered. In recent years, the concepts behind re-useable space transportation systems have changed considerably when compared with the Space Shuttle. This has been largely driven by private companies such as SpaceX and their development of Falcon 9 and Starship. The first stages of these vehicles follow sub-orbital trajectories after stage separation and perform a series of retro-propulsion burns, which as mentioned in Sec. III, are responsible for complex thermochemical interactions between freestream air and the hot exhaust plumes. This shifts the main challenges for ground-based facilities from the high velocity and enthalpy constraints for re-useable vehicles returning from LEO to the influence of exhaust plumes that surround the vehicle as it performs burns. Factors such as post-combustion due to the presence of atmospheric oxygen and flow stagnation from the hypersonic Mach number counterflow result in the vehicle being immersed in high temperature reacting flows. At high altitudes where the nozzles are operating in an under-expanded state, additional factors such as post-expansion and plume–plume interactions begin to influence the aerothermal environment due to the impingement of the flow on the baseplate, creating regions of locally high heat flux. Finally, hydrocarbon fuels like kerosene eject soot particles, which contribute to radiative heat fluxes and can make up a large percentage of the total heat loads.⁹⁰ This means that to simulate the aerothermal environment, wind tunnel models need to be capable of combusting fuel and oxidizer, ideally at ratios representative of full-scale propulsion systems. This presents significant challenges to the model designer, as the slender cylindrical shapes of sub-scale booster models now not only need to have space for instrumentation and the associated cabling but also thermal management systems due to the high combustion temperatures and internal ducting from a combustion chamber to the nozzle exits. These challenges are not only limited to boosters but apply to all models that are utilized for retro-propulsion tests. Figure 18 highlights the complexity of constructing a sub-scale model with an internal flowpath. Model support systems need to be modified to allow the flow of exhaust gases to the model while a plenum feeds the eight separate jets on the underside of the body. Due to the scaling requirements of wind tunnel models, nozzle exit diameters are typically measured in millimeters, making the manufacture and integration of these flowpaths extremely difficult.⁹¹ 

Thermal stresses on the internal flowpath can be reduced if combustion chamber temperatures are lowered through a reduction in the oxidizer-to-fuel ratio. Unfortunately, this is not a perfect solution, as the excess fuel has been shown to react with atmospheric oxygen, causing significant post-combustion emissions that and not only shift the requirement of thermal management from the interior of the model to the exterior but also result in excessive, non-flight representative heat loads on the vehicle surface.⁹² The post-combustion is visualized in Fig. 19, where subsonic aerothermal tests were conducted during the RETPRO project³⁸ in the vertical test section (VMK) at DLR Cologne.

To create a simplified test campaign, the option to replace a chemically reacting plume with air is also a viable option. The physics of having pressurized air create an exhaust jet representative of a plume is similar to that of a blow-down wind tunnel, where flow expansion can lead to condensation (see Fig. 20). This is a topic that has been investigated in detail in Refs. 93–95, where it was shown that condensation is often present in cases where the air is not heated, but its influence on plume structure and the aerodynamic loads are minimal.

The treatment of model blockage when simulating high altitude burns also needs to be considered. Conventional wisdom states that a model placed within a wall-bounded test section reduces the cross-sectional area, in turn changing the velocity field. In the case of a first-stage booster exhausting an expanding plume in the direction of the wind tunnel nozzle, it is the penetration length and the cross-sectional area of the plume that govern blockage and model placement. If a highly expanded plume with a large penetration length is located too close to the wind tunnel nozzle, it could cause such a large effective area reduction that choking occurs. It has also been shown that a single active nozzle produces a static plume, while three active nozzles have plume–plume interactions, resulting in oscillatory behavior between long and short penetration modes (see Fig. 21). This highlights the importance of supplementary design tools like CFD in the model development phase to determine the proper scale of the model as well as the correct mounting location in the tunnel. Finally, accounting for the plume shape and how it is influenced by scaling is critical for accurate plume modeling. This has led to the introduction of two additional similarity parameters for retro-propulsion experiments, namely the ambient pressure ratio and the thrust coefficient.⁹⁶ These parameters need to match the full-scale vehicle in order to obtain representative plume lengths and radial displacements.

Despite these challenges associated with model design, facility selection for re-useable first stage boosters is rather simple, as the trajectory falls well within the performance envelope of blow-down tunnels. This means that sufficient time is available to have the plume activated and achieve a steady state before taking measurements or varying the angle of attack.

Experimentally testing re-useable second stages shifts the challenges from the model design to facility selection. While in the past, re-useable spacecraft has only been associated with Earth orbit, the future is focusing on lunar and interplanetary travel, which significantly changes the simulation requirements in terms of both altitude and velocity. Manned vehicles returning from a lunar mission experience velocities up to around 12 km/s (see Fig. 12), while missions from Mars have been estimated to have Martian re-entry speeds in the range of 5.8–11 km/s, with Earth re-entry velocities in the range of 14–22 km/s.⁹⁷ These speeds mean that a large portion of the entry trajectory falls within the operating range of expansion tubes. An additional consideration for designing Martian entry experiments, or any other foreign atmospheres, is the additional challenges associated with replicating the gas compositions. For example, the Martian atmosphere is composed of ∼95% Carbon Dioxide (CO₂) and 5% Nitrogen (N₂) and Argon (Ar). This results in a lower speed of sound due to the reduced ratio of specific heats (∼1.29 on Mars vs 1.4 on Earth), which also influences flow characteristics such as shock density gradients and enthalpy. For all non-continuously run experimental facilities, exchanging the test gas mixture is generally possible but can require modifications to the tunnel infrastructure. For shock or expansion tubes, the process is as simple as charging the driven section with the desired gas.⁷³ For other facilities such as arc jet plasma wind tunnels, it has been shown that the exact composition of a Martian atmosphere cannot be realized due to unacceptable levels of cathode erosion.⁹⁸ Continuously run facilities draw ambient air, meaning that it is not possible to alter the test gas, and as such, the Mach number and dynamic pressure from the predicted Mars entry trajectory cannot be reproduced. Variable density facilities have the advantage of being able to alter Reynolds number and dynamic pressure, but these still fall short of replicating the flight envelope on Mars.⁹¹ 

The final point of discussion relates to flight testing and its importance for use as design verification. The nagging question associated with ground-based testing is always centered around the scalability of results to flight conditions and which facilities are best able to mimic real-world, full-scale aerothermodynamics. The unfortunate fact about flight test programs is that they are expensive. This mostly limits research institutions to conducting testing using vehicles flying sub-orbital trajectories. Sounding rockets are the launchers of choice, and due to their payload capacities, the size of the test articles is limited. Although larger scale models can be realized compared to wind tunnel models, they are still significantly smaller than full-scale manned vehicles. Despite the relatively small scales and Reynolds numbers, these types of flight tests provide engineers and scientists with a wealth of data for correlating numerical results and wind tunnel data with flight test data for a range of hypersonic specific technologies.^99–102 Unfortunately, in-flight measurements for operational launch vehicles are scarce, particularly now in the age of private space companies, which have no interest in or motivation to make these data available to the public. Some minor exceptions to this rule have appeared in the form of a collaboration between NASA and SpaceX that focused on understanding the aerodynamics and thermal environment of the Falcon 9 during retro-propulsion maneuvers, with a particular focus on the validation of numerical tools.^103,104 However, the publicly released information did not contain any useful data for outsiders attempting to validate their own numerical tools or experimental results. This overall unavailability of data is in stark contrast to the days of the Apollo or Space Shuttle programs, which were funded through taxpayer money and NASA, therefore, had the obligation to publish. As a result, there is a superabundance of literature relating to all facets of the Apollo and Space Shuttle designs, including experimental and flight results. Due to the simple and well-documented shape of the Apollo capsule, along with the wind tunnel and flight test data available, recreating experiments or performing code validations is easily achievable.^105,106 Despite the large repository of Space Shuttle flight results, including the critical field of hypersonic boundary layer transition,¹⁰⁷ the data cannot be replicated independently as the geometry is not available in the open literature. Recreating the geometry from freely available drawings and diagrams is possible, but it introduces unquantifiable uncertainties. Figure 22 presents an example of the available flight data pertaining to in-flight measurements, highlighting the hypersonic boundary layer transition from the Space Shuttle.¹⁰⁸ More recently, during the development of the Orion capsule, which is used to house the crew on the NASA Space Launch System, experimental data as well as geometrical details have been published,^109,110 which may form the foundation for the publication of future flight data.

In summary, ground-based testing is faced with multiple challenges associated with facility selection, instrumentation, and the scalability of results to flight conditions. However, given the scarcity of flight test data, results collected from wind tunnel tests typically form the foundation for the validation of numerical methods. In turn, this allows vehicle designers to build a basis of trust in high fidelity numerical tools, which are becoming more capable of accurately modeling complex problems related to re-useable spacecraft.

The development of simulation tools has so far been characterized by a continuous development toward higher fidelity and accuracy. Their history started with simple panel methods in the 1970s, followed by linearized and non-linear potential methods and numerical solution schemes for the Euler and Navier–Stokes equations, which first used statistical Reynolds averaged Navier–Stokes equation (RANS) methods for the treatment of turbulence and then later also scale resolution schemes up to direct numerical simulation (DNS). The main enabler is the constant rise of available central processing unit (CPU) resources, which corresponds to a million-fold increase during the last 20 years. This has been accompanied by algorithmic advances such as the introduction of higher order solution methods. However, low fidelity methods are still applied for pre-design applications, which require a fast exploration of large parameter spaces. A general trend in numerical simulation techniques is a constant evolution from empirical to physics based models. However, major challenges still persist such as the quantification of the effects of physical modeling errors, discretization errors, and epistemic uncertainties from unknown model parameters.

This also allows the derivation of closed-form expressions for the lift and drag of simple body shapes. The so-called modified Newtonian method uses an estimate for the stagnation pressure to scale the results [replaces the 2 in Eq. (3) with an estimate for the pressure coefficient at the stagnation point] and yields more realistic results. In the shock-expansion approach,¹¹³ a combination of local compression and expansion waves caused by the surface shape is used for an evaluation of the local surface pressure in a streamwise direction. This method works particularly well for 2D configurations. The tangent-wedge and tangent-cone methods¹¹⁴ use a similar approach. Here, the surface properties of any point on the vehicle surface are given by an equivalent wedge or cone flow that is constructed tangent to the point, with the tip being fixed at the vehicle's leading edge. This works for slender configurations where all resulting cone or wedge angles are below the limit of shock detachment.

An alternative to surface inclination methods for pre-design purposes is the application of blast wave theory.^25,115 This method is based on Haye’s equivalence principle:¹¹⁶ The hypersonic flow over a slender body is equivalent to an unsteady flow in one less space dimension. A blast wave is initiated by a point or line energy source, and the unsteady propagation in time is mapped to a steady propagation in the mean flow direction. Analytical solutions for blast waves can then be used to derive the aerodynamic properties of the mapped body.

On the next level of complexity are potential methods that are based on the compressible Euler equations. As for incompressible flows, the assumption of irrotational or potential flow greatly simplifies the calculation procedure. However, contrary to the incompressible case, the resulting potential equation is non-linear, and there is no general analytical solution method. A noteworthy approximate method is the slender body theory.^25,117 Here, the potential equations are linearized for small perturbations, and analytical solutions are possible. Famous results of this theory are the Sears–Haack body (minimum drag for given volume and length) and the von Kármán body (minimum drag for given length and base area).

Another established and widely used solution method for 2D compressible supersonic flows is the method of characteristics.¹¹⁸ The idea is to reduce partial differential equations to a family of ordinary differential equations along which the solution can be integrated. This method is applicable to the Euler equations for supersonic flows due to the hyperbolic nature of the governing equations. The result is that along Mach lines (characteristic lines), the sum of the Prandtl–Meyer angle (a function of the local Mach number and the ratio of specific heats) and the local flow angle is constant. This is the so-called Riemann invariant and represents the solution of the reduced ordinary differential equations along the set of characteristic lines. Hence, the problem is reduced to a non-linear algebraic system of equations. This model is widely used in the design of supersonic nozzles and general isentropic compression or expansion surfaces. It cannot treat discontinuities such as shock waves or embedded regions of subsonic flow.

The Boltzmann equation is theoretically valid for the entire range of Knudsen numbers, from free molecular to continuum flow. For mono-atomic gases at a steady state, it is a non-linear integro-differential equation for the six-dimensional probability density function of gas particle positions and momenta.^119,120 Additional dimensions are added in cases of unsteady problems or inelastic collisions, e.g., due to the exchange of kinetic energy with the internal energy of excited molecular states. Popular solution algorithms are based on stochastic particle methods like the Direct Simulation Monte-Carlo (DSMC) technique.¹²¹ The DSMC method employs a large number of Lagrangian particles over a computational domain. The evolution of the flowfield is represented by decoupled particle advection and probabilistic binary collisions. An important requirement of accurate DSMC simulations is the proper resolution of molecular scales: the time step size shall be smaller than the local mean collision time, and the grid cell size shall be smaller than the mean free path. These restrictions motivated the development of other particle-based stochastic methods that are more efficient for high density flows like the BGK (Bhatnager–Gross–Krook) or Fokker–Planck methods. Their main idea is to replace the collision operator in the Boltzmann equation by simplified expressions using a convection–diffusion (Fokker–Planck) or relaxation-to-equilibrium (BGK) ansatz.^122–124 The computational cost of these methods thus becomes essentially independent of the collision rate. The statistical noise of stochastic methods generally scales with the inverse Knudsen number and the inverse Mach number or inverse squared Mach number for shear stress and heat flux, respectively.¹²⁵ Hence, continuous approaches such as the discrete velocity method¹²⁶ are attractive as an alternative to particle based methods for low flow speeds and denser gases. However, in the case of high speed flows, stochastic methods require fewer CPU resources and do not suffer from the limitations of continuous methods to resolve distribution functions in strong non-equilibrium.

Chapman and Enskog developed a solution to the Boltzmann equation in the case of slightly non-Maxwellian distribution functions based on the assumption that the collision operator is dominant¹¹⁹ in the limit of high densities. The distribution functions used are perturbations of the Maxwell function. A zero-order perturbation results in the Euler equation, a first-order perturbation in the Navier–Stokes equation, and a second-order perturbation in the Burnett equations. Although the validity range of the Burnett equations extends to the higher Knudsen numbers compared to the Navier–Stokes equations, they are very difficult to solve, require ill-defined supplementary boundary conditions, and can potentially violate the second law of thermodynamics in arbitrary non-equilibrium situations.¹²⁷ Hence, there is very little practical application to date. The validity range of the different sets of governing equations with different solution methods is schematically summarized in Fig. 23.

By far, the most commonly used set of governing equations for the aerothermal analysis of space vehicles are the Navier–Stokes equations. They are generally valid for Knudsen numbers below around an order of magnitude of 10⁻², which corresponds approximately to an altitude of 80 km for the analysis in Fig. 11. The range of validity of Navier–Stokes can be slightly extended by the application of slip boundary conditions, which allow velocity and temperature jumps at walls.¹²⁸ However, this corrects only errors that are caused by surface boundary conditions and does not account for general deficiencies in the rarefied flow region such as the presence of non-physical barrel shocks that affect the back-flow predictions of flow expansions into the vacuum¹²⁹ or false predictions of continuum shocks on the leeward side of slender entry vehicles at high angles of attack, which result in an over-prediction of static pressure and hence an under-prediction of vehicle lift and drag. Yet, in some cases, Navier–Stokes equations have surprisingly good results, even for moderately large Knudsen numbers. One example is the internal structure of a shock wave¹³⁰ with an associated Knudsen number of 0.1. Note that the derivation of the Navier–Stokes equation using the Chapman–Enskog approach also provides rigorous expressions for transport coefficients such as viscosity, thermal conductivity, and diffusion coefficients, along with mixture rules for multi-component gases. The required input parameters are usually a set of binary collision integrals for the gas mixture under consideration. The expressions are obtained by identifying the resulting terms of the collision operator in the Chapman–Enskog solution with the corresponding gradient terms in the Navier–Stokes equations (e.g., the term in front of velocity gradients corresponds to the viscosity).¹²⁰ 

At present, the majority of applied aerothermal analyses of spacecraft rely on the numerical solution of the Navier–Stokes equations by numerous methods of CFD. CFD provides solution methodologies for the conservation equations of mass, momentum, and energy for viscous compressible flows, and the necessary building blocks for CFD-based solution techniques are schematically shown in Fig. 24.

The first and central building block concerns the numerical method (gray in Fig. 24). Generally, a discretization scheme is used to consistently transfer the governing continuous partial differential to a discrete system of algebraic equations that depends on a finite number of parameters. Consistency implies that the discretization process can be reversed through a Taylor series expansion to recover the governing equations. Then, an equation solver is used to find an approximate solution to the algebraic system. This solver needs to be stable which means that numerical errors tend to decrease during the solution process. A convergent method consists of a consistent discretization and a stable solver. In this case, the approximate and exact solutions converge in the limit of an infinite number of parameters. This concept is schematically outlined in Fig. 25.¹³¹ 

The discretization can either be discrete (the numerical solution consists of a set of point values or solution samples) or functional (the solution is expressed as a function depending on a finite number of parameters). Popular examples are finite difference methods for discrete discretizations and finite-element or spectral methods for functional discretizations. Often, hybrid discretizations are used in which discrete representations are combined with piece-wise functional representations in control volumes or elements. They are characterized by the fact that the functional representations are not continuous over control volume or element boundaries. Examples are the discontinuous Galerkin finite-element method¹³² or the finite volume (FV) method with gradient reconstruction inside control volumes.

There are a few general challenges related to CFD for hypersonic flows that concern numerical discretization and solution schemes. Flow fields are characterized by the presence of very large gradients. For example, the accurate representation of temperature gradients close to the wall that is required for heat load prediction leads to the requirement of a very tight wall normal grid spacing¹³³ in the order of 10⁻⁶ m, which also results in unfavorable large cell aspect ratios above 10⁴. Furthermore, strong shock waves can produce large errors when the shock is not aligned with the grid. A well-known example is the accumulation of errors in the stagnation region of a blunt body downstream of the strong bow shock, which leads to a complete breakdown of the numerical solution and meaningless results.¹³⁴ This is the so-called carbuncle phenomenon. Hence, one of the key concerns in the development of numerical techniques is to add dissipation for stability without adverse effects on the actual flow physics or the accuracy of the solution scheme.

At present, by far the most widely used schemes are finite volume methods. Alternative approaches such as discontinuous Galerkin appear to be a promising candidate for future solvers as they provide a higher order of accuracy on arbitrary elements. However, they still need further development and validation for hypersonic applications to find a solution for unresolved issues such as sub-element shock capturing.

The idea of finite volume (FV) approaches is to use flux functions to compute mass, momentum and energy fluxes across cell boundaries and construct the iterative cell update (residual) from their balance. FV methods are based on the integral formulation for conservation laws and are, therefore, particularly attractive for Navier–Stokes equations in the presence of discontinuities such as shocks. Generally, different approaches are used for the inviscid and viscous terms, as also schematically shown in Fig. 24. The inviscid fluxes are computed from reconstructed states at both sides of the cell boundaries. This reconstruction scheme uses available cell averages and determines the order of accuracy of the numerical method. Slope limiters are used to limit gradients during the reconstruction process and to stabilize the solution around discontinuities and in regions of steep gradients with associated local extrema. This adds additional numerical dissipation in critical regions, and the optimization of limiters is one of the primary concerns in the development of CFD methods. The fluxes across the cell boundaries are then computed from the states at both sides of the cell boundaries by approximate Riemann solvers, of which a wide range has been developed.¹³⁵ This approach captures the hyperbolic characteristics of Euler equations at Mach numbers above 1 and, therefore, avoids the non-physical upstream propagation of disturbances. However, a concern is the applicability in the limit of low Mach numbers, where standard upwind-based solution schemes suffer from accuracy deficiencies and usually require dedicated modifications.¹³⁶ The viscous terms are elliptic and dissipative by nature. They can simply be discretized with a computationally efficient central stencil without destabilizing the solution procedure.

As already outlined in Sec. III, hypersonic flows are energetic, and high temperatures occur. When the flow is brought to rest (e.g., at viscous walls) or passes through a shock wave, most of its significant kinetic energy is converted to internal energy, which leads to the excitation of vibrational and electronic energy modes and to chemical reactions. In many cases, the time scales of these processes are comparable to flow residence times, which requires the modeling of thermal relaxation and chemical kinetics (see also Fig. 11). As an example, assuming a calorically perfect gas without vibrational and electronic excitation and without endothermal dissociation reactions would result in a gas temperature in the shock layer in Fig. 10 of around 8000 K instead of the physically correct value of 3800 K. This results in the need for accurate models for high temperature thermodynamics, chemistry, and viscous transport coefficients, as shown by the red and blue boxes in Fig. 24. Generally, the fluid is treated as a chemically reacting mixture with different components or species. The mixture properties are evaluated based on the values of the individual species using mixture rules. For thermodynamics, this is the simple mass weighted sum of the individual contributions. The treatment of transport coefficients is more complex and is described below. The species thermodynamic properties are evaluated either from partition functions and spectroscopy constants (which enable the separate treatment of internal energy modes and the modeling of thermal non-equilibrium) or from curve fitted data. Curve fits are mostly applied for polyatomic molecules that are present, e.g., in the exhaust plumes of rocket engines. Here, the partition function approach, with the resulting need for a comprehensive set of spectroscopy constants, is often too cumbersome.

If thermodynamic and chemical relaxation time scales are small, the flow can be considered to be in equilibrium, and no thermal or chemical rate models are required. The fluid mixture can generally be described by two independent state variables, which also define the local mixture composition and the vibrational and electronic excitation states.

For both the thermodynamic properties and the equilibrium composition, which is a function of the chemical potentials of the participating species, the uncertainties in the available data are very small. An example is given in Fig. 26 for the distribution of specific heat at constant volume for equilibrium air. It compares the results from a partition function approach¹²⁰ and curve fitted data from NASA polynomials¹³⁷ and European Space Agency (ESA) lookup tables.¹³⁸ The results for the specific heat are not distinguishable. The local peaks caused by the occurrence of dissociation and ionization are also clearly visible.

Note that the effect of anharmonic vibration is often negligible in equilibrium situations. This is because it becomes important only at temperatures where the corresponding molecules are already dissociated, as shown by the curves for the nitrogen specific heat for harmonic and anharmonic oscillator assumptions in Fig. 26.

Chemical and thermodynamic relaxation rates are often in a similar order of magnitude as fluid time scales. In this case, chemical and relaxation kinetics have to be considered in the modeling approach. This is performed by additional transport equations for the mass conservation of mixture species and for the vibrational (and electronic) energies. The modeling happens in the source terms of those additional equations. Chemical non-equilibrium is modeled by the source terms in the species mass conservation equations, which represent the chemical production or destruction rates. In many cases, it is sufficient to calculate them using the law of mass action together with a suitable Arrhenius-based reaction mechanism. Here, detailed mechanisms can be used if pure air flows are considered. The presence of reacting exhaust gases in retro-propulsion scenarios of hydrocarbon-fueled rockets results in very large, stiff, and complex reaction mechanisms, and often simplified global or skeletal schemes are used. Despite the relatively simple approach, there are still uncertainties in the modeling of chemical non-equilibrium effects in air shock layers. An example is given in Fig. 27 ¹³⁹ which shows the shock standoff distance over the ratio of average energy exchange by chemical reactions and free stream energy flux on the stagnation streamline of a cylinder flow. Different popular air chemistry models (Park, Dunn and Kang, and Gupta)¹³⁹ are compared for air and nitrogen flow at different total enthalpies and still show a variation of predicted shock standoff in the order of 20%. Generally, the shock standoff scales with the dissociation state of the molecules. The endothermal reactions decrease the temperature, which increases the density in the shock layer and results in a smaller shock standoff.

Even larger sources of uncertainty are introduced by the modeling of thermal non-equilibrium effects, with the vibrational excitation of the molecules being the most important one. As previously stated, the standard approach is to solve additional transport equations for the vibrational energy of the considered molecules. The associated source terms are governed by vibrational–translational and vibrational–vibrational relaxation processes. Associated relaxation times are evaluated using semi-empirical correlations.¹⁴⁰ However, a comparison to detailed calculations from gas kinetic theory shows large uncertainties even for pure nitrogen plasmas.¹⁴¹ An appropriate modeling strategy also has to account for the strong coupling of vibrational excitation and dissociation chemistry. Highly excited molecules tend to dissociate with a higher probability, which is called preferential dissociation. Furthermore, the vibrational energy of molecules formed by recombination reactions tends to be greater than the equilibrium value. This may significantly impact local flow properties in the shock layer; nevertheless, the prediction of overall thermal loads at cold and catalytic walls is often insensitive to detailed relaxation phenomena.¹⁴² A more complex issue of thermal non-equilibrium modeling is the treatment of potential non-Boltzmann populations at large altitudes. This requires state-to-state kinetics with a very large set of transition rates. This approach often suffers from incomplete sets of collision cross sections, and due to the very heavy CPU demand, it is usually limited to generic 1D analysis of, for example, a stagnation streamline.¹⁴³ It is generally used to assess and calibrate simpler models, which are more applicable to large scale analyses of real configurations.

The last class of required rate models refers to gas-surface processes such as pyrolysis, catalysis, and ablation. A detailed discussion is left out of this paper. Ablative thermal protection systems are not attractive for re-useable spacecraft, and for catalysis it is often sufficient to consider fully or non-catalytic cases. This is performed by the appropriate boundary conditions in the species mass conservation equation (prescription of the local equilibrium composition for fully catalytic and vanishing wall normal gradients of the partial density for non-catalytic simulations). An excellent, detailed review of gas-surface rate processes is given by Candler.¹⁴⁰ 

Another important building block in Fig. 24 are the viscous transport models. Transport coefficients for momentum, heat, and mass transfer (viscosity, heat conductivity, and diffusivity) need to be approximated, and the assumption of a Newtonian fluid results in the linear stress–strain relationship. A widely used approach is the application of curve fits for the viscosity of the different species (e.g., Sutherland or Blottner) and the application of appropriate mixture rules (e.g., the Wilke rule).^137,144 The heat conductivity of the mixture is then either evaluated using a Prandtl number assumption or Eucken’s approximation, which takes into account different contributions from internal energy modes.¹⁴⁴ The diffusion coefficient is constant for all species due to mass conservation constraints from the application of Fick’s law and is computed using a Schmidt number assumption.¹²⁰ A more accurate approach is to use the Chapman–Enskog method.^119,120,145 It relies on data for binary collision integrals, which are available for all interaction pairs of a mixture. Comprehensive data are available for air, but in most cases, not for complex hydrocarbon mixtures in exhaust plumes. Because the evaluation of the mixture coefficients involves the solution of linear systems, which results in heavy computational loads in CFD simulations, the application of approximate mixture rules, e.g., the Yos-method,¹⁴⁴ is popular. The Yos formula is remarkably accurate for the viscosity of neutral plasmas and for the heavy particle contribution to heat conductivity.¹⁴⁵ The contribution of electrons can be computed separately using methods such as Devoto’s approach.¹⁴⁶ A similar problem is the treatment of diffusion, which involves the solution of the Stefan–Maxwell equations.¹¹⁹ In addition, here, approximations have been developed to reduce the computational effort. The most popular option is the Ramshaw model.¹⁴⁷ 

The uncertainty for viscous transport coefficients, which is still present for high temperature air, is illustrated in Fig. 28. It compares results for the viscosity and thermal conductivity of high temperature air for three popular modeling approaches, including Blottner curve fits and the Wilke mixture rule,¹⁴⁴ NASA-polynomial curve fits,¹³⁷ and data from the evaluation of binary collision integrals with the Yos mixture rule.¹²⁰ 

The last important building block of CFD methods in Fig. 24 refers to the modeling of turbulence. Turbulence is a three-dimensional broadband multiscale problem in space and time with a highly non-linear coupling between those scales. Generally, a deterministic description of the complete detailed velocity field is impossible. The most important effect for engineering applications is the additional diffusivity being introduced by the vortical motion of the flow in turbulent eddies. This results in an increase in momentum, heat, and mass transfer. Due to the existence of coherent structures, characteristic scales, and spatial and temporal correlations, turbulent flows are not fully random. A cascade process, which is present in all turbulent flows, continuously transfers turbulent energy from larger to smaller scales within the so-called inertial subrange. Dissipation to heat due to the effect of molecular viscosity occurs at the smallest eddies around the Kolmogorov scale, κ_K. The turbulent energy production is driven by the large eddies and their interaction with the mean flow at the so-called integral length scale, κ_l. Usually, the integral scale is very large compared to the Kolmogorov scale. The difference increases with increasing Reynolds numbers, which makes the simulation of the full spectrum increasingly difficult for high Reynolds flows typical of atmospheric entry problems. An overview of the turbulent spectrum and associated modeling approaches is shown in Fig. 29.

The main modeling approaches are statistical turbulence models based on the Reynolds averaged Navier–Stokes (RANS) equations, a mixture of resolved large scales and modeled small scales called large Eddy simulation (LES), and the resolution of all turbulent scales in a so-called direct numerical simulation (DNS).

DNS corresponds to a deterministic representation of all the turbulent scales of the flow, from the largest ones to the smallest dissipation scales (Kolmogorov scale). In that case, all the turbulent scales are associated with the resolved part, while there is no unresolved part. Due to the extremely high CPU demand, this technique is limited to very simple flow fields and low Reynolds numbers, which imply a reduced width of the turbulent spectrum.

LES techniques rely on a decomposition of the aerodynamic field between the large and small scales of the flow, the largest ones being directly resolved, while the effect of the small ones is represented through the use of a model. This is useful because, due to their much more homogeneous, isotropic, and self-similar nature, the small scales are considerably easier to model. The production of turbulence occurs through the interaction of the large scales with shear in the mean flow, and the potential anisotropic properties of the large scales are directly resolved. The computational time scales with the Reynolds number, Re, as t ∝ Re^2.5. This is due to the increasing size of the inertial subrange close to the wall in boundary layers due to the associated reduction of the dissipation length scale. To remedy this problem, hybrid models are developed, which principally aim at modeling the near-wall region rather than resolving it. The most popular approaches are detached eddy simulation (DES) and wall modeled LES. In DES, a statistical RANS model is used in the boundary layer, and wall modeled LES relies on the application of local wall models that exploit self-similarities in turbulent boundary layer profiles. Here, the computational time scales only as t ∝ Re^0.5, which highlights its usefulness for engineering applications.

The RANS framework is based on a purely statistical approach. The idea is to decompose the flow field into a mean and a fluctuating part and average the governing equations (Reynolds or Favre averaging of the Navier–Stokes equations). The resulting additional terms contain averages of fluctuating flow quantities and represent the momentum, mass, and heat transport phenomena caused by the effect of turbulence. The terms need modeling as the averaging does not provide any closed-form solutions. The most popular options are two-equation turbulence models, where one transport equation covers the turbulent kinetic energy and the other its dissipation rate. Additional modeling assumptions (e.g., linear relation of turbulent stress and strain rate tensor from the Boussinesq hypothesis, constant turbulent Prandtl and Schmidt numbers) then allow us to use these quantities for the modeling of turbulent mass, momentum, and heat transport in the Navier–Stokes equations. An excellent summary of turbulence modeling in CFD is given in the book by Wilcox.¹⁴⁸ 

Another significant challenge is the modeling of the laminar-to-turbulent transition process, which is an active research field on its own. Different approaches include semi-empirical models based on transport equations for characteristic quantities like intermittency and transition Reynolds number, models based on stability analysis, and the direct computation of instabilities in scale resolving LES and DNS frameworks. Due to the complexity of the problem, all approaches are still subject to major uncertainties. A detailed review of current methodologies is given by Pasquale et al.¹⁴⁹ 

Space missions may include atmospheric entries from hyperbolic orbital velocities or highly radiating exhaust plumes from rocket propulsion systems. The effects of gas radiation in the shock layer and wake flows may become significant at such flight conditions; hence, an accurate computational fluid dynamics analysis of these flows also requires adequate modeling of the local absorption and emission properties and the energy transport by radiation (solution of the radiative transport equation). The range of numerical techniques for the integration of radiation in entry vehicle simulations includes Monte Carlo models, which enable accurate and computationally efficient predictions of radiative energy transport in participating (absorbing and emitting) media for complex geometries. Alternatives are continuous methods like spherical harmonics or discrete ordinates. Re-entry-type flow properties are characterized by a large variation in local absorption and emission properties in the spectral range. Therefore, and due to the absence of general averaging rules, a fine spectral discretization is needed for practical computations. In some cases, narrow-band or opacity binning methods offer the potential to reduce the required resolution while maintaining a high level of accuracy.¹⁵⁰ 

An example of a combined analysis is shown in Fig. 30, where the convective and radiative heat fluxes on the front heat shield of the hypervelocity flight experiment FIRE II are plotted. In addition, the heat flux spectrum at the stagnation point clearly shows the footprint of bound–bound (atomic lines) and bound-free radiative processes.¹⁵⁰ 

In summary, the numerical simulation of high temperature flows around re-useable spacecraft presents a challenging multi-physics problem. It involves models for chemical reactions, thermal excitation processes, multi-component transport phenomena, turbulence and transition, as well as radiation. Many of these phenomena are too complex to be accessible for computations based on first principles, and simulation tools still rely heavily on empirical correlations or largely simplified models. This introduces significant uncertainty to the results, and a close interaction with ground based and flight experiments is needed to overcome these issues.

Current research in numerical techniques largely focuses on the reduction of model uncertainties.

Concerning turbulence, there is no approach that is capable of predicting all relevant phenomena and details for the entire range of flight conditions. RANS models continue to remain the standard for engineering applications. Although significant improvements in their predictive capability have been achieved in the past decades, their accuracy in predicting complex flow fields is still limited. This is the case for complex pressure gradients at shock and expansion wave systems, strong streamline curvature, separated flows, anisotropic turbulence (e.g., downstream of shock waves), and treatment of heavily cooled walls.¹⁵¹ LES methods demonstrated a dramatic accuracy improvement with minimal (or no) requirement for calibration of model constants. However, LES methods are still limited to low Reynolds numbers. Further research is currently focusing on the development of approximate models for the near-wall flow to enable simulations at higher Reynolds numbers. One example is detached eddy simulation (DES). Further progress in this methodology recently led to the development of different delayed-DES (DDES) or improved delayed DES (IDDES) schemes, which especially aim at improving the treatment of the RANS/LES interface zones.¹⁵² 

Some recent research aiming to improve RANS models focuses on shock turbulence interaction. A systematic deficiency of current turbulence models, which are based on RANS, is their inability to correctly predict the interaction of turbulence with shocks. This is because RANS models do not account for the unsteady motion or fragmentation of the shock wave within the interaction zone. These fundamental weaknesses also translate into more complex situations like the interaction of shock waves with a boundary layer where RANS models consistently fail to predict the heat flux, especially in the re-attachment zone.¹⁵³ Because the typical shock thickness is about one order of magnitude less than the smallest turbulent scale (Kolmogorov), the turbulence causes the shock to locally oscillate and, hence, the relative shock Mach number is reduced. Furthermore, entropy, vorticity, and sound waves are emitted from the oscillating shock and dissipate into turbulence amplification downstream. Typically, significant over-prediction of turbulent energy amplification occurs without dedicated adjustment of the applied turbulence model. The development of correction models for different RANS models has become possible due to the availability of recent DNS results.^154,155 The main idea is to develop a statistical correction for the coupling of shock movement and upstream velocity fluctuations.^156,157 Example results for the development of the shock normal and tangential components of the Reynolds stress tensor, R₁₁ and R₂₂, downstream of a Mach 6 shock wave are shown in Fig. 31. It can be seen that uncorrected RANS models can over-predict the turbulent amplification across a shock wave by several orders of magnitude.

The development of models to predict the laminar-to-turbulent transition is still at an early stage, and major uncertainties are present. Robust transport equation-based models suffer from many empirical correlations that need calibration. Models based on local stability analysis impose very high requirements on the quality of the baseline solution. Recent progress in experimental techniques such as the introduction of quiet supersonic tunnels and the further development of relevant measurement techniques for hypersonic applications such as temperature sensitive paint, laser differential interferometry, and micro-particle image velocimetry (PIV) form the foundation of further progress.^158–160 

Countless improvements can be performed in the area of physico-chemical modeling, which is an area of continuing research and, rightfully, the domain of physical chemists. Most uncertainty is introduced by thermo-chemical rate models and, most importantly, by the treatment of thermal non-equilibrium effects. Again, progress is limited due to the sparse availability of experimental data. Current approaches which are based on multi-temperature models are almost all calibrated to flows under compression in shock layers, and significant challenges still exist for expanding flows like those seen in wind tunnel nozzles. A multitude of unresolved problems are related to the influence of contamination on vibrational relaxation, for example, water vapor or the coupling of vibrational excitation and dissociation reactions. Due to the high CPU-demand and the limited availability of cross-section data, state-to-state models are still limited to generic academic applications.

Another area of active research is uncertainty quantification for numerical simulation techniques.^161,162 This covers several sources like the model due to its mathematical form, numerical solution and discretization schemes, or model parameters. Uncertainty quantification is the process of characterizing all major sources of uncertainty, quantifying their propagation within the simulation method and, therefore, their effect on analysis outcomes. Rigorous uncertainty quantification suffers from the large number of samples required, and recent research focuses to remedy this problem by either using surrogate models based on, for example, neural networks¹⁶¹ or improved propagation methods with reduced sample requirements such as the polynomial chaos expansion technique.¹⁶³ 

Currently, the field of machine learning is extremely dynamic and is showing that it can provide algorithms and means to learn from data without explicit mathematical models. A particularly successful concept for science and engineering applications is deep learning based on neural networks with multiple layers. These are able to learn simple relationships in the initial layers, and subsequent layers can combine this information to model more abstract behavior. Since many physical problems have a similar hierarchical structure, they can be efficiently modeled using neural network approaches.¹⁶⁴ There are several reviews about the potential applications of machine learning in computational fluid dynamics simulation tools.^165,166 Present applications include models for turbulence closure, feature extraction, dimensionality reduction, reduced order modeling, and error estimation. In the field of turbulence, neural network-based RANS closures correct the Reynolds stress tensor in complex flow situations, and LES models of small (subgrid) turbulent scales are improved, especially for regions of critical spatial discretization. Other applications include surrogate models for particularly stiff and CPU-demanding sub-models like combustion chemistry.^167,168 The challenge here is to enforce physical consistency and avoid temporal accumulation of errors and drift phenomena in both the physical models and mathematical algorithms.^169,170

Both algorithms and hardware for quantum computing are reaching a critical stage in their development and could potentially offer significant advancements in computational power and capabilities.¹⁷¹ Recent advances in the field of fluid mechanics include the successful application of quantum circuits as differentiable function approximators to solve the one-dimensional Navier–Stokes equations.¹⁷² However, other studies identify decreasing quantum advantages for large Reynolds-number (large geometrical scale) problems based on lattice-Boltzmann methods.¹⁷³ The potentials for CFD are also recognized in NASA’s CFD-Vision 2030 report.¹⁶² Nevertheless, quantum computing is still at an early stage, suffering from high error rates and the need for algorithmic developments. The extreme performance gains are limited to specific problems such as factorization. The transfer from laboratory experiments to robust technology for scalable applications still needs to be completed.

The environmental impact of space utilization is a subject of growing attention.¹⁷⁴ Due to the recent surge in re-entering debris and reuseable components, nitrogen oxides from re-entry heating contribute measurably to stratospheric ozone depletion. Rocket engine emissions occur throughout the entire flight trajectory at all altitudes, which results in a radiative forcing potential per unit mass emitted that is around 500 times greater than surface and aviation sources. Here, especially black carbon emissions are critical. The increasing ozone damage and climate effects are likely to give rise to regulation of the rapidly growing space industry. On top of that, there are considerable emissions from destructive entry (e.g., aluminum oxide). Satellite re-entries from mega-constellations are likely to become the dominant source of high-altitude alumina in the near future, with presently unknown effects.¹⁷⁵ The prediction of nitric oxide due to entering debris or re-useable components is relatively easy as the high temperature air chemistry is well understood. On the contrary, the effect and amount of soot and black carbon emissions are principally unknown, and no validated predictive tools exist. Modeling is a significant challenge; even for aero-engines, large uncertainties are still present,¹⁷⁶ and there is insufficient experimental data available. Due to the harsh environments in rocket plumes, measurements are extremely difficult.

Ground based testing of re-useable vehicles with active retro plumes has forced wind tunnel facilities to be operated in non-traditional ways, creating previously unseen challenges. Current research is focusing on accurate scaling of the retro plume as well as creating wind tunnel models capable of surviving both internal and external thermal stresses that result from replicating the correct exhaust gas composition and temperature. This is critical for understanding heat fluxes, for the structural design of the full-scale vehicle, as well as for the sizing of thermal protection systems. To date, there has been no such campaign that has successfully achieved this goal.

Advances in diagnostic tools for impulse facilities are primarily focused on overcoming the challenges associated with the short test times and high enthalpy flow conditions. The response times of diagnostic tools play an extremely important role in the ability to obtain meaningful data from these facilities. The field of free flight testing and optical tracking for drag, or in the case of models with onboard propulsion systems, thrust measurements, continues to evolve. Recently, there have been advances in the area of TSP, where the coating formulations have been refined to have acceptable response times. However, challenges arise under high enthalpy conditions, where the emission spectra of the paint overlap with the test gas luminosity, making data evaluation difficult. This particular problem is the subject of current research, as this technique can be used for the effective determination of heat fluxes as well as for diagnosing laminar to turbulent boundary layer transitions.

In the context of a multi-faceted approach to vehicle design, where numerical tools are used alongside experimental data, CFD solutions are superior to experimental datasets in terms of the volumetric resolution of the flowfield. The majority of experimental techniques, including spectroscopy and Schlieren imaging, are categorized as line-of-sight tools. This means that a single integral result is provided for the flowfield, which leaves some amount of uncertainty in the local conditions. Using Schlieren photography as an example, this can manifest in the form of unknown shock origins. Laser based techniques such as particle image velocimetry have been successfully used in high speed compressible flow environments, achieving volumetric and planar resolution of the flow field. Under extreme conditions of high enthalpy flow, issues can arise with effective particle seeding, meaning flow seeding systems will vary according to the facility test conditions.⁸⁹ Combining volumetric flow diagnostic tools with surface coating techniques in the future would significantly improve the ability of numerical validations to be effectively conducted.

Finally, the generation of flight data for validation purposes is largely influenced by the cost of test flights, which for the majority of research institutions is a severe limiter. As launch activity is predicted to surge over the coming decades, this may cause a trickle-down effect in making flight testing more affordable, creating an environment where the acquisition of flight test data is no longer as challenging as it has been in the past.

The aerothermal design of space vehicles is a challenging multi-physics task. In ground-based testing, it is impossible to fully replicate the flight conditions. The harsh flow conditions and short test times make the application of many diagnostic tools difficult. On the other hand, simulation methods require a wide range of models that still suffer from significant uncertainties in their input parameters (e.g., thermal relaxation), significant simplification assumptions (e.g., turbulence), or a large amount of empiricism (laminar–turbulent transition). Flight tests, while technically possible, are often limited by high costs, launcher availability, and publication restrictions. High fidelity simulation methods such as scale resolved turbulence models or state-to-state thermodynamics require large computational resources and are not yet applicable to comprehensive parametric studies.

Hence, for practical design tasks, a combined application of pre-design methods, high fidelity numerical simulation, and wind tunnel tests is required to manage uncertainties. Pre-design methods can be used to explore larger design spaces. Their results can be calibrated and assessed based on a combination of high-fidelity simulations and ground testing. Numerical simulation is also a useful means to extrapolate ground testing results to flight conditions.

The new concept of retro-propulsion based re-useable rockets also imposes major challenges. This involves complex chemistry in the plume interaction zone as well as highly unstable and fluctuating flow fields in the launcher wake, which require expensive time resolved investigations both in experiments and simulation.

The transition regime between rarefied and continuum flow remains critical for numerical aerothermodynamic analysis as the application limits of Navier–Stokes based simulation tools are reached and, yet, the density is still too large to apply DSMC economically. Here, ground based experiments are still inevitable. Recent progress in Fokker–Planck and BGK methodologies shows a promising potential to remedy this problem. Further current research focuses on the implementation of short response non-intrusive measurement techniques such as TSP, PSP, or absorption spectroscopy, which specifically aim at a better understanding of surface heating patterns and local flow conditions. On the numerical side, the introduction of discontinuous Galerkin finite-element methods for higher order accuracy and improvements in physico-chemical modeling, e.g., vibration–chemistry interaction or state-to-state kinetics, are likely to improve predictive capabilities in the near future.

The authors have no conflicts to disclose.

Sebastian Karl: Conceptualization (lead); Data curation (lead); Writing – original draft (lead). Tamas Bykerk: Conceptualization (supporting); Data curation (supporting); Investigation (supporting); Writing – original draft (supporting).

Data available on request from the authors.