Myriam Bastide

Robert Bernhard

Berthold Remi Kempf

Flavien Denis

Keywords: hypersonic flow, thermal heating, ablation, shock tunnel, CFD hypersonic flow, thermal heating, ablation, shock tunnel, CFD

After three years so quickly passed, the end of the PhD is finally in sight. It is now time to look back at how far we have come.

Résumé

La dernière décennie a vu l'essor de nouvelles applications dans le domaine du vol hypersonique. Le contrôle de la désintégration des satellites en fin de vie est devenu nécessaire pour les activités spatiales et la nouvelle course à l'armement a relancé le développement de missiles et projectiles hypersoniques. Dans ce contexte, l'extension du domaine de vol des véhicules hypersoniques nécessite de nouvelles études sur l'échauffement aérodynamique et sur l'ablation des matériaux.

La présente étude porte sur la détermination des profils de flux thermique et d'ablation par fusion pour des matériaux métalliques en vol hypersonique en atmosphère dense. Les profils de flux thermique et d'ablation sont mesurés en soufflerie à choc. Un protocole expérimental, combinant métal ablatif à bas point de fusion et visualisation optique, permet la détermination des profils d'ablation dans des écoulements hypersonique en quelques millisecondes. L'étude expérimentale met en évidence un profil caractéristique de l'ablation par fusion en écoulement turbulent, caractérisé par un maximum d'ablation localisé à 30 • du point d'arrêt.

L'étude numérique présente et valide le solveur ablationFOAM dédié à la prédiction de l'échauffement et de l'ablation en écoulement hypersonique. L'utilisation de ce dernier pour la simulation des essais menés lors de l'étude expérimentale démontre que le profil d'ablation caractéristique est la conséquence du développement de la couche limite turbulente en amont du maximum d'ablation. Il est montré également que l'ablation par érosion est négligeable par rapport à l'ablation par fusion. Enfin, l'utilisation des modèles de turbulence k -ω SST et k -ω SST γ -Re θ est également discutée.

Mots-clés : écoulement hypersonique, échauffement, ablation, tunnel à choc, CFD

Abstract

The last decade has seen the development of new applications in the field of hypersonic flight. The control of the disintegration of satellites at the end of their life has become necessary for the space activities and the new arms race has re-started the development of hypersonic missiles and projectiles. In this context, the flight domain extension of hypersonic vehicles requires new studies on aerothermal heating and on the thermal ablation.

The present study focuses on the determination of wall heat flux and ablation profiles for metallic materials in dense hypersonic flow. The heat flux and ablation profiles are measured in a shock tunnel. An experimental protocol, combining low melting point ablative metal and optical visualisation, allows the ablation measurement with a hypersonic impulse facility within the milliseconds. The experimental study highlights a characteristic profile for the ablation by melting in turbulent flow. It is characterised by an ablation maximum localised at 30 • from the stagnation point.

The numerical study presents and validates the ablationFOAM solver dedicated to the prediction of the aerothermal heating and of the ablation by melting in hypersonic flow. The simulation of the tests carried out during the experimental study demonstrates that the characteristic ablation profile is the consequence of the development of the turbulent boundary layer upstream of the point of ablation maximum. Furthermore, it is shown that the ablation by erosion is negligible compared to the ablation by melting. The use of the turbulence models k -ω SST and k -ω SST γ -Re θ is also discussed.

Zusammenfassung

Im letzten Jahrzehnt werden neue Anwendungen im Bereich des Hyperschallflugs entwickelt. Die Kontrolle des Zerfalls von Satelliten am Ende ihrer Lebensdauer ist für die Raumfahrtaktivitäten notwendig geworden und das neue Wettrüsten hat die Entwicklung von Hyperschallraketen und -geschossen wieder in Gang gebracht. In diesem Zusammenhang erfordert die Erweiterung der Flugdomain von Hyperschallflugkörpern neue Studien zur aerodynamischen Erwärmung und zur Hyperschallablation von Materialien.

Die vorliegende Studie konzentriert sich auf die Bestimmung von Wärmestrom-und Schmelzabtragungsprofilen für metallische Werkstoffe im Hyperschallflug in dichter Atmosphäre. Die Wärmestrom-und Ablationsprofile werden in einem Stoßtunnel gemessen. Ein experimentelles Protokoll, das niedrigschmelzendes ablatives Metall und optisches Messverfahren kombiniert, ermöglicht die Messung der Ablation in Hyperschallströmungen in Millisekunden. Die experimentelle Studie zeigt ein charakteristisches Profil der Ablation durch Schmelzen in turbulenter Strömung, das durch ein Ablationsmaximum gekennzeichnet ist, das 30 • vom Staupunkt lokalisiert ist.

In der numerischen Studie wird der ablationFOAM Solver vorgestellt und validiert, der für die Vorhersage von Erwärmung und Ablation in Hyperschallströmungen entwickelt wurde. Die Simulationen der in der experimentellen Studie durchgeführten Versuche zeigen, dass das charakteristische Ablationsprofil die Konsequenz der Entwicklung der turbulenten Grenzschicht vor dem Punkt des Ablationsmaximums ist. Außerdem wird gezeigt, dass die Ablation durch Materialabtragung im Vergleich zur Ablation durch Schmelzen vernachlässigbar ist. Die Verwendung der Turbulenzmodelle k -ω SST und k -ω SST γ -Re θ wird ebenfalls diskutiert.

Stichwörter: Hyperschallströmung, thermische Erwärmung, Ablation, Stoßtunnel, CFD v Que de dangers, que de fausses routes dans l'investigation des sciences ! Par combien d'erreurs, mille fois plus dangereuses que la vérité n'est utile, ne faut-il point passer pour arriver à elle ! Le désavantage est visible : car le faux est susceptible d'une infinité de combinaison; mais la vérité n'a qu'une manière d'être. Qui est-ce d'ailleurs qui la cherche bien sincèrement ? Même avec la meilleure volonté, à quelles marques est-on sûr de la reconnaître ? Dans cette foule de sentiments différents, quel sera notre critère pour en bien juger ? Et ce qui est le plus difficile, si par bonheur, nous la trouvons à la fin, qui de nous en saura faire un bon usage ? Jean-Jacques Rousseau, Discours sur les sciences et les arts, 1750 What a variety of dangers surrounds us! What a number of wrong paths present themselves in the investigation of the sciences! Through how many errors, more perilous than truth itself is useful, must we not pass to arrive at it? The disadvantages we lie under are evident; for falsehood is capable of an infinite variety of combinations; but the truth has only one manner of being. Besides, where is the man who sincerely desires to find it? Or even admitting his good will, by what characteristic marks is he sure of knowing it? Amid the infinite diversity of opinions where is the criterion by which we may certainly judge of it? Again, what is still more difficult, should we even be fortunate enough to discover it, who among us will know how to make right use of it? Jean-Jacques Rousseau, Discourse on the Arts and Sciences, 1750

List of Symbols

a = Sound velocity [m • s -1 ] = Thermal diffusivity [m 2 • s -1 ] b = b-factor [K -1 ] C H = Stanton number [-] C x = Axial force coefficient [-] C p = Specific heat capacity at constant pressure [J • kg -1 • K -1 ] D = Deformation rate tensor [s -1 ] D = Model diameter [m] f v = Volumetric force [N • m -3 ] d = Surface displacement [m] e = Specific total energy [J • kg -1 ] = Specific energy [J • kg -1 ] = Thermal effusivity [W • K -1 • m -2 • s -0.5 ] H f = Enthalpy of fusion [J • kg -1 ] h = Altitude [km] = Enthalpy [J • kg -1 ] = Plank constant [m 2 • kg • s -2 • K -1 ] I = Electrical current [A] I = Identity tensor [-] k = Turbulence kinetic energy [J • kg -1 ] = Boltzmann constant [J • K -1 ] k = Body force [P a • m -1 ] K n = Knudsen number [-] L = Characteristic length [m] Le = Lewis number [-] M = Mach number [-] m = Mass [kg ] ṁ = Mass flow rate [kg • s -1 ] n = Surface normal [m] p = Pressure [P a] P 0 = Freestream total flow power [W • m -2 ] P r = Prandtl number [-] Q = Source term [[Φ] • kg • m -3 • s -1 ] = Thermal heat [W ] q = Heat flux vector [W • m -2 ] q = Heat flux [W • m -2 ] q = Thermal source term [W • m -3 ] q = Dynamic
[K -1 ] = Phase fraction [-] = Thermal diffusivity [kg • m -1 • s -1 ] β = Ballistic coefficient [kg • m -2 ] Γ = Diffusion coefficient [kg • s -1 ] γ =
W • m -1 • K -1 ] θ = Position angle [ • ] µ = Dynamic viscosity [P a • s] ν = Kinematic viscosity [m 2 • s -1 ] = Fundamental vibrational frequency of a molecule [H z] ξ = Thermal effusivity [J • m -2 • K -1 • s -0.5 ] = Transformed coordinate [-] ρ = Density [kg • m -3 ] σ = Stefan-Boltzmann constant [W • m -2 • K -4 ] = Mechanical

Objectives and Outline of the thesis

In the introduction of [START_REF] Van Driest | The problem of aerodynamic heating[END_REF], van Driest formulates the following four questions regarding the hypersonic ablation:

• "What is the rate at which heat enters the surface of a high-speed vehicle?

• Where does the heat go?

• What happens when the heat gets there?

• How can the designer live with it?"

The present work shall try to answer these questions for the ablation by melting, and more specifically for the ablation of metallic parts of vehicles or projectiles. The present document constitutes a synthesis of the work done at the French-German Research Institute of Saint-Louis (ISL) in the group for Aerodynamics, Measurements & Simulations (AMS) between 2018 and 2022. The topic of ablation was investigated in the 60's at ISL through hyperballistic tunnel experiments and analytical studies. The experimental methods developed and the results obtained are compiled in the PhD thesis of Luneau [START_REF] Luneau | Contribution à l'étude des phénomènes aérothermiques liés au vol hypersonique décéléré[END_REF] and in several ISL reports. The goal of the current work is to present an experimental method for the study of ablation effects in the hypersonic shock tunnel and finally to develop a numerical model for the simulation of aerothermal heating and the relevant ablation effects.

This document is organized in three parts:

• The Part I introduces the topic hypersonic flight and liquid ablation. In Chapter 1 a brief history of the hypersonic flight and the TPS is given. An introduction to the modern issues concerned with these topics is provided. The state of the art is presented in Chapter 2.

• The Part II describes the experimental studies performed in the framework of this thesis. First, the wall heat flux measurements on a sphere in the hypersonic shock tunnel are presented in Chapter 3. This part presents the wall heat flux measured along a sphere in a dense flow and without ablation. The second experiment is presented in Chapter 4 describes the measurement of several ablation profiles in the shock tunnel.

The ablation profiles and ablation start times are presented. Benefits and drawbacks of this new measurement method are discussed.

• The part III is dedicated to the numerical modelling. The Chapter 5 describes the numerical model used in the numerical solver for hypersonic flow and ablation. The verifications and validations of the solver are presented in Chapter 6. Then the Chapter 7 analyses the CFD results obtained by the simulations of the corresponding wall heat flux measurements. Finally the simulations of the ablation experiments are presented in Chapter 8.

• The Part IV presents the conclusion of the thesis.

Part I

Introduction

Chapter 1

Overview of hypersonic flight and thermal protection

1.1 Historical perspective on thermal protection "...re-entry...is perhaps one of the most difficult problems one can imagine...It is certainly a problem that constitutes a challenge to the best brains working in these domains of modern aerophysics...one can use mass-transfer cooling, thus keeping the temperature under the allowed limit..."

Theodore von Kármán, 1956 [START_REF] Von Kármán | Aerodynamic heating -The temperature barrier in aeronautics[END_REF] The problem of aerothermal heating has emerged with the beginning of space research in late the 1950s. After the success of the first orbiting of the Soviet satellite Sputnik 1 on the October 1957, the challenge of returning the Earth arose. Indeed, during the reentry into the Earth's atmosphere, the high velocity of the reentry vehicle combined with compressibility effects increases the temperature of air surrounding the vehicle up to several thousands of Kelvin. Thus, the reentry vehicle must withstand the high thermal load it encounters during the several minutes of reentry.

The reentry problem can be solved by two design solutions. First, the reduction of the maximum wall heat flux. Indeed the peak wall heat flux decreases as the surface radius increases ( [START_REF] Van Driest | The problem of aerodynamic heating[END_REF] & [START_REF] Fay | Theory of stagnation point heat transfer in dissociated air[END_REF]). That is why the vehicle surfaces facing the flow, also called the stagnation surfaces, are designed to be blunt and rounded. It is to note that the shape of the stagnation surface strongly impacts the overall aerodynamics of the reentry vehicle, since the bow shock is formed in the stagnation region and generates an entropy layer along the flying body. Hence the radius of the curvature is also dependent on the flight mission. For example, a reentry capsule dedicated to human space flight and designed to decelerate slowly with a low entry path angle features a large radius of curvature. Figure 1.1a shows a typical profile of a manned reentry capsule. This keeps the wall heat flux as low as possible. On the contrary, the flight mission of an ICBM reentry vehicle requires a short reentry duration, and hence a high reentry velocity. This is achieved with a steep reentry angle in combination with a high ballistic coefficient (β = m/(S r e f C x )). Maximizing this parameter implies the use of a blunted cone, as shown in Figure 1.1b with a low nose tip radius, which is only limited by its thermal resistance.

The second design solution for the reentry problem lies in covering the stagnation surfaces with a Thermal Protection System (TPS). Its function is to isolate the payload and the structure from the heat load and to manage the thermal energy introduced into the heat shield. Two kinds of TPS can be distinguished and are described first: the heat sink and the ablative heat shield. [START_REF] Woollard | Aerodynamic and aerothermal TPS instrumentation reference guide[END_REF] The heat sink absorbs the wall heat flux the vehicle is exposed to. It is composed of a massive metallic part, which heats up during the reentry and also melts partially. This technology was used exclusively for TPS of warheads. As an example, the Air Force Mk I was designed with a 500kg copper heat sink while the Navy Mk I and Air Force Mk II used a sink made of beryllium. A beryllium heat sink had also been considered initially for the Mercury program but was eventually abandoned in favour of an ablative heat shield due three main disadvantages: the heat sink mass was judged prohibitive due to the weight of the massive metallic part, the very high residual temperature of the heat sink after reentry which resulted in a fire hazard, and the manufacturing difficulties of the beryllium parts. The last use of a beryllium heat sink was for scientific purpose on 22 May 1965 during the Fire II mission. For the measurement of the wall heat flux, three beryllium heat sinks separated by phenolic layers and used as calorimeters were placed in the front of the reentry capsule, as shown in Figure 1.2. When the first heat sink in contact with the flow heats up and then begins to melt, it is jettisoned and the next virgin heat sink faces the flow. This way, the heat sink temperature can be used to measure the wall heat flux during the reentry. While being the only TPS system based on metal heating and melting, the heat sink is nowadays no longer in use.

The second type of TPS is the ablative heat shield. Unlike the heat sink, it is not designed to absorb and store the heat but rather to consume it through chemical reactions, sublimation and gas blowing. The aerothermal heating induces depolymerisation, sublimation and pyrolysis of the heat shield material. These reactions are endothermic and thus have a cooling effect. The resulting gaseous products are blown into the boundary layer. As a result, the boundary layer thickness increases which decreases the temperature gradient at the wall and therefore also decreases the wall heat flux. The gas blowing also transports the heat from inside the heat shield into the flow, thus it thermally isolates the heat shield from the surrounding air flow. The research focused on Teflon, which only sublimates, as well as on composite materials. For example Gagarin's spacecraft Vostok (1961)(1962)(1963) was equipped with a heat shield made of asbestos fibre impregnated with resin and the Mercury spacecraft (1958)(1959)(1960)(1961)(1962)(1963) was protected with an head shield composed of aluminium honeycombs covered by fibre glass impregnated with phenolic resin. Ablative TPS have been used for all reentry space capsules and also remain the only TPS employed for the high speed entry into other planetary atmospheres.

The 1970s saw the design of reusable space orbiters: the US Space Shuttle Orbiter and the soviet Buran orbiter, which led to the development of reusable TPS, i.e. a non ablative heat shields. Here, the protection is based on the radiative cooling of the TPS material. The thermal energy is re-emitted in the form of thermal radiation into the surrounding environment. The TPS material should be heat resistant, with a high melting and sublimation point, and also should have favourable thermal radiation properties. For the Space Shuttle Orbiter, a ceramic-coated carbon-carbon composite protected the nose and the leading edges. Reusable TPS are effective at moderate thermal heating and are considered for reentry from Low Earth Orbit (LEO).

In parallel to the research on ablative and reusable TPS, alternative thermal protections have been investigated by the aeronautics and turbomachinery industries. During the cold war, the need for supersonic reconnaissance aircraft flying faster then Mach 3, such as the Lockheed A-12 and the SR-71 Blackbird, led to the use of new metallic alloys, which provide good heat resistance and mechanical strength. For example, Titanium alloy was used for the external structure of these aircraft. In the turbomachinery industry, the temperature increase at the first stage of the gas turbine necessitated the development of supper-alloys, which are commonly nickel-based, such as Inconel.

Finally to be mentioned is the technique of active cooling, which is often used for the thermal protection of rocket nozzles and turbine blades. In case of internal cooling, a coolant fluid is convected within the blade structure, whereas for external cooling, a cold gas is injected into the boundary layer, thermally isolating the wall surface from the hot flow. Injection is achieved through methods such as film cooling or transpiration cooling. By combining a nickel-based superalloy and an active cooling, the actual first stage turbine can work in hot gas up to 1800K at take-off, thus several hundred Kelvin above its material melting point.

Modern issue for the TPS at hypersonic speed

As shown in the section 1.1, ablative and reusable TPS have been used for space reentry since the 1970's. Extensive studies have investigated the complex phenomena which occur simultaneously around and inside the TPS. Special effort has been paid for the modelling of ablation in porous media, such as for the PATO program [START_REF] Lachaud | A generic local thermal equilibrium model for porous reactive materials submitted to high temperatures[END_REF] and also for the modelling of the ablation-radiation coupling occurring during high speed reentries [START_REF] Erb | Investigation of Galileo Probe Entry Heating with Coupled Radiation and Ablation[END_REF]. The surface recession rate during high speed reentry remains difficult to predict, as also proven by the flight data of the Galileo probe during its jovian entry. The stagnation point recession was overestimated by 30% while the recession at the frustum region was underestimated by up to 50%. Improvement of the knowledge of the ablation and oxidation processes will allow the TPS to be adequately designed, thus reducing its mass for the benefit of the payload mass.

Complementary to the heat shield ablation, the ablation of materials developing a liquid layer becomes a significant issue for two hypersonic applications, namely the mitigation of space debris risks and the design of future hypersonic weapon systems. Regarding the space debris mitigation, the recent regulations ( [3] & [START_REF]Loi n°2008-518 du 3 juin 2008 relative aux opérations spatiales[END_REF]) require the mitigation of the reentry risks for impact on people and property and thus recommend the controlled de-orbiting of LEOsatellites and rocket debris and also their total destruction during reentry. One of the main risk comes from metallic parts, such as fuel tanks or rocket nozzles, which are large and heat resistant enough such that a significant fraction of them is not destroyed during the reentry. For these, both the ablation by melting and the oxidation reaction play a major role for the total mass loss and also for the aerodynamic shape changes.

On the battlefield, a current trend is to increase the flight speed of missiles and projectiles. The flight domains for these new hypersonic applications are presented in Figure 1.3.

The development of hypersonic missiles is motivated by the requirement to keep the weapon efficient despite a contested airspace and thus to be faster than the enemy's defensive systems. Two hypersonic missile categories can be defined: Hypersonic Cruise Missiles (HCM) and the Hypersonic Glide Vehicles (HGV). The HCM are a further evolution of current supersonic cruise missiles. They are designed for a flight speed between Mach 5 and 10 at an altitude below 50km. They accelerate with a rocket booster and then use a scramjet for the cruise phase. Only Russia has officially claimed the first operational use of a HCM in March 2022 with the nuclear-capable air launched Kinzhal [START_REF]Kinzhal" missile threat[END_REF]. The real operational performance remains unclear, however Russia definitely has the capability to perform HCM flight tests, as demonstrated by the test flight of the second russian HCM project, the submarinelaunched 3M22-Tsirkon [START_REF] Dahlgren | Russia launches Zircon hypersonic missiles from submarine. Missile Threat[END_REF] in October 2021. In the 2010's, France has started studies for the AMSP/A next generation with the PEA CAMOSIS and Prométhée [START_REF] Boulaud | L'avenir des forces nucléaires françaises. Rapport d'information n°668[END_REF]. A main technical issue associated with HCM development is the scramjet design with the challenge of controlling the internal chemical reaction kinetics. Moreover the missile structure must sustain the heat load and protect the payload, such as warhead and electronics, from high temperatures.

Similar thermal load issues arise for the HGV. This weapon system aims to replace ballistic warheads of MRBM and ICBM. Its gliding and manoeuvring capabilities make its de-tection and its interception more difficult. The HGV flight domain is quite different from the one of a HCM. It is launched by a ballistic missile and is then released at high altitude, finally gliding to the target. During the glide phase, the HGV has no propulsion and controls its trajectory by the use of gravity and aerodynamic forces only. The HGV capabilities depend on the launcher performance (ICBM, MRBM), with the glide speed between Mach 10-20. Although it can carry both nuclear and conventional warheads, it appears to be limited to nuclear strikes due to the high system cost and the low precision on target (probable CEP ≈ 100m). Russia, China and the USA have already performed HGV flight tests. The operational use of the Russian Avanguard [START_REF] Robinson | From Sänger to Avangard -hypersonic weapons come of age[END_REF] will most likely be available within the next years whereas the chinese DF-17 system has already been declared operational. American HGV are still at the demonstrator level, but several projects are conducted simultaneously such as the US Navy-Conventional Prompt Strike or as the US Air Force AGM-183 Air Launched Rapid Response Weapon. On the European side, only France develops a HGV demonstrator V-MAX for which the first flight is expected in 2022. HGV will be protected from the thermal load with a TPS. The hypersonic gliding phase will last longer than a classical ballistic reentry and thus the long duration heating will be more severe than usual.

Finally to be mentioned, another hypersonic application lies within future hypersonic projectiles. They constitute a new type of ammunition system mainly dedicated to conventional and electromagnetic naval guns [START_REF] Hundertmark | A scenario for a future European shipboard railgun[END_REF]. Electromagnetic systems reach a high muzzle velocity, which extends the maximum range up to 185km for the HVP projectile [START_REF]HVP Hypervelocity Projectile. Datasheet[END_REF]. In addition the terminal velocity is high enough to destroy armored targets by kinetic effects only. The initial acceleration is obtained with a railgun or with a reinforced powder gun. Hypersonic projectiles and their associated launchers are still under development, even if the demonstrator level is reached as demonstrated by the HVP of BAE systems for the USA, the Pilum project for Europe, and the multiple railgun projects in China, Turkey and India. Their commissioning is yet not planned because short range missiles are actually considered more flexible. However recent research on medium calibre, high rate-of-fire railguns shows a potential use for the defence against anti-ship missiles [START_REF] Gallant | Design Considerations for an Electromagnetic Railgun Firing Intelligent Bursts to Be Used Against Antiship Missiles[END_REF]. The HVP flight particularity is that the maximum flight speed corresponds to the muzzle velocity, so it is reached at low altitude in the dense atmosphere. That results in high thermal heating during the first seconds of the flight. For cost and mechanical strength reasons, we can expect that the HVP will not be equipped with a heat shield and will be composed mainly of metallic materials, which results in strong thermal heating. While their external aerodynamics can be altered by liquid ablation, the internal conduction can lead to overheating of the payload.

Finally, it is important to underline that the majority of the known military hypersonic projects are offensive weapon systems and only the USA have publicly announced a research program for a HGV interceptor within the Glide Breaker programm [START_REF] Saylor | Hypersonic Missile Defense: Issues for Congress[END_REF]. Point defence interceptors are commonly used by but their speed is actually around Mach 3-4, thus a new generation of hypersonic ABM seems therefore needed for the defence against HCM and HGV threats.

Chapter 2

Theory and state-of-the-art of liquid ablation

In contrast to the ablation by pyrolysis and gas blowing and to the ablation by surface oxidation, the ablation by melting has been the subject of only a few studies in the past. This is due to the fact that the heat sink TPS was early abandoned for aerospace applications and that the melting phenomenon is only a secondary effect in modern TPS. In this chapter, a review of the phenomena associated with the ablation by melting and a literature survey is presented.

Hypersonic flow field

Figure 2.1 shows a typical hypersonic flow around a sphere. Characteristic for this type of flow is the proximity of the bow shock to the stagnation region. The region between the shock and the sphere surface is called the shock layer. Within this region, the gas temperatures and pressures reach high levels. Conversely the gas velocity is low compared to the freestream velocity, resulting in a region of subsonic flow centred on the stagnation point. Depending of the freesteam conditions, the boundary layer can become turbulent on the forebody part of the sphere. For the case of an incompressible flow around a sphere, the boundary layer transition is expected to occur for a Reynolds number Re D larger than around 3 • 10 5 -4 • 10 5 . An increase of the Reynolds number shifts the transition forward. The boundary layer is finally fully turbulent for a Reynolds number above 10 6 [START_REF] Oertel | Strömungsmechanik: Methoden und phänomene[END_REF]. The nature of the boundary layer is important for the ablation at low altitude, since the local wall heat fluxes and local wall shear stresses are much higher for a turbulent boundary layer than for a laminar one.

Ablation due to excess heating primarily occurs on all stagnation regions. In the case of the ablation by melting, the material facing the flow heats up and finally melts. It generates a liquid layer of melted material, which is sheared off by the flow and is convected downstream along the surface of the flying body. The rear parts of the body are not directly ablated, however the hypersonic flow pattern includes expansion and recompression regions, especially in the wake flow, such that the liquid layer can re-solidify in the expansion regions [START_REF] Moody | Tungsten and Molybdenum Ablation Modeling for Re-entry Applications[END_REF].

Conservation equations

The conservation equations describe the hypersonic flow ( [START_REF] Ferziger | Computational methods for fluid dynamics[END_REF], [START_REF] Moukalled | The Finite Volume Method in Computational Fluid Dynamics: An Advanced Introduction with OpenFOAM® and Matlab[END_REF] & [START_REF] Anderson | Hypersonic and high-temperature gas dynamics[END_REF]). A general form for a conservation equation of a scalar quantity φ is given by: For a non reactive, compressible flow, the conservation equations are the following:

∂ρφ ∂t t r ansi ent t er m + ∇ • ρφ u convec t i on t er m = ∇ • Γ φ ∇φ d i f f usi on t er m + Q φ sour ce t
• Mass conservation equation:

Dρ D t = ∂ρ ∂t + ∇ • ρ u = 0 (2.2)
• Momentum conservation equation:

D ρ u D t = ∂ ρ u ∂t + ∇ • ρ u ⊗ u = ∇ • T + k (2.3)
Where the stress tensor T can be split in two terms such that the pressure is defined as:

p = - T xx + T y y + T zz 3 = - t r (T ) 3 (2.4)
And the viscous stress tensor τ is:

τ = T + p I (2.5)
The gas considered in this study is a Newtonian fluid, so the viscous stress tensor τ is a function of the strain rate ∇ u only:

τ = µ ∇ u + (∇ u) T + λ (∇ • u) I (2.6)
With the Stokes' friction law, the bulk viscosity λ is set to: λ = 2 3 µ. The momentum conservation equation can be written as:

D ρ u D t = ∂ ρ u ∂t + ∇ • ρ u ⊗ u = -∇p + ∇ • τ + k (2.7)
• Energy conservation equation in terms of total specific energy:

D ρe D t = ∂ ρe ∂t + ∇ • ρe u = -∇ • p u + ∇ • uτ + ∇ • q + q (2.8)
In the case of a chemically reacting flow, a species-specific mass conservation equation should be added for each species.

Characteristic parameters

Three characteristic parameters describe the main properties of a non-reacting hypersonic flow: the Mach number, the Knudsen number and the Reynolds number. Additionally, the flow enthalpy, the stagnation pressure and the velocity gradient the stagnation point also determine the magnitude of surface heating and thus ablation intensity.

The Mach number, M = u/a, is the ratio of the flow velocity with respect to the sound velocity. It is a measure of the gas compressibility. The onset of the hypersonic flow regime is not clearly defined but commonly estimated between Mach 4 and 5, where the thermal effects become noticeable.

The Knudsen number, K n = λ/L, indicates whether the flow should be considered dense or rarefied. Below approximately 0.03, the flow is a continuous medium and the Navier-Stokes equations with slip-free wall conditions are valid. Above 0.03, non-continuum flow effects appear, such as wall slip effects or non-equilibrium molecular velocity distributions. From 0.2 onwards, the Boltzmann equation must be used ( [START_REF] Moss | Direct simulation of transitional flow for hypersonic reentry conditions[END_REF]) to describe the flow physics. All flows subjected to this work are dense flows and hence treated as continuum flow.

The Reynolds number, Re = ρuL/µ, quantifies the ratio of inertial to viscous forces. It is also an indication of the local turbulence state in the flow, e.g. if the boundary layer is fully laminar, fully turbulent or in a transitional state. For a Blasius boundary layer, the transition to turbulent flow is expected around a critical Reynolds number of Re c = ρUx c /µ = 5•10 5 . For a hypersonic flow the transition occurs over a wide range of critical Reynolds from Re c = 10 5 to 10 8 ([100] & [START_REF] Anderson | Hypersonic and high-temperature gas dynamics[END_REF]). The prediction of the critical Reynolds, i.e. the location of the laminarturbulent transition, is one of the most challenging problems for the modern aerothermodynamics. Indeed it is a function that depends on many parameters related to the fluid flow, such as the Mach number, the Reynolds number, the temperature, the pressure and the velocity gradients. It is also dependent on the type of hypersonic vehicle, such as its shape or its surface roughness. Figure 2.2 shows the dependence of the heat transfer coefficient on the Reynolds number. The thermal heating can be several times higher for a turbulent boundary layer compared to a laminar one.

High-temperature gas properties

Across the bow shock, the pressure and temperature increase drastically. The bulk kinetic energy of the flow is first converted into translational and rotational energies, and then the gas molecules begin to vibrate. If the energy is high enough, the electrons transit from their ground state to higher energy states. The energy transfer rates between all energy modes depend on the flow density, velocity and temperature and effectively equilibrate the flow over time, such that the translational, rotational, vibrational and electronic temperatures tend towards the same value. This final state is called thermodynamic equilibrium. The thermodynamic state of the flow is important since it relates the internal energy distribution of the gas to its temperature. (2.10) In addition, the gas molecules dissociate and react chemically. Similarly to thermodynamic equilibrium, a chemical equilibrium can be defined, which is also dependent on the temperature and pressure. Reaction rate coefficients have been proposed in the literature for hot air and nitrogen and can be found in [START_REF] Park | Nonequilibrium hypersonic aerothermodynamics[END_REF] or [START_REF] Dunn | Theoretical and experimental studies of reentry plasmas[END_REF].

In the present work, nitrogen flow is used as test gas and the maximum stagnation temperature in the experiments stays below 3200K . Therefore the translational and rotational are fully and the vibrational modes are partially excited. Furthermore, a weak, non-measurable effect of molecular dissociation and no effect of electronic excitation are expected. The thermodynamic and transport properties of nitrogen are expressed as a function of the temperature. The polynomial coefficients can be found in [START_REF] Mcbride | Coefficients for calculating thermodynamic and transport properties of individual species[END_REF] and are valid up to 6000K . The coefficients for nitrogen are given in Table 2.1 for the heat capacity and in Table 2.2 for the transport variables.

C p = R a 0 + a 1 • T + a 2 • T 2 + a 3 • T 3 + a 4 • T 4
(2.11)

N 2 A B C D E F G H For 300K < T < 1000K
0.604 -0.436E2 -0.884E3 0.190E1 -0.943 0.123E3 -0.118E5 -0.107 For 1000K < T < 5000K 0.651 0.285E2 -0.167E5 0.152E1 0.651 -0.151E3 -0.137E5 0.218E1

Thermal heating

The thermal heating of a hypersonic vehicle is caused by the heat transfer between the flow and the vehicle surface. It is a convection problem with the addition of thermal radiation and chemical reactions. For a chemically reacting flow with radiation, the total thermal conductivity in the flow is expressed as:

q = q cond ,g + q d i f f ,g + q r ad ,g (2.14)
The term q cond ,g quantifies the heat conduction and defines the energy flux caused by the temperature gradient. It is given by:

q cond ,g = -λ ∇T (2.15)
The term q d i f f ,g is the energy flux caused by the diffusion of all species in the flow. It is given by:

q d i f f ,g = i ρ i u i h i (2.16)
The last term q r ad ,g is the energy flux caused by the gas radiation in the shock layer. It is given by [START_REF] Luneau | Contribution à l'étude des phénomènes aérothermiques liés au vol hypersonique décéléré[END_REF]:

qr ad ,g = σT 4 0 ≈ 2.44 • 10 -23 R ρ ∞ ρ 0 σT 10 0 (2.17)
Below 10000K , the contribution of the thermal radiation to the total heating can be neglected, since the radiation flux is negligible compared to the energy advected by the flow [START_REF] Anderson | Hypersonic and high-temperature gas dynamics[END_REF].

In the present work, q d i f f ,g ,w and q r ad ,g ,w are set to zero since no chemically reacting flows are considered and the flow radiation is neglected. Finally, the total thermal conductivity is equal to: q = q cond ,g = -λ ∇T (2.18)

From Equation 2.18, the wall heat flux at the vehicle surface is given by: qw = -λ g ∂T ∂n w,g (2.19) Once the wall heat flux penetrates the vehicle surface, it is distributed between the internal thermal conduction, the surface radiative cooling and the phase change caused by ablation: qw = qcond,w,s + qr ad ,w,s + qabl a,w,s (2.20) With the heat conduction in the solid at the wall given by: qcond,w,s Stagnation wall heat flux q w, 0 (W m 2 ) q w, 0, Oertel q w, 0, Sutton q w, 0, Seiler q w, 0, Smeets

= 1.2kg • m -3 & R = 0.01m
It is dependent on the 4th power of the surface temperature, so the radiative cooling can be neglected if the surface temperature is low enough, as for the current work.

The term qabl a,w,s appears only if the wall surface melts and is ablated. It is a function of the ablation velocity and the latent heat. In case of the ablation by melting, the latent heat of fusion is used and the wall temperature will not exceed the melting temperature of the material.

qabl a,w,s = ρ s v a ∆H f ,s (2.23)

Stagnation wall heat flux

The estimation of the stagnation wall heat flux has been an important issue since the beginning of space exploration. Thus its prediction is a key factor for the design of hypersonic vehicles. In a laminar flow, the stagnation point is expected to match the position of maximum wall heat flux on the vehicle surface. Several methods from analytical calculations and empirical correlations have been proposed for the estimation of the stagnation wall heat flux.

A classical form for the stagnation wall heat flux correlation is derived from the Fay & Riddel theory [START_REF] Fay | Theory of stagnation point heat transfer in dissociated air[END_REF] and is expressed as:

q0 = C ρ α v β R -0.5
(2.24)

The coefficients C , α and β are dependent on the correlations chosen. Some values are given in Table 2.4. The evolution of the stagnation wall heat flux in function of the flight velocity is computed in Figure 2.3 with a density ρ = 1.2kg • m -3 and a vehicle nose radius of R = 0.01m. It is independent of the freestream temperature since the freestream enthalpy is negligible compared to the total enthalpy. As long as the surface temperature is negligible compared to the recovery temperature, it does not influence the stagnation wall heat flux.

The correlations are derived from boundary layer analysis. Examples of such an analysis can be found in [START_REF] Van Driest | The problem of aerodynamic heating[END_REF] for flate plate and blunt body problems, [START_REF] Fay | Theory of stagnation point heat transfer in dissociated air[END_REF] and in [START_REF] Lees | Laminar heat transfer over blunt-nosed bodies at hypersonic flight speeds[END_REF] for flows including chemical reactions at thermal equilibrium or non equilibrium. For a blunt body problem, van Driest [START_REF] Van Driest | The problem of aerodynamic heating[END_REF] proposed the following expression for the wall heat flux along a sphere with a laminar boundary layer: For a laminar boundary layer on a sphere in thermal equilibrium, Fay & Riddell [START_REF] Fay | Theory of stagnation point heat transfer in dissociated air[END_REF] proposed the following expression for the stagnation wall heat flux: q w,0 = 0.76P r -0.6 ρ w µ w 0.1 ρ e µ e 0.4 1 + Le

q w,0 = 0.763P r -0.

Local wall heat flux

For the wall heat flux along a sphere, Luneau proposed the following correlation ( [START_REF] Luneau | Contribution à l´étude de la répartition du flux de chaleur par convection s´exerçant sur une maquette de révolution en vol hypersonique. Formules pratiques[END_REF], [START_REF] Luneau | Contribution à l'étude des phénomènes aérothermiques liés au vol hypersonique décéléré[END_REF] & [START_REF] Luneau | Déformation d´un corps métallique par ablation. Équations de base conduisant à la réalisation d'une ablation théorique sur machine à calculer[END_REF]):

qw = C 1 R n ρ N ∞ u M ∞ f 1 (θ) h 0 -h w h 0 -h ∞ (2.28)
With: The evolution of the wall heat flux repartition function f 1 (θ) is shown in Figure 2. 4. It shows that the wall heat flux in a laminar boundary layer along a sphere monotonously decreases, whereas the wall heat flux in a turbulent boundary layer increases up to a maximum located at the angle 36 • . Moreover, the ratio between the turbulent wall heat flux to the laminar wall heat flux is proportional to the 1/3 power of the Reynolds number ( [START_REF] Luneau | Contribution à l'étude des phénomènes aérothermiques liés au vol hypersonique décéléré[END_REF] & [START_REF] Michel | Frottement et flux de chaleur turbulents en écoulement bi et tridimensionels[END_REF]).

f 1 (θ) = sin 1-n (θ) cos m (θ) θ n (2.29)

Liquid ablation

Liquid layer model

Figure 2.5 illustrates the phenomena that occur during the ablation by melting. Analytical models for the liquid layer can be found in several publications ( [START_REF] Bethe | A Theory for the Ablation of Glassy Materials[END_REF], [START_REF] Chang | Analysis of the aerodynamic ablation of a metal sphere[END_REF], [START_REF] Luneau | Contribution à l'étude des phénomènes aérothermiques liés au vol hypersonique décéléré[END_REF], & [START_REF] Duffa | Ablative thermal protection systems modeling[END_REF]). Initially developed for the ablation of silica, it is assumed that the viscous force contribution in the liquid layer is much higher than the inertial effect. The liquid layer is therefore modeled as a creeping flow and more specifically as a Couette flow: the tangential velocity at the solid/liquid interface is zero and non zero at the gas/liquid interface. Due to the low viscosity of the liquid metals, the creeping flow assumption is not fully justified for the ablation of metal, however this assumption is convenient for the simplification of the boundary layer equation.

The mass conservation in the liquid layer is expressed by:

∇ • u = ∂ (r u) ∂s + ∂ (r v) ∂y = 0 (2.30)
With the Couette flow assumption, the tangential velocity u follows a linear repartition along the wall normal direction:

u u e = 1 - y δ e (2.

31)

The integration of Equation 2.30 with Equation 2.31 gives:

ṁ f -ṁv = ρ l 2r ∂ r ρ l u l ∂s (2.32)
Assuming the sublimation to be negligible at the liquid-gas interface, the ablation mass flux is:

ṁ f = ρ l 2r ∂ r ρ l u l ∂s (2.33)
The momentum conservation equation along the y-axis in the liquid layer is:

ρ l C p l v ∂T ∂y = ∂ ∂y λ ∂T ∂y (2.34)
Introducing ∂T /∂x = -1/λ l q, Equation 2.34 is integrated along the y-axis and gives the ratio of wall heat flux to liquid-gas heat flux as:

ql qw = 1 - 2 3 C p l λ l ṁ f δ l + 1 6 C p l ρ l δ l u l λ l ∂δ l ∂s (2.35)
And:

T l = T f + qw δ l λ l 1 - 1 4 C p l λ l δ l ṁ f + 1 12 C p l ρ l u l λ l ∂δ l ∂s (2.36)
In a Couette flow, the shear stress is constant inside the liquid layer: The analytical solution for the ablation by melting is obtained by solving the following system of equations with five equations and the five unknowns δ l , u l , ṁ f , ql and T w :

τ = u l µ δ (2.
                   ṁ f = ρ l 2r ∂(r ρ l u l ) ∂s τ l = u l µ l δ l ql qw = 1 -2 3 C p l λ l ṁ f δ l + 1 6 C p l ρ l δ l u l λ l ∂δ l ∂s T l = T f + qw δ l λ l 1 -1 4 C p l λ l δ l ṁ f + 1 12 C p l ρ l u l λ l ∂δ l ∂s ṁ f = ql ∆H f (2.38)
A method for solving the system of equations and approximate solutions are proposed in [START_REF] Luneau | Contribution à l'étude des phénomènes aérothermiques liés au vol hypersonique décéléré[END_REF].

Effect of the oxidation

The ablation by oxidation, also called ablation by combustion, occurs when oxygen reacts with the wall surface or the ablation layer. Similar to a combustion reaction, it is an exothermic reaction, thus the ablation speed is generally enhanced in presence of oxidation reactions. The ablation profiles of duralium spheres ablated in nitrogen or in air [START_REF] Luneau | Contribution à l'étude des phénomènes aérothermiques liés au vol hypersonique décéléré[END_REF] show the influence of the test gas on the ablation. These differences are shown in Figure 2.6. The oxygen diffuses trough the ablation layer and oxidises the ablation products. The oxidation is enhanced in case of a turbulent liquid layer, since it increases the oxygen diffusion. As a result, the oxidation increases the wall heat flux, especially downstream of the location of laminar-turbulent transition. The ablation in this region is much more intense than without oxidation. As consequence, the conical forebody shapes of ablated spheres obtained in air are much sharper than those in nitrogen [START_REF] Luneau | Contribution à l'étude des phénomènes aérothermiques liés au vol hypersonique décéléré[END_REF].

Almost all metals have an oxidation layer at their surface. The melting point of the oxide layer is generally higher than those of the pure metal. Bridges & Pinchok [START_REF] Bridges | An experimental investigation of the re-entry behavior of metals and alloys[END_REF] showed with

) (u 0 = 5200m• s -1 & ρ ∞ = 1.7kg • m -3 )
pre-oxidised models made of titanium, aluminium, stainless steel, Inconel X and beryllium that only a thick oxidation layer in case of titanium or stainless steel models delays the ablation onset. They observed also that spallation and cracks on the oxidation layer negate the protective effect.

Effect of the erosion

The ablation by erosion is a secondary phenomenon as part of the liquid ablation process. It is caused by the wall shear stress induced by the flow viscosity. This effect is not related to any thermal phenomenon but is rather a purely mechanical effect. Very few experiments quantified the contribution of erosion to the total ablation. As example, for a duralium sphere, Luneau [START_REF] Luneau | Contribution à l'étude des phénomènes aérothermiques liés au vol hypersonique décéléré[END_REF] estimated the contribution of the erosion to 15% of the total ablation.

Experimental facilities

The ablation is a complex phenomena, in which fluid mechanics, thermodynamics, chemistry, heat transfer and solid state physics interact with each other. All these disciplines involve different time and space scales. Consequently, different types of experimental facilities are required for the study of each aspect of the ablation process. Furthermore, the facility design parameters, such as type, size, and available input power determine the maximal experiment duration and the flow enthalpy.

Low enthalpy wind tunnel

Low enthalpy wind tunnels are classically used for the aerodynamic characterisation of hypersonic vehicles. Typical test times range from several seconds to a continuous operation. The test gas can be preheated, so the stagnation temperature can reach several hundreds of Kelvins. These types of facilities generally do not reach the high temperatures required for the ablation of metal. However, ablation studies are possible by using substitutive ablation materials, such as wax [START_REF] Stock | Hypersonic low temperature ablation an experimental study of cross-hatched surface patterns[END_REF], water ice [START_REF] Silton | Ablation onset in unsteady hypersonic flow about nose tip with cavity[END_REF], dry ice ( [START_REF] Kohlman | Experiments on the use of dry ice ablating windtunnel models[END_REF] & [START_REF] Callaway | Photogrammetric measurement of recession rates of low temperature ablators in supersonic flow[END_REF], naphthalene ( [START_REF] Combs | Technique for studying ablation-products transport in supersonic boundary layers by using PLIF of naphthalene[END_REF], [START_REF] Combs | Heat-shield ablation visualized using naphthalene planar laser-induced fluorescence[END_REF]& [START_REF] Agostinelli | Investigation of Hypersonic Flow in the VKI H3 Wind Tunnel: from Facility Characterization to Boundary-Layer Interaction over Low-Temperature Ablators[END_REF]) or camphor [START_REF] Stock | Hypersonic low temperature ablation an experimental study of cross-hatched surface patterns[END_REF]. Water ice and dry ice are classically used for recession rate measurements in supersonic flow, since the first melts and the latter sublimates. Wax and camphor allow the visualization of ablation surface patterns such as cross-hatched surface patterns. The fluorescence property of naphthalene allows the tracking of the ablation products in the flow using Laser-Induced Fluorescence (LIF) methods.

Plasma torch wind tunnel

Plasma torch experiments are commonly used for the investigation of material degradation of objects subjected to high wall heat flux ( [START_REF] Rashis | Qualitative Measurements of the Effective Heats of Ablation of Several Materials in Supersonic Air Jets at Stagnation Temperatures up to[END_REF], [START_REF] Bridges | An experimental investigation of the re-entry behavior of metals and alloys[END_REF], [START_REF] Steg | Hypersonic ablation[END_REF], [START_REF] Winter | Characterization of ablation product radiation signatures of PICA and FiberForm[END_REF] & [START_REF] Duan | Ablation of C/SiC-HfC composite prepared by precursor infiltration and pyrolysis in plasma wind tunnel[END_REF]). The test durations, typically from several seconds to minutes, allow long duration ablation experiments, which provide useful results for the validation of ablation models. Investigations of chemical reactions, pyrolysis, gas blowing and mass loss can be conducted in such facilities. The operating performance expressed through mass flow rate or fluid temperature are set by the energy and power available to these facilities. Usually, high Mach number flows can be generated, however with low flow density. On the contrary a operation under dense flow conditions is possible if the Mach number is reduced. Due to these limitations plasma torches are not well suited for ablation studies in dense flow, i.e. at low altitude. The area of interest is generally restricted to the stagnation point region.

High enthalpy shock tunnel

High enthalpy shock tunnels are not classically used for ablation experiments since the typical test duration is of several milliseconds only. However ablation phenomena can be investigated with a low melting point ablator, such that the material is ablated within the test duration. The flows generated are realistic in terms of Reynolds and Mach numbers, so the aerothermodynamics around a model of a flight vehicle can be reproduced and investigated by applying similarity laws.

Free flight

Free flight experiments comprise all the tests carried out with hypersonic models, demonstrators, aircrafts or spacecrafts in hyperballistic facilities or during experimental rocket launches and atmospheric reentries.

Ablation data obtained during reentries are scarce but very valuable. Indeed, they serve as reference data for other ablation studies. Most of the real flight data concerns the ablation of the heat shield of reentry capsules. For instance, temperature and char-thickness were meseasured during the Apollo flights AS202 and AS501 Apollo 4 [START_REF] Curry | Apollo ablator thermal performance at superorbital entry velocities[END_REF] and AS502 Apollo 6. Thermal radiation and liquid ablation recession rate were measured during the Fire II flight. Heating data for a complex body can be found for example in the NASA report summarizing the flight data of the Columbia Space Shuttle Orbiter during its five first flights ( [START_REF] Williams | Columbia: The first five flights[END_REF], [START_REF] Williams | Columbia: The first five flights[END_REF], [START_REF] Williams | Columbia: The first five flights[END_REF], [START_REF] Williams | Columbia: The first five flights[END_REF], [START_REF] Williams | Columbia: The first five flights[END_REF] & [START_REF] Williams | Columbia: The first five flights[END_REF]).

Ablation tests with sounding rockets have been also performed, for example for the ablation measurements of a Teflon nose cone at a Mach number between 9 and 13 [START_REF] Hopko | A Free-Flight Investigation of Ablation of a Blunt Body to a Mach Number of 13[END_REF].

A hyperballistic tunnel is a ballistic range for free flight tests of hypersonic objects. The acceleration is commonly achieved with a two stage light gas gun. The free flight distances range from several tens to hundreds of meters. It allows the ablation study in a dense atmosphere, generally the amtospheric density. The benefits of such facilities are the excellent flow similarity with real reentry flows. Moreover the model moves in a quiet atmosphere, so that the turbulence upstream the bow shock is low. In the case of a closed test tunnel, the density and atmospheric composition are set to match with the required altitude or planetary atmosphere. Because the model moves at several kilometres per second, measurement techniques need to be fast measurement methods, such as X-ray flash-radiography or optical velocity measurements. In [START_REF] Chang | Analysis of the aerodynamic ablation of a metal sphere[END_REF], comparisons between analytical and free fight data are presented. Luneau presents the measurement methods and the ablation results obtained within the ISL hyperballistic tunnel in [START_REF] Luneau | Contribution à l'étude des phénomènes aérothermiques liés au vol hypersonique décéléré[END_REF] .

Numerical modelling

Most of the recent studies focused on the numerical modelling of the ablation process. The liquid ablation phenomenon is often considered as a secondary phenomenon and is rarely modelled. In addition the transport properties for the high temperature liquid layer are still uncertain, resulting in uncertainty about the effective wall heat flux passing through the liquid layer. Numerical issues are still unsolved for the modelling of liquid ablation. The liquid layer thickness is very small compared to the overall domain around the hypersonic vehicle. As a result, the direct computation of this layer using CFD models is still a challenge while keeping the computation duration acceptable. Obviously the shape change caused by the ablation can affect the gas flow around the vehicle significantly, which further alters the wall heat flux and eventually also the ablation speed.

Two categories of numerical modelling tools are described in the literature. The first category includes spacecraft-oriented ablation codes. They are low-to-medium fidelity tools, designed for the analysis of a complete system. They are dedicated to the fast prediction of heating, ablation and disintegration for satellites or hypersonic vehicles. For that purpose, they are based on correlations and on a fast solution of simplified conservation equations. Several such codes can be cited as those used in Europe for reentry problems: PAMPERO (CNES), SCARAB (ESA), FAST/MUSIC (ONERA) and DEBRISK (CNES) ( [START_REF] Omaly | DEBRISK, a tool for re-entry risk analysis[END_REF], [START_REF] Annaloro | Space debris atmospheric entry prediction with spacecraft-oriented tools[END_REF], [START_REF] Annaloro | Aerothermodynamics modelling of complex shapes in the DEBRISK atmospheric reentry tool: Methodology and validation[END_REF] & [START_REF] Annaloro | Rebuilding with PAMPERO of destructive hypersonic tests on honeycomb sandwich panels in the T-117 wind tunnel[END_REF]).

The second category includes the high fidelity tools based on CFD or Direct Simulation Monte Carlo (DSMC) solvers. The simulation results are obtained by either solving the Navier-Stokes equations or by statistically emulating rarefied flows described by the Botzmann equations. Here, the physics of the ablation phenomena have to be modelled. These tools provide better quality results than the low-to-medium fidelity tools but are slower and require a significant amount of computational resources. Navier-stokes based solvers are generally used for simulations of up to 80km of flight altitude, assuming continuum flows. The DSMC method and the Boltzmann equations are used for low density or rarefied flow simulation above 60km of altitude. Furthermore these codes generally include modules for the high temperature gas effects such as thermochemical effects, non equilibrium effects, multi-temperature models and radiation. Theses code are generally dedicated to hypersonic flow problems but not specialized for the simulation of ablation. Relevant tools are for example hy2Foam, dsmcFoam ( [START_REF] Casseau | A two-temperature open-source CFD model for hypersonic reacting flows, part two: multi-dimensional analysis[END_REF], [START_REF] Casseau | Hypersonic simulations using open-source CFD and DSMC solvers[END_REF] & [START_REF] Casseau | An open-source CFD solver for planetary entry[END_REF]), MISTRAL [START_REF] Annaloro | Space debris atmospheric entry prediction with spacecraft-oriented tools[END_REF] and NSMB [START_REF] Vos | ARD Post Flight Analysis[END_REF]. For specialized tools including ablation, the PATO code can be cited ( [7] & [START_REF] Lachaud | Porous-material analysis toolbox based on OpenFOAM and applications[END_REF]). It is an OpenFOAM based toolbox dedicated to the ablation in porous media such as the ablation of reentry heat shields.

Part II

Experimental investigation

Two experimental campaigns are conducted as part of this thesis. The first aims to measure the incoming wall heat flux along a spherical model instrumented with thin-film sensors. The obtained wall heat flux profiles allow the investigation of the flow field around the model and are used for the validation of the numerical CFD solver. The second experiment aims to investigate the liquid ablation phenomenon in a shock tunnel. A measurement method dedicated to impulse facilities is developed for ablation profile measurements in only 2 milliseconds.

All experiments are conducted within the high-enthalpy shock tunnels STA and STB in the "Aerodynamics, Measurements and Simulations" group (AMS) of the French-German Institute of Saint-Louis (ISL). A schematic and pictures of the facilities are presented in Figure 2.7 and 2.8. These facilities are able to generate hypersonic flow from Mach 3.5 to Mach 14 up to a total enthalpy of 8M J •kg -1 . Categorized as hypersonic impulse facilities, their steady flow test duration does not exceed several milliseconds.

Each shock tunnel is composed of a shock tube, a nozzle, a measurement section and a dump tank. The shock tube generates the total pressure and temperatures required for hypersonic flow creation. It is divided by a steel membrane into two sections of 100mm inner diameter, respectively a driver gas tube, also named high-pressure part, and a test gas tube, named driven part. The high pressure part measures 2.5m in length and is filled with a mixture of hydrogen and nitrogen up to 500bar . The driven part is 18.4m long and contains the test gas, which is nitrogen in the present work. The steel membrane is calibrated to rupture at a defined pressure, which depends of the required flow conditions. When the membrane breaks, the high-pressure part and the test gas tube are suddenly brought into contact, which generates a shock wave that propagates into the test gas, heating and compressing it. On reaching the driven part end wall, the shock wave is reflected backward, compressing the test gas a second time. The temperature and enthalpy reached after the second compression correspond to the stagnation conditions of the hypersonic flow.

The shock tube is separated from the nozzle inlet by a thin plastic membrane, which is generally a 50µm thick. The high pressure in the shock tube breaks this membrane and thus the test gas enters the nozzle and successively expands and accelerates. The model is fixed to a sting in the measurement section and is mounted about 65mm downstream of the nozzle end. Before the flow starts, the pressure in the measurement section is about 5Pa. Optical accesses for the visualizations is provided on each side of the test section.

After 1ms the hypersonic nozzle flow reaches a steady state. This time is considered as the starting point for the measurements. The shock wave velocity, the stagnation, or nozzle inlet, pressure and the nozzle exit static pressure are measured, such that the freestream flow conditions can computed from these data using a one-dimensional shock-tunnel code. The typical steady flow duration for this work is 2ms, during which the total pressure varies by less than 15%.

The facilities are described in more details in [START_REF] Gnemmi | Shock tunnels at ISL[END_REF].

Solid Edge

Echelle :

1: 100

Wall heat flux measurements on sphere

This experimental campaign aims to measure the wall heat flux on a sphere under hypersonic flow conditions. Two aims are pursued: firstly, the measured wall heat fluxes will be used for the assessment of the numerical wall heat flux computed with the CFD solver ablationFOAM. Secondly, the wall heat flux distribution will provide information about the flow close to the model surface, i.e. information about the boundary layer.

The wall heat flux measurements are performed in the shock tunnel STB for two Mach numbers and different flow densities. The data measurements were performed in the context of the bachelor thesis of Mielke [START_REF] Mielke | Wärmestromdichtemessung an Flugkörpern im Hyperschall-Hochenthalpie-Stoßwindkanal. Bachelorarbeit, Institut[END_REF]. The sensors are thin-film thermometers produced by the Shock Wave Laboratory of the RWTH Aachen University [START_REF] Olivier | Thin film gauges and coaxial thermocouples for measuring transient temperatures[END_REF]. The sensor body consists of a ceramic cylinder of diameter 2.3mm composed of zirconium dioxide (Z r O 2 ), a good thermal insulator. A thin layer of CHAPTER 3. WALL HEAT FLUX MEASUREMENTS ON SPHERE nickel, below 1µm, covers the sensor face, which faces the flow. The sensor is able to measure a surface temperature change of less then 0.1K with a response time of several microseconds. The temperature at the sensor surface is measured via the electrical resistance change of the nickel layer given by: R = R 0 (1

Experimental setup and processing

+ α∆T s ) (3.1)
By supplying the sensor with a constant current, the voltage variation during the test is given by Ohm's law:

∆u(t ) = ∆R(t )I (3.2)
Thus the temperature variation of the nickel layer is computed as:

∆T s (t ) = ∆u(t ) αR 0 I (3.3)
With the assumption of constant thermal properties within the sensor body, the wall heat flux is obtained by solving the heat conduction equation in the sensor body. Because the shock tunnel experiment is short, the temperature deep inside the sensor core does not change and the one-dimensional semi-infinite heat transfer theory can be used. Thus, the heat flux at the sensor surface is linked to the temperature variations through ( [START_REF] Smeets | Berechnung der Wärmeströme aus Endwandgrenzschicht und laminarer Stoßgrenzchicht zur Eichung von Oberflächenthermometern mit Stoßwellen[END_REF] and [START_REF] Kinnear | Design, calibration and testing of transient thin film heat transfer gauges[END_REF]):

qs (t ) = ξ π ∆T s (t ) t + 1 2 t 0 ∆T s (t ) -∆T s (τ) (t -τ) (3/2) d τ (3.4)
Introducing Equation 3.3 in 3.4:

qs (t ) = ξ παR 0 I ∆u(t ) t + 1 2 t 0 ∆u(t ) -∆u(τ) (t -τ) (3/2) d τ (3.5)
Equations 3.4 and 3.5 are exact only if the thermal properties of the sensor body are not temperature dependent. Walenta [START_REF] Walenta | Analog networks for high heat-transfer-rate measurements[END_REF] proposed a correction to this model to take the temperature dependence into account. Assuming a linear relationship of the thermal properties with temperature, a corrected temperature can be defined as:

λ = λ 0 (1 + b∆T s ) (3.6) ∆T s,cor r (t ) = ∆T s (1 + b∆T s ) (3.7)
Finally the corrected temperature is inserted into Equation 3.4. The wall heat flux across the sensor surface is computed by solving this equation with an algorithm developed at ISL. The initial time for the heat flux evaluation corresponds to the time when the sensor signals are stabilized after flow start. All constants, excepted the R 0 from Equation 3.3 and the bcoefficient in Equation 3.7, are set through a calibration of each sensor by measuring the wall heat flux on the surface of a shock tube end wall after shock impingement. The initial resistance R 0 is measured prior to each run, since its value can significantly vary between individual runs due to ageing and damage of the thin-film. The b-coefficient is set to 0.035 for the Mach 10 flows and to 0 for the Mach 4.5 flow respectively. The uncertainty of the thin-film wall heat flux measurement is estimated to 15%.

The value of the b-coefficient between 0 and 0.035 has been adopted from previous experiments. The value of 0.035 is used for the Mach 10 experiments (see part 3.1.2) as it had been determined empirically in previous similar ISL experiments ([8] and [START_REF] Biegger | Wärmestrommessungen an Kugeln bei Machzahl 10[END_REF]). A comparison with wall heat fluxes computed alternatively setting b = 0 leads to a change in heat flux of less than 10%, as shown in Figure 8.6 in Appendix V. The weak dependence on the bcoefficient indicates that the increase of the wall surface temperature ∆T s is very low, only a few Kelvin. Thus, Equation 3.7 can be simplified: ∆T s,cor r (t ) = ∆T s . The small variation of surface temperature is explained by the low flow density and short flow duration at the Mach 10 conditions 3.1. This entails a low heat transfer coefficient and consequently a low wall heat flux as expressed in Equation 2.24. This argumentation also supports the use of a value of 0.035 for the b-coefficient for the Mach 10 conditions.

On the other hand, for the Mach 4.5 conditions the b-coefficient plays a significant role in the wall heat flux calculation, as shown in Figure 8.7 in Appendix V. The wall heat fluxes computed with b = 0.035 are up to two times higher than those computed with b = 0. So the wall heat flux measurement in dense flow is highly dependent on the b-factor value. The comparison with the CFD simulations presented in Section 7.1.2 is used for the determination of the correct b-coefficient of zero. In other words, the sensor temperature increase is still low enough such that the sensor material properties stay nearly constant, supporting the correction term in the Equation 3.7 at a value of zero.

It is to note that the thin film sensors are fragile and do not withstand impacts of the membrane debris erratically expelled from the shock tube. For the Mach 4.5 tests, a minimum of one sensor located in the stagnation point region was destroyed per test because of an impact. The missing data manifest themselves by a discontinuity in the wall heat flux profiles presented in Appendix V.

Flow conditions

Six flow conditions using nitrogen as test gas are employed for the tests. They are summarized in Table 3.1 and cover two Mach numbers, Mach 10 and Mach 4.5, respectively (see Appendix V for the altitude computation). The wall heat flux profiles obtained for the Mach 10 flows are used as validation cases, since results for the same conditions and the same model were already measured previously at ISL ( [8] & [START_REF] Biegger | Wärmestrommessungen an Kugeln bei Machzahl 10[END_REF]). The Mach 4.5 flows correspond to the flow conditions used for the ablation study, presented later in the Chapter 4. According to the Reynolds number, the boundary layer along the sphere forebody is expected to be laminar for the Mach 10 flows and partially turbulent for the Mach 4.5 flows. Table 3.1: Flow conditions for the wall heat flux measurements

Condition Number M ∞ h Re D u ∞ T ∞ P ∞ T 0 h 0 n • of runs (-) (km) (-) (m • s -1 ) (K ) (P a) (K ) (J • kg -1 )

Results

The raw data are available in Appendix V. The averaged wall heat flux profiles for each flow condition are presented in Figure 3.2. The errorbars show the statistical uncertainties and correspond to a confidence interval of 95% as described in Appendix V. The high uncertainties computed for locations near the stagnation point for the Mach 4.5 conditions are not only the results of run-to-run variations but also the results of missing data due to broken sensors caused by membrane impacts. The wall heat flux profiles measured at Mach 10 are symmetrical centred on the stagnation point. Except near the stagnation point, the differences between individual measurement points are within the uncertainty bonds. This indicates a flow repeatability between individual runs. In the range [-12.5 • ; 12.5 • ], the differences between the runs are much higher and cannot be explained by the sensor accuracy only. Deviations for the stagnation wall heat flux of up to 20% are found and attributed to run-to-run variations. Also, the position of the wall heat flux maximum is not measured at the stagnation point but rather at a location ±12.5 • off-axis from the model symmetry axis. The exact position, however, cannot be precisely assessed due to the measurement uncertainty and also the rather large sensor spacing.

The second group of wall heat flux profiles were obtained for the Mach 4.5 conditions. The asymmetries of the wall heat flux profiles are more pronounced upstream of the 40 • position. As for the Mach 10 conditions, the maximums of the wall heat flux are not located on the stagnation point but rather near 25 • from the symmetry axis. Figure 3.3 summarizes the position of this maximum for all flow conditions as a function of the Reynolds number. The measurement accuracy is limited by the sensor spacing, which is 12.5 • for this model. Figure 3.4a shows the dimensionless wall heat flux profiles qw / qw,0 . The profiles obtained for the conditions at Mach 4.5 and altitudes 3km and 6km clearly differ from others. Indeed, the ratio of the maximum wall heat flux to the stagnation wall heat flux is much higher for these conditions, reaching a value of 1.75 at a position around 25 • . This high value could be explained if boundary layer becomes turbulent upstream of the position of the maximum wall heat flux. As an example, the difference in heat transfer coefficient between laminar and turbulent boundary layers can be found for a flat plate in Figure 2.2. For a Reynolds number between 10 6 -10 7 , the ratio of the maximal turbulent wall heat flux to the laminar one is found between 2 and 4 in [START_REF] Luneau | Contribution à l'étude des phénomènes aérothermiques liés au vol hypersonique décéléré[END_REF] & [START_REF] Michel | Frottement et flux de chaleur turbulents en écoulement bi et tridimensionels[END_REF]. The possibility of a laminar-turbulent transition is further supported by Figure 3.4b which shows the local Stanton number as function of the local Reynolds number. For the Mach 10 flows and for all sensors positions, the local Reynolds number never exceeds a value of 10 5 . On the contrary, for the Mach 4.5 flows the local Reynolds number is always above 10 5 for all sensors except the one at the stagnation point. The shock tunnel is not considered as quiet wind tunnel, so a transition is likely under the Mach 4.5 conditions. We should mention that the NACA technical note [START_REF] Stine | Theoretical and Experimental Investigation of Aerodynamic-Heating and Isothermal Heat-Transfer Parameters on a Hemispherical Nose With Laminar Boundary Layer at Supersonic Mach Numbers[END_REF] reports a similar shift in the position of the wall heat flux maximum at Mach 1.97 and Re D ≈ 10 6 . The authors indicate that the boundary layer was laminar, based on optical Schlieren visualisations.

By using CFD simulations, it is shown later in Chapter 7 that although the freestream turbulence intensity of the shock tunnel is not known, the laminar-turbulent transition in the boundary layer is most likely responsible of the characteristic wall heat flux profiles of the Mach 4.5 flows at high flow density.

Discussions about the Camel effect

Wall heat flux measurements on a sphere at Mach 10 were already performed earlier at ISL ( [8] & [START_REF] Biegger | Wärmestrommessungen an Kugeln bei Machzahl 10[END_REF]). The comparison between the previous results and those obtained during the present work for the same flow condition are presented in Figure 3.5. The previous results obtained for the same geometry in the same facility [8] also show a similar shift of the maximal wall heat flux positions towards both sides of the stagnation point. It is referred to as the "Camel effect" since the wall heat flux profiles obtained show two maximums making the curve appear like the two humps of a camel. No explanation for this effect is given in the literature. Adeli et al. [8] pointed out the inability of the CFD simulation to reproduce this phenomenon. The Reynolds number of all tests suggests that the Camel effect is not related to the turbulence. Moreover, in Figure 3.6, the comparison of all dimensionless wall heat fluxes with a theoretical Newtonian distribution does not highlight a specific condition in terms of flow density or model diameter for the onset of the Camel effect.

The Figure 3.5a highlights that the reproducibility of a low density shock tunnel flow is somewhat problematic. For example, the wall heat flux for Mach 10 and h = 70km of the present work is much closer to the one obtained for Mach 10 and h = 60km in the previous work. In spite of this deviation, the local Stanton numbers of the old and new experiments are all much closer together for each flow density as shown in Figure 3.5b. Dimensionless wall heat flux q w /q w, 0 ( )

M10

Flow conditions

The Mach number 4.5 flow is selected for the ablation study. Even if the flow enthalpy for this Mach number is not the maximum achievable with the ISL shock tunnel, the flow density is high enough to maximize the energy flux while keeping the steady flow duration sufficiently long. The flow is generated by expanding the test gas in the Mach 4.5 contoured nozzle. Nitrogen is used as test gas in order to avoid an oxidation reaction during the ablation process. It should be mentioned that the dissociation of molecular nitrogen is negligible since the maximal stagnation temperature is below 1600K. Therefore, the model ablation is expected to be solely a combination of material melting and erosion.

Four flow conditions, summarized in Table 4.1 are used for the tests. The freestream conditions are computed from the following experimentally measured parameters: initial driver and test gas pressures, the shock propagation velocity in the shock tube and the stagnation pressure after shock reflection at the nozzle throat. The steady flow duration for these conditions is about two milliseconds, during which the stagnation pressure varies by 15%. For all flow conditions the Knudsen number is much lower than 0.03, so that the flow can be considered as a continuum flow without surface slip, as described in Section 2.1.2. Due to the combination of the noisy shock tunnel flow and the freestream Reynolds number [START_REF] Oertel | Strömungsmechanik: Methoden und phänomene[END_REF], the boundary layer is expected to become turbulent in the vicinity of the stagnation point region.

Table 4.1: Flow conditions for the ablation experiments Condition

M ∞ h Re/x u ∞ T ∞ P ∞ T 0 h 0 P 0 n • (-) (km) (m -1 ) (m • s -1 ) (K ) (P a) (K ) (J • kg -1 ) (W • m -2 )

Low temperature ablator

Due to the short steady flow duration of only two milliseconds, metals commonly used for aerospace applications such as steel or titanium cannot be used in this ablation study, as they do not have enough time to melt. For this reason, a low-temperature ablation material has to be used. Materials used in previous studies are not well-suited for liquid ablation experiments in shock tunnels. Naphthalene [START_REF] Combs | Heat-shield ablation visualized using naphthalene planar laser-induced fluorescence[END_REF] and camphor [START_REF] Stock | Hypersonic low temperature ablation an experimental study of cross-hatched surface patterns[END_REF] are ablated by sublimation, while models made of water ice [START_REF] Silton | Ablation onset in unsteady hypersonic flow about nose tip with cavity[END_REF] or CO2 ice [START_REF] Callaway | Photogrammetric measurement of recession rates of low temperature ablators in supersonic flow[END_REF] melt and also evaporate during the test preparation.

Gallium has been chosen as substitute ablation material for the following reasons: it is a non-toxic metal with a low melting temperature of 302.9K and a high boiling point of 2477K, so a purely liquid ablation can be expected during the test. Gallium properties at 10 5 Pa are given in Table 4.2. No significant difference to these values is observed at lower pressures. As also the case for other metals, a thin oxide layer covers the surface of the gallium model. Despite a high melting temperature of 2070K for the gallium oxide β-Ga 2 O 3 [START_REF] Víllora | Electrical conductivity and carrier concentration control in β-ga2o3 by si doping[END_REF], the experiments described later demonstrate that the oxide layer has a minor influence on the ablation process and is quickly stripped off by the flow.

With the exception of its low melting point, the thermal properties of gallium are of the same order of magnitude as those of other commonly used metals. The gallium ablation can thus be considered similar to theirs. Table 4.2: Gallium properties for the solid phase (from [2] and [START_REF] Ayrinhac | Compilation of thermodynamic and elastic data of liquid Ga[END_REF])

ρ G a C p G a λ G a ∆H f ,G a T f ,G a T s,G a σ el σ t s (kg • m -3 ) (J • kg -1 • K -1 ) (W • m -1 • K -1 ) (J • kg -1 ) (K ) (K ) (M P a) (M P a)
5904 370 41 80200 302.9 2477 8-25 15-40

Model geometry

Three different model geometries are used for the ablation experiments: a 8mm-hemispherecone-cylinder, a 10mm-hemisphere-cone-cylinder and a 14mm-hemisphere-cylinder. The geometries with their dimensions are shown in First, gallium is melted and heated up to 333K. At the same time, 3D-printed molds are assembled on the model sting and preheated to 333K for 15min. Then the liquid gallium is cast into the molds and left to cool down at ambient temperature during 15min. Thereafter, the model is cooled down to near 273K until the gallium solidifies. Finally, the molds are removed and the model is polished before being used in the experiment. This procedure ensures that the model is free of blowholes and is mechanically smooth. Pictures of the model production tools are presented in Figure 4.2. It is also to note that the gallium ablation experiment is difficult to conduct in summer, because the laboratory temperature may reach the gallium melting temperature and so the model begins to melt before the test.

A summary of the performed runs with the corresponding Reynolds number is given in Table 4.

Optical setup

A direct visualization system, as schematically shown in Figure 4.3, is used for the determination of the model position and for the tracking of the surface recession. This optical setup produces model images without any interference by the light refraction induced by the bow shock wave. A view of the optical setup is presented in Figure 4.4. The light source consists of a 6kJ flash lamp illuminating a diffuser, which provides nearly constant background illumination for as much as 10ms.

Two high speed cameras Photron SA-X record the test at 12500 frames per second with an exposure time of 0.293µs, ensuring negligible motion blur. Using a beam splitter and two different camera lenses, different magnification levels for each camera are achieved. The first camera records a low magnified image of the model with a first camera lens of 300mm focal length combined with a 1.4 extender placed before the beamsplitter. The second camera records a highly magnified image through an additional camera lens of 105mm focal length. With this setup, both cameras share the same line of sight with the approximate magnifications of 0.02mm/px for the first and 0.005mm/px for the second. The exact magnifications are measured before each test, as they are further required for the image post-processing. The low magnified camera records an image of the entire model and part of the sting, while the high magnified camera view only a part of the hemisphere corresponding to a region delimited by an angle of approximate ±40 • from the symmetry axis. A picture of the cameras and the beam splitter is presented in Figure 4.5a and an example of images obtained with the two magnification levels is shown in Figure 4.5b. • Filtering, grey-scale dilatation and threshold: Prior to the surface contour extraction, each frame is normalized by applying successively a median blur, a gray scale dilatation and a threshold. The threshold level is set so that the light reflections from the test section windows do no affect the following surface extraction.

• Ablation profile extraction: The contour of the model is extracted for each video using the algorithm [START_REF] Suzuki | Topological structural analysis of digitized binary images by border following[END_REF] implemented in the function findContours() included in the OpenCV library [START_REF]Open Source Computer Vision[END_REF]. Finally the ablation profiles are computed by determining the minimum distance between the reference image contour and the actual frame contour for each point of the contour.

Experimental results

Qualitative description

Ablation is observed for all flow conditions and for all models, primarily in the vicinity of the stagnation point region on the hemispherical part. The generated liquid layer is then sheared off along the model surface by the flow. No ablation by erosion was observable in any video, confirming that the main ablation process in this experiment is the ablation by melting.

The model images presented in Figure 4.7 do not show the ablated region but rather the region affected by the liquid layer, which re-solidified once the steady flow had passed. They do, however, provide a first indicator of the ablation intensity. As expected, the ablation is more important for the conditions with the highest flow enthalpy and flow power, namely flow conditions h = 1km and h = 3km. For condition h = 13km, on the contrary, the region contaminated by the re-solidified ablation layer is restricted to the model hemisphere only, thus indicating only weak ablation.

For the measurements at high enthalpy, it is to note that debris of the membrane, which initially separates the test section from the shock tube, sometimes impact and damage the model surface. In this case, the test result is rejected.

The Figure 4.8 shows the overlays of the last and first frame of the video during the steady flow time. They are done manually and the white regions in the front of the model correspond to the ablated regions. First observations show the ablation occurs only on the hemisphere part and not on the cone or on the cylinder. Moreover the ablation seems more intense on the side of the stagnation point.

The conclusion of the qualitative results is that the generation and the study of ablation is possible within the ISL shock tunnel. However the flow condition is restricted in terms of flow density. If the flow density is too high, the risk of membrane impact and thus damage to the model increases. If the flow density is too low, the aerothermal heating is also too low and the ablation intensity is too weak to be quantified.

In the following part, only quantitative results obtained for flow conditions h = 1km and h = 3km are presented, because the surface recessions for the conditions h = 6km and h = 13km are too small compared to the measurement uncertainties.

Effect of the flow condition on the heating phase duration

The heating phase duration is defined as the period between the steady flow start and the ablation start. During this phase the model wall temperature increases from the initial temperature to the melting temperature of the material. The steady flow time is experimentally defined when the pressure signal at the nozzle end reaches a nearly constant value. The ablation start time corresponds to the first time when a liquid layer is visible on the test videos.

The heating phase durations extracted from the test videos are presented in Table 4.4. The measurement uncertainty of about 0.2ms is mainly caused by the uncertainty on the nozzle steady flow starting point. For the high enthalpy conditions h = 1km and h = 3km, the ablation starts earlier compared to the lower enthalpy conditions due to the increase of the wall heat transfer.

Ablation start time and heat transfer in a semi-infinite medium

Due to the short test duration, the aerothermal heating affects only the surface of the gallium model. With the assumption of a thin thermal layer affected by the heating, the temperature increase near the surface can be described by the unsteady heat transfer theory for a semi-infinite medium, which is the case in the present experiment. Then, the temperature profile in the model is provided by the following heat transfer model with a fixed wall heat flux as boundary condition ( [START_REF] Knikker | Transferts Thermiques[END_REF] & [START_REF] Landau | Heat conduction in a melting solid[END_REF]):

ξ = x 2 at (4.1) T (x, t ) = T i + 2 qw,s λ at e -ξ 2 π + ξer f (ξ) -ξ (4.2)
First the thin thermal layer assumption is tested with the Equation 4.2. During a period of 0.5ms and by assuming a fixed wall heat flux of 9 • 10 6 W • m -2 , the dept of penetration is 0.24mm for a temperature change of 1K and 0.35mm for a change of 0.1K . Compared to the model diameter, the assumption of a semi-infinite medium is realistic for the determination of the ablation start time.

The theoretical ablation start time is given by: T (x = 0, t abl at i on st ar t ) = T f (4.3)

t abl at i on st ar t = π 4 ρC pλ T f -T i qw,s 2 (4.4)
In Figure 4.9a, the measured ablation start times are compared with the theoretical values.

The uncertainty of the analytical start time is represented by the green area. The parameters used with their uncertainties are available in Appendix V. Because the temperature difference between initial and melting temperature is only around 10K , a small initial temperature uncertainty leads to a significant uncertainty on the ablation start time. However, analytical and experimental start times are of the same order of magnitude and the evolution of the ablation start time in function of the wall heat flux follows the same trend. All measured ablation start times are higher than the theoretical ones. Two effects likely explain this difference: Firstly, the initial temperature of the model is lower than expected. During the test preparation, the test section is evacuated down to 5P a and the surface temperature decreases in spite of the high thermal inertia of the sting and the test section vessel.

Secondly, the real stagnation wall heat flux is lower than calculated using the correlations above. In Figure 4.9b, the experimental stagnation wall heat flux, computed from the ablation start time using Equation 4.4, is compared to the analytical stagnation wall heat flux obtained from the correlations presented in Section 2.2.1. It is the average of the values obtained with the Sutton-Graves [START_REF] Duffa | Ablative thermal protection systems modeling[END_REF], Oertel, Seiler, Smeets , Anderson and van Driest [START_REF] Van Driest | The problem of aerodynamic heating[END_REF] correlations. A difference of up to 50% is found between the experimental and analytical wall heat flux. It should be emphasized that the thin oxidation layer on the model surface can also delay the ablation start due to the higher melting temperature of the gallium oxide.

Ablation profiles as a function of time

The present method allows the tracking of the ablation profiles as a function of time for the high enthalpy flow cases. An example of ablation profiles as a function of time is given in Figure 4.10. At the early phase of ablation the measured profile is somewhat asymmetric. From 1.5ms onwards, the profile asymmetry is less pronounced. Even though some local asymmetries remain, the observations show that the ablation globally preserved the model symmetry. The ablation profiles for angle above 60 • are not correctly resolved, because of the lower precision of the low magnification video (the high magnification one is limited to a region less with than 40 • ) and because the ablation layer thickness at the edge of the hemisphere represents up to one pixel on the low magnification video.

From the ablation start time onwards the profiles are similar in their overall shape, featuring two local maximums per half profile (profiles with positive or negative angles): the first located at the stagnation point and the second near 20 • -40 • . These maximums are even more pronounced at 2ms. Similar ablation profiles were observed on water ice models in supersonic wind tunnel tests [START_REF] Silton | Ablation onset in unsteady hypersonic flow about nose tip with cavity[END_REF] as well as on aluminium spheres in a hyperballistic free flight tunnel [START_REF] Luneau | Contribution à l'étude des phénomènes aérothermiques liés au vol hypersonique décéléré[END_REF]. In these experiments, the location for maximum ablation was related to the position of the laminar-turbulent transition of the boundary layer.

Effect of the hemisphere diameter on the ablation profiles

The ablation profiles obtained at 2ms at flow condition h = 3km for the three different hemispheres are presented in Figure 4.11. Despite some scattering in the data and a surface displacement of less than 0.1mm, the characteristic ablated shape with the two ablation maximums is clearly observable for each tested geometry.

The first ablation maximum is located at the stagnation point, where the laminar wall heat flux is maximal and the wall shear stress tends to zero. The second ablation maximum corresponds to the position of maximal ablation and is always observed near 30 • from the symmetry axis. The higher ablation rate results from an increase in the melting rate due to an increasing wall heat flux. The Reynolds number regime and the noisy shock tunnel flow suggest that the boundary layer may become turbulent somewhere along the hemisphere surface. This would lead to an increase of the turbulent wall heat flux downstream of the laminar-turbulent transition location. According to [START_REF] Luneau | Contribution à l'étude des phénomènes aérothermiques liés au vol hypersonique décéléré[END_REF] & [START_REF] Michel | Frottement et flux de chaleur turbulents en écoulement bi et tridimensionels[END_REF], the maximum turbulent wall heat flux should be located at 36 • . What can be observed in Figure 4.12a, however, is a shift of the ablation minimum located between 10 • and 20 • . It is shifted from 15 • -20 • for the 8mm-hemisphere to 10 • for the 14mm-hemiphere. It is in agreement with the expected direction of the displacement of the laminar-turbulent transition and the location of the ablation minimum could be defined as the location of the laminar-turbulence transition.

Ablation profiles as a function of the flow conditions

The dependency of the ablation profiles on the flow condition are presented exemplarily for the 10mm-hemisphere-cone model. The test videos show that the gallium model at the lowest enthalpy condition h = 13km is never ablated within the 2ms. For the condition h = 6km, the model starts to be ablated by the flow but the surface displacement is too small compared to the measurement uncertainty and thus the ablation profile cannot be quantitatively determined. On the contrary, ablation profiles can be successfully measured for the flow conditions h = 1km and h = 3km. Due to frequent impacts of shock tunnel membrane debris on the model, only two valid ablation measurements are available for condition h = 1km. Thus for better comparability, the surface displacement d is normalized by the surface displacement at the stagnation point d 0 , yielding d * :

d * = d d 0 (4.5)
The ablation profiles are presented in Figure 4.12b. They clearly show the same trend for both conditions h = 1km and h = 3km. Indeed they are similar in shape with the two ablation maximums located at nearly the same position on the hemisphere: one at the stagnation point and one near 30 • , thus suggesting no change in the overall physics of the ablation process. The scattering in the data prevents the determination of the precise location in the ablation minimum and maximum, making it impossible to determine whether the location of the laminar-turbulent is shifted as a function of the Reynolds number.

Discussions

The flow conditions for the shock tunnel experiment correspond to a hypersonic flight at low altitude, so the flow density is high, particularly if it is compared to the flow generated in a plasma torch facility or flow encountered during high-altitude reentry of spacecraft. For this reason, the measured ablation profiles are characteristic ablation profiles for a dense flow. They differ from the measured ablation profiles in rarefied or high altitude flows, because the position of the ablation maximum is not located at the stagnation point. The wall heat flux profiles presented in the Chapter 3 and the ablation profiles show the same evolution and the same maximum locations.

In addition, the measured ablation profiles in the present work are very similar to the ablation measurements obtained for similar Reynolds numbers and available in the literature. By measuring the ablation profile on an aluminum sphere in free flight, Luneau [59] found a similar ablation profile. In a nitrogen atmosphere, he demonstrated that the ablation maximum without combustion (or oxidation reaction) is located near 30 • ± 5 • . For Luneau, the laminar-turbulent transition is located in this region and increases the local wall heat flux, thus also increasing the local ablation intensity. A second similar ablation measurement was presented by Silon & Goldstein [START_REF] Silton | Ablation onset in unsteady hypersonic flow about nose tip with cavity[END_REF] on a water ice model placed in a Mach 5 air flow generated by a supersonic wind tunnel. The location of the ablation maximum was found near 45 • and the authors also cited the laminar-turbulent transition to be the cause. Pictures of the aluminium sphere and of the water ice model are shown in Figure 4.13.

The Tables 4.5 and 4.6 compare the freestream conditions and the critical Reynolds number Re c = U ∞ • R • θ abl a,max /ν ∞ = 0.5 • Re D • θ for the ablation experiments. The critical Reynolds number is computed with the assumption that the location of laminar-turbulent transition coincides with the location of the ablation maximum. All considered freestream Reynolds numbers are of the same order of magnitude. The location of maximal ablation does not appear to be correlated directly with the local Reynolds numbers. The critical Reynolds numbers cover two orders of magnitude from 10 4 to 10 6 . No experimental data for the transition location over an hemisphere in hypersonic flow have been found in the literature. For an incompressible flow, however, the critical Reynolds number for a Blasius boundary layer is 5 • 10 5 and data for hypersonic transition are available for cones [START_REF] Anderson | Hypersonic and high-temperature gas dynamics[END_REF]. For a sharp cone at Mach 4.5, the critical Reynolds is measured between 2 • 10 6 and 2 • 10 7 . Assuming that the ablation maximum is caused by the transition to the turbulence, it shows that the transition over a hemisphere occurs much earlier than for a sharp cone despite the strong negative pressure gradient along the sphere. Eventually the three experiments lie in a grey zone where the critical Reynolds numbers appears to be low compared to the expected one. However, for the shock tunnel and the wind tunnel facilities, the elevated freestream turbulence could explain the quick transition after the stagnation point despite the low Reynolds number.

(a) Ablated aluminium sphere after a free flight test [START_REF] Luneau | Contribution à l'étude des phénomènes aérothermiques liés au vol hypersonique décéléré[END_REF] (b) Water ice model ablated in a supersonic wind tunnel [START_REF] Silton | Ablation onset in unsteady hypersonic flow about nose tip with cavity[END_REF]

Part III

Numerical investigation

This part presents the tool for the prediction of the wall heat fluxes and the ablation along the forebody parts of a flight vehicle. The validation of the tool is described as well as the simulations of the experiments presented in the Part II. The analysis of the numerical results further supports the conclusions on the turbulence drawn based on the experimental results.

The solver for heat transfer and ablation has been developed within the OpenFOAM framework. It is compatible with OpenFOAM v2012 and is validated on Ubuntu 18.04.6 LTS with the GCC compiler. The OpenFOAM framework was chosen because it is a free, open source CFD code under GPL license. In addition numerous modules, such as ones for mesh manipulation or for the thermodynamic modelling, are already available. Finally, the OpenFOAM environment has been used in recent research projects for the modelling of multiphase porous reactive materials, as for example for simulations of heat shield ablation (PATO project [7]) or for the modelling of electric discharges in hypersonic flows [START_REF] Nekris | Numerical Modelling of Electric Discharges Generated in High-Speed Flows[END_REF].

From the programmer's point of view, OpenFOAM is a C++ programmed toolbox dedicated for solving partial differential equations with the finite volume method. It is an objectoriented code, so the modules are programmed inside classes. Solvers and utilities can easily be created by extending the existing modules.

The solver for the aerothermal heating and ablation is named ablationFOAM. It is based on merging the pressure based solver sonicFOAM [START_REF] Opencfd | OpenFOAM API guide v2012: sonicFOAM.C[END_REF] with the heat transfer solver chtMul-tiRegionFOAM [START_REF] Opencfd | OpenFOAM user guide v2012: chtMultiRegionFOAM[END_REF]. sonicFOAM is a transient solver for compressible, trans-/supersonic and turbulent flows. It has been preferred to the density based solver rhoCentralFOAM [START_REF] Opencfd | OpenFOAM user guide v2012: rhoCentralFOAM[END_REF] because ISL-internal studies had indicated that the supersonic boundary layer is more accurately simulated with the solver sonicFOAM. chtMultiRegionFOAM [START_REF] Opencfd | OpenFOAM user guide v2012: chtMultiRegionFOAM[END_REF] is a solver dedicated to transient simulations of heat transfer problems. It supports simulations with multiple regions, which can be solid or fluid domains.

It should be noted that the chemical reactions in air or in nitrogen and at the solid/gas interface are not yet implemented in the ablationFOAM solver, since no experiment in the present work presented includes such effects. However modules for combustion and chemical reactions are already available in OpenFOAM and Nekris [START_REF] Nekris | Numerical Modelling of Electric Discharges Generated in High-Speed Flows[END_REF] and Casseau [START_REF] Casseau | An open-source CFD solver for planetary entry[END_REF] proposed implementations of thermochemical models for high speed flows in OpenFOAM.

Chapter 5

Mathematical model

Solver structure

The numerical solver ablationFOAM is a multi-region, segregated, pressure-based solver for compressible and turbulent flows. It also simulates the internal heat conduction in solids. Some features are listed and explained below:

• Transient solver: the solver is a transient solver, so the time derivative terms of the conservation equations are included into the matrix system to solve. The size of a time step is described via the Courant-Friedrich-Lewy (CFL) number. For a transient simulation, the maximum CFL number practically possible is 1.5-2. When steady flows are simulated, the solver is used as a pseudo-transient solver with an implicit relaxation of the matrix system. It stabilizes the simulation and a CFL number up to 100 is possible for coarse grids.

• Multi-region solver: multiple regions can be defined with an associated mesh for each one. The solution of each equation system is independent between the different regions and the data communication occurs at the domain interface only. The number of allowed domains is virtually unlimited. In the present work, two domains are defined: one representing the gas domain and the second the solid domain.

• Segregated solver: the set of conservation equations is solved iteratively. One matrix system is created and solved for each conservation equation. During the solution of a matrix system, all fields excepted the one being currently solved are set constant. This solution method minimizes the memory requirements but the convergence rate is slower compared to a coupled solver.

• Pressure-based solver: the pressure is used as a main variable and is solved with a pressure equation. The density is directly computed from the pressure and temperature with an equation of state. The solver structure follows the PIMPLE algorithm.

• Compressible flow: the solver is suited for non-constant density flows and works for the entire Mach number spectrum. In the limit of incompressible flow, a pressure variation leads to a velocity change, while in highly compressible flow a pressure variation leads to a density change.

• Turbulent flow: Turbulent flow can be simulated with the ablationFOAM solver. Taking advantage of the object-oriented programming style provided by the C++ language, the turbulence models available in OpenFOAM can be used in this solver.

• Internal heat conduction in solid: the energy conservation equation is the only conservation equation solved for the solid domain. In a solid, this equation is reduced to an equation for the heat conduction.

• Mesh motion capability: the mesh of the solid domain moves and adapts in order to follow the melting front during the ablation process.

The post-processing is carried out with the OpenFOAM post-processing functions, the Paraview software [START_REF] Ahrens | Paraview: An end-user tool for large data visualization[END_REF], the Jupyter Lab Notebooks [6] and the Python libraries Numpy [START_REF] Harris | Array programming with NumPy[END_REF], Pandas [START_REF] Mckinney | Data Structures for Statistical Computing in Python[END_REF] and Matplotlib [START_REF] Hunter | Matplotlib: A 2d graphics environment[END_REF].

The solver procedure is shown in Figure 5.1. Each iteration starts with the detection of the solid domain cells affected by the ablation. The mesh is then moved and adapted according to the position of the melting front, as presented in Section 5.5. Next, the conservation equations for the fluid domain are solved with a PIMPLE algorithm. Finally the energy equation of the solid domain is solved. The equation solving process is repeated as many times as specified in the user-defined parameters. Then, a new time step begins and the the procedure is repeated until the simulation end time.

Model for the fluid phase

The fluid phase is treated as a Newtonian fluid and as a thermally perfect gas. The set of equations are:

• The mass conservation equation:

∂ρ ∂t + ∇ • ρ v = 0 (5.1)
• The momentum conservation equation:

∂ ρ v ∂t + ∇ • ρ v ⊗ v = ∇ • T + f v (5.2)
With:

-The stress tensor:

T = -p + λ (∇ • v) I + 2µD
-The deformation rate tensor:

D = 1 2 ∇ ⊗ v + (∇ ⊗ v)
T the bulk viscosity λ given by the Stokes hypothesis: λ = - 2 3 µ no volumetric force considered:

f v = 0
The momentum conservation equation can be written as:

∂ ρ v ∂t + ∇ • ρ v ⊗ v = -∇p + ∇ • τ ( 5.3) With:
-The viscous stress tensor: τ = T + p I

• The energy conservation equation:

∂e ∂t + ∇ • e + p v = -∇ • q + ∇ • τ • u (5.4)
In the present work, the viscous dissipation term ∇ • τ • u is neglected. The contribution of this term compared to the thermal conduction term -∇ • q is negligible, as shown in the following chapters with the comparisons of the numerical wall heat fluxes with the experimental data.

• The equation of state used for the density update:

ρ = p RT (5.5)
This set of equations is solved by using the PIMPLE algorithm. The temporal discretization in this work uses the 1st-order, implicit Euler scheme. Unless explicitly mentioned as different in the following parts, the default spatial discretizations used are the following:

• for the gradient operator: the cellLimited Gauss linear 0.5 scheme.

• for the divergence operator: the Minmod divergence scheme.

• for the laplacian operator: the Gauss linear limited 1 scheme.

• for the interpolation: the linear interpolation method.

The pressure equation is solved via a multigrid solver GAMG, while momentum and energy equation are solved with a preconditioned bi-conjugate gradient PBiCGStab (Newton-Krylov solver). In the steady state simulations, implicit under-relaxation is applied.

The fluid phase uses a single-species thermodynamic model. The temperature of the nitrogen flows in the shock tunnel experiments is much lower than 4000K, so no gas dissociation is expected. Moreover the gas flow is assumed to be in equilibrium. This assumption is confirmed by the comparisons of experimental and numerical data, except for the shock tunnel tests at the Mach 10 flow condition with the lowest flow density equivalent to an altitude of 70km. The vibrational excitation of the gas molecules is taken into account with a temperature-dependent specific heat capacity C p:

h = C p(T ) • T (5.6)
C p is evaluated with a forth order polynomial function. The coefficients are computed from those proposed in [START_REF] Mcbride | Coefficients for calculating thermodynamic and transport properties of individual species[END_REF]:

C p(T ) = R a 1 + a 2 T + a 3 T 2 + a 4 T 3 + a 5 T 4 (5.7)
The transport properties, the viscosity and the thermal conductivity, are modelled with eight degree polynomials constructed from the coefficients proposed in [START_REF] Mcbride | Coefficients for calculating thermodynamic and transport properties of individual species[END_REF]:

µ = 8 0 a i T i (5.8) λ = 8 0 b i T i
(5.9)

Model for the solid phase

The solid domain represents the structure of a hypersonic vehicle. It includes models for the thermal conduction, the phase change and the shape modification. Only the energy conservation equation is solved in the solid.

Energy conservation equation

The energy conservation equation uses the enthalpy as main variable. It includes terms for the thermal conduction, thermal convection and heat sources, it is written as:

∂ρh ∂t + ∇ • ρh ( v -v CV ) = -∇ • q + q (5.10)
The heat conduction term -∇ • q is treated explicitly based on the temperature q = -λ∇T . A peculiarity of this solver is that the convection term cannot be removed as it is for the classical heat equation in solid. However the mass convection is zero in the solid ( v = 0), thus the energy transport because of the mass transport is also zero. Also the mesh of the solid domain moves and is deformed in order to follow the melting front. This motion implies that each control volume, represented by the mesh cells, has is own displacement speed v CV , which creates an artificial energy flux passing through the cell faces ∇ • ρh v CV . Finally, the energy equation to solve is:

∂ρh ∂t -∇ • ρh v CV = λ∇ 2 T + q (5.11)
The equation is solved with a Newton-Krylov PBiCG solver. The temporal dicretization uses the 1st-order implicit Euler method, the same as for the fluid region. The spatial discretizations used are the following:

• for the gradient operator: the cellLimited Gauss linear 0.5 scheme.

• for the divergence operator: the linear scheme.

• for the laplacian operator: the Gauss linear limited 1 scheme.

• for the interpolation: the linear interpolation method.

Phase change model

The solid is considered as a single species and does not sustain any chemical reaction. The melting process implies that the thermodynamic state of the species does not only depend exclusively on the temperature T but also depends on the actual phase. Indeed the phase change occurs at a fixed temperature, as shown in Figure 5.2, thus the melting process cannot be described only by the temperature. The phase change is modelled by introducing the phase fraction α field. It is bounded between zero and one. One means that the cell is fully solid, zero that the cell is completely melted and the values in between correspond to the intermediate states. Figure 5.3 shows the dependence of the phase fraction on the enthalpy. After solving the energy conservation equation, the phase change α is directly computed from the new enthalpy value through the relation:

α(h) = max mi n 1 + C pT f -h H f ; 1 ; 0 (5.12)
After the calculations of the enthalpy and the phase fraction, the temperature field of the solid domain has to be updated. The classical method implemented in OpenFOAM is based on a Newton-Raphson method for iteratively finding the new temperature according to the enthalpy and the heat capacity. This method is not suited for a melting single species. Indeed during the melting process, the temperature is not related to the enthalpy and the derivative ∂h/∂T is infinite. To overcome this limitation, the temperature is computed directly by the following equation: With the combination of Equations 5.12 and 5.14, the temperature remains constant during the melting process and is determined as following:

h(T, α) = C pT + (1 -α) H f (5.13) ⇒ T = h(T, α) -(1 -α) H f C p (5.
     α = 0 ⇒ T = h+H f C p α ∈ ]0; 1[ ⇒ T = T f α = 1 ⇒ T = h C p
(5.15)

Gas-Solid boundary

The multi-domain capability in the simulation requires a data communication method at the interface of two adjacent domains. For the ablation simulations, the wall temperature and the wall heat flux are shared by the solid and the fluid phases, as presented in Figure 5.4a.

φ f = α e 2 φ nbr + (1 -α e 2 )(φ P + ∇φ nbr ∆x)
(5. 16)

∇φ f = α e 2 φ nbr -φ P ∆x + (1 -α e 2 )( ∇φ nbr ) (5.17)
With:

• the temperature used as wall value: φ = T

• the temperature gradient used as gradient at wall:

∇φ nbr = qnbr λ P = λ nbr λ P ∇T nbr • the mixing coefficient α e 2
• the first cell height ∆x

If α e 2 → 0, the boundary condition behaves as the Neumann boundary condition:

φ f = φ P + ∇φ nbr ∆x (5.18) ∇φ f = ∇φ nbr (5.19)
If α e 2 → 1, the boundary condition behaves as the Dirichlet boundary condition:

φ f = φ nbr (5.20) ∇φ f = φ nbr -φ P ∆x (5.21)
The preceding equations require a mixing coefficient. Originally, the OpenFOAM boundary condition 'turbulentTemperatureCoupledBaffleMixed' [START_REF] Opencfd | OpenFOAM API guide v2012: turbulentTemperatureCoupledBaffleMixed-FvPatchScalarField.C[END_REF] uses a mixing coefficient based on the mesh geometry. It is defined as the ratio of the first cell heights at the interface. It does not reflect any effect of the thermal property of the domain and it forces the use of a similar first cell height for the two domains in order to maintain the simulation stability. This method is inefficient because the height of the first cell in the fluid domain must be very small to ensure a sufficiently small Courant number. In comparison, the height of the first cell in the solid domain can be much larger, without any consequence on the results.

For this reason, a mixing coefficient based on a thermal property of the domains is preferred. The thermal effusivity is then used for this purpose:

e = λρC p (5.22)
The thermal effusivity denotes the capability of a material to exchange energy with its surrounding environment. At the interface of two domains with different thermal diffusivities, the one with the higher value will dominate the contact temperature, close to its own internal temperature. Based on this observation, the new mixing coefficient is constructed as the ratio of the square of the thermal effusivity: For the simulations of ablation, the thermal effusivity of the gas flow is negligible compared the thermal effusivity of the metal. Thus, the solid domain sets the contact temperature at the interface, which is equal to the metal melting temperature after the start of ablation. In the same time, the wall heat flux is set by the temperature gradient at the wall in the gas domain.

Mesh motion

The shape change of the solid due to the ablation is simulated via the deformation of the mesh, which adapts itself to the melting front in the solid domain. The procedure for the mesh motion consists of two steps:

• the detection of the ablated cells at the surface of the solid domain.

• the determination of the surface displacement for the ablated cells.

• the propagation of the mesh surface motion inside the solid domain.

As shown in Figure 5.1, the detection of the ablated cells is the first operation within a time step. On the fluid-solid interface, the enthalpy is read. If it exceeds the value after melting defined as C pT f + H f , the interface is considered as completely melted and the new position for the melting front is determined. This position corresponds to the location where the enthalpy is equal to C pT f + H f . It is found by using the enthalpy value and the enthalpy gradient normal to the face according to:

∆d abl at i on = C pT f + ∆H f -h i nt er f ace
∇h cel l cent er (5.24)

Figure 5.5 illustrates this procedure. According to the surface displacement caused by the ablation, the cell face displacement is then converted into a cell face velocity: v abl at i on = ∆d abl at i on ∆t (5.25)

Once the mesh velocity at the boundary is known, this velocity is propagated into the mesh domain by solving the Laplacian equation:

∇ • γ∇ v cv = 0 (5.26)
The term γ represents the diffusivity factor, which is used for the control of the mesh motion near the wall. For the ablation simulation, it is set to 'uniform' (same mesh motion diffusion into the entire domain) or to 'inverseDistance' (the mesh motion diffusion decreases inversely with the distance to the wall. The prescribed mesh motion should satisfy the Space Conservation Law (SCL) [START_REF] Jasak | Dynamic mesh handling in OpenFOAM[END_REF]:

∂ ∂t V ∂V - S n • u cv ∂S = 0 (5.27)
As seen in Equation 5.11, the mesh motion v cv is included in the convection term. The mesh motion and the validation of the 5.27 condition are handled by the dynamic mesh library of OpenFOAM, which also corrects the flux φ of the energy conservation equation according to the mesh motion.

Solver verifications and test cases

The verification of the solver is carried out with test cases for which the analytical solutions are known:

• The propagation of discontinuities and the treatment of strong gradients are investigated with a Riemann problem: the 1D-shock tube problem. (Section 6.1)

• The shock wave modelling is validated with the post-shock values computed for super-/hypersonic flows around a 15 • -wedge. (Section 6.2)

• The heat transfer in the solid is validated with a collection of test cases. They are used for the validation of the heat conduction models and of the data exchange between two mesh domains. (Section 6.3)

• The mesh motion algorithm is checked via several mesh deformations. (Section 6.4)

• Finally the model for ablation is validated with Stephan problems. (Section 6.5)

Riemann problem

The Riemann problem is a classical numerical test case for the verification of the accuracy and the numerical stability of a solver. It provides an analytical solution to the Eulerequation for the propagation of a shock wave in an infinite medium. At the initial time, a 1D-domain is divided into two sub-domains: a high pressure part and and low pressure one. At the simulation start, the high pressure part expands through a propagating shock wave into the low pressure part.

The initial conditions are given in Table 6.1 and are taken from [START_REF] Saurel | Exact and approximate Riemann solvers for real gases[END_REF] (these conditions are also used in [START_REF] Henneaux | Towards a High-Fidelity Multiphase Solver with Application to Space Debris Aerothermal Ablation Modeling[END_REF]). In addition, the gas thermodynamic properties are :

R = 231.11J •kg -1 •K -1 and C p = 808.91J • kg -1 • K -1 .
The flow is considered inviscid. The results are taken at a time t = 8•10 -6 s. The simulations are performed with meshes of 100, 1000 and 2000 cells. For each mesh, simulations are performed with Courant numbers of 10 -1 , 10 -2 and 10 -3 . The Figures 6.1, 6.2 and 6.3 present the simulated profiles at t = 8 • 10 -6 s and the Figures 6.4, 6.5 and 6.6 a zoom on the shock wave located near the position x = 0.35m. The pressure profile shows an expected profile for a shock tunnel case with expansion waves in the high pressure part and a shock wave in the low pressure one. The computed pressure and velocity levels are in good agreement with the exact solution. The post-shock temperature level shows significant variation in function of the Courant number chosen. The correct temperature is found for Courant numbers below 10 -2 only. The zoomed-in profiles show that with a Courant number of 10 -1 the shock wave propagation is too slow compared to the exact solution, such that it leads to a wrong post-shock temperature. This observation is true independently of the mesh discretization. The mesh discretization has a minor impact on the solver accuracy. Even with only 100 cells, the solution matches with the exact solution. As expected with a coarse mesh, a slight numerical diffusion is observed at the profile corners.

In conclusion, the ablationFOAM solver is able to solve a Riemann problem. The result accuracy depends significantly of the time step. For transient simulations, a maximum Courant number should not be exceed 0.01 for a good match with the shock dynamics. The sharp capture of discontinuities is controlled by the spatial discretization. A cell density of 2500cells/m appears as a good compromise between computational time and result accuracy.

Super-/Hypersonic flows around a 15 • -wedge

The second verification test case is the super-/hypersonic wedge flow. It aims to validate the post-shock values and the oblique shock angle simulated by the solver. The analytical solutions are well known and can be found in many books on compressible fluid dynamics ( [START_REF] Anderson | Hypersonic and high-temperature gas dynamics[END_REF], [START_REF] Hirschel | Basics of Aerothermodynamics[END_REF]). The wedge angle is 15 • and the tested Mach number are 2, 3, 10 and 25. The flow conditions and properties are presented in Table 6.2 and correspond to a nitrogen flow. The perfect gas hypothesis is used and the flow is considered laminar. The viscosity is modelled with the Sutherland equation.

The computational meshes are shown in Figure 6.7. Three mesh sizes are used: a coarse (2479 cells), a medium (9916 cells) and a fine mesh (39249 cells). A view of the Mach fields for the Mach 3 flow is presented in Figure 6.7. The coarse mesh leads to a thick shock caused by the coarse spacial discretization. The successive mesh refinements reduce the shock thickness as expected. The simulated post-shock values and the wave angles are in very good agreement with the theory, as presented in Table 6.3. Even the results obtained with the coarse mesh show discrepancies below 2.5% despite the poor discretization of the shock. This demonstrates that the simulated post-shock values are nearly independent of the discretization across the shock. It is an important observation, since the maximum Courant number is reached in the refined cells just upstream of the shock. Computation time can be saved by reducing the shock refinement, without degrading the result It is important to note that convergence of the Mach 25 case is difficult to obtain. This is caused by the mesh, which is not well adapted for the small wave angle, and is also caused by the strong gradient across the shock, which leads to strong oscillations with the Min-Mod scheme. This is fixed by computing the divergence operator with the 1st-order upwind scheme. It produces enough numerical diffusion to stabilize the simulation. The Table 6.3 shows that the post-shock values are still correctly computed.

During these simulations, the final solution of one simulation is reused as initial state for the next simulation on a finer mesh. This method saves computational time by skipping the flow establishment in the convergence process, a period which is also often numerically unstable.

Finally, the 15 • -wedge case shows the ability of the solver to simulate highly compressible flows. Good results can be quickly obtained with a coarse mesh then more accurate results are achieved with finer meshes. Table 6.2: Flow conditions

M ∞ p ∞ T ∞ R C p (-) (P a) (K ) (J • kg -1 • K -1 ) (J • kg -1 • K -1 ) (
M ∞ Mesh size θ w ave error M 2 M 1 error p 2 p 1 error ρ 2 ρ 1 error T 2 T 1 error (-) (-) (%) (%) (%) (%) (%) 2 coarse < 1 2.5 < 1 < 1 < 1 2 fine & 3 shock refinements < 1 < 1 < 1 < 1 < 1 3 medium 1.1 1.7 2.5 < 1 < 1 3 fine < 1 < 1 < 1 < 1 < 1 10 fine & 1 shock refinements < 1 < 1 < 1 1.3 < 1 10 fine & 2 shock refinements < 1 < 1 < 1 < 1 < 1 25 fine & 2 shock refinements < 1 < 1 < 1 < 1 < 1

Heat transfer in solid

This section describes the verification of the heat transfer model in the solid domain. The test cases and the analytical solutions are taken from [START_REF] Knikker | Transferts Thermiques[END_REF] and [START_REF] Taler | Solving direct and inverse heat conduction problems[END_REF]. They are summarized in Table 6.6 with the corresponding heat transfer conditions. The solid properties are constant. The detailed simulation conditions are presented in Appendix V. Figure 6.9 compares the numerical temperature profiles in the solid domain with the analytical ones for the steady heat conduction problem in a single domain. All numerical results match with the analytical solutions even with a coarse mesh, as shown in Figure 6.9e. In conclusion, the model for the heat transfer in solid with multi-domain capability is validated and no discrepancy with the analytical solutions are observed. The resolution of the heat transfer equation is simple. The only limiting factor for the simulation speed is the size of the mesh cells, which naturally limits the time step size.

Mesh motion

This section is dedicated to the validation of the mesh motion algorithm. Five test cases are presented in Table 6.5. The full set of initial and boundary conditions is available in Appendix V. The tested domain is a 2D-beam with a uniform initial temperature. All boundary conditions are set to zero gradient, so that no energy transfer through the boundary is allowed. Thus the temperature must remain constant during the mesh motion.

A view of the deformed mesh is presented in Figure 6.11. No temperature variation is detected during the mesh motion. The flux correction, which modifies the convection term according to the mesh motion, is required.

Finally, the mesh motion test cases demonstrate the validity of the mesh motion model.

Stefan problem

The Stephan problem is a boundary value problem with a moving discontinuity inside the domain. It describes phase change phenomena such as melting. It was formulated first by Stefan [START_REF] Stefan | Über die Theorie der Eisbildung, insbesondere über die Eisbildung im Polarmeere[END_REF] in 1889 who investigated the ice formation on the polar sea. As shown in Table 6.6, two test cases are used taken from [START_REF] Taler | Solving direct and inverse heat conduction problems[END_REF]:

• the first case simulates the melting of a semi-infinite solid. The melting front moves inside the computational domain but the mesh remains static.

• the second case simulates the ablation by melting of a semi-infinite material. The mesh is dynamic and its boundary follows the melting front.

The initial and boundary conditions as well as the thermodynamic properties of the solid can be found in Appendix V. Melting 1D-beam Fixed temperature and phase change 2 Ablation 1D-beam Fixed temperature, fixed gradient and surface recession Figure 6.12 shows the results obtained for the melting case. The temperature in the solid never exceeds the melting temperature, which validates the temperature constrains of the melting expressed in Equation 5. [START_REF] Ayrinhac | Compilation of thermodynamic and elastic data of liquid Ga[END_REF].

The simulated melting front moves in agreement with the theoretical one given by the equation:

s(t ) = 2β λ C p • ρ t (6.1)
Where β is the solution of the equation:

βe β 2 erfβ = 1 π C p T w,2 -T f ∆H f (6.2)
These results validate the phase change model presented in Section 5.3.2. The results for the ablation test case are presented in Figure 6.13. The wall temperature in the solid increases until the melting temperature is reached. The theoretical ablation start time is given by:

t s = π 4 ρC pλ T f -T w,2 qw 2 (6.3)
It is equal to t = 1.96 for the considered ablation case. This time is represented by the red vertical line on the Figure 6.12b. The ablation start time obtained with the numerical simulation is in agreement with the analytical one.

A theoretical expression for the steady state ablation speed is given by:

v = d s(t ) d t = qw ρ C p T f -T w,1 + ∆H f (6.4)
The theoretical value is 0.01m •s -1 and is represented in Figure 6.13b by the black oblique line. After approximately 12s, the simulated ablation speed reaches its steady state value, which matches with the theoretical one. The ablation test case demonstrates the correct behaviour of the mesh motion module when operating together with the phase change module.

Simulations of the wall heat flux on a sphere

The wall heat flux represents the energy transfer between the fluid and the solid domain. The aims of this chapter is: first the validation of the ablation solver for the wall heat flux prediction in hypersonic flows and secondly the simulations of the ISL wall heat flux measurements presented in the Chapter 3.

Therefore three sets of simulations are performed and presented in the present chapter:

• The two first sets are related to the wall heat flux measurements in low density hypersonic flows. The flow is assumed laminar, so the simulations do not include any turbulence model. Two series of experiments are reproduced: the Hollis experiments [START_REF] Hollis | Experimental and computational aerothermodynamics of a Mars entry vehicle[END_REF] and the ISL shock tunnel experiments for low density flows presented in Chapter 3.

• The last set of simulations reproduces the ISL shock tunnel experiments for the high density flows presented in Chapter 3. The flow is treated as a turbulent flow and the effects of the turbulence models are investigated.

In general, the wall heat flux is dependent on the properties of both the solid and fluid domains. In the case of experiments performed with impulse facilities, the solid surface temperature stays almost constant during the considered period, as shown in Chapter 3. Therefore the wall heat flux is only controlled by the fluid domain. For this reason, the computational domain for the wall heat flux simulations is composed of the fluid domain only. It is also convenient because the thermal state of the model is often unknown or undocumented in the wall heat flux experiments.

Laminar flow

Hollis experiments [41]

The simulations presented in this section reproduce the experiments of Hollis described in the Appendix C of [START_REF] Hollis | Experimental and computational aerothermodynamics of a Mars entry vehicle[END_REF]. Wall heat fluxes are measured on hemispheres of 12.7mm-diameter (1/2 ), 19.05mm-diameter (3/4 ) and 25.4mm-diameter [START_REF] Denis | Standard Atmosphere. National Oceanic and Atmospheric Administration[END_REF]. The flow conditions are described in Table 7.1. The conditions n • 307 are generated with the hypersonic wind tunnel NASA Langley 31-inch Mach 10 air tunnel. The other flows studied are generated with the NASA HYPULSE expansion tube. For all conditions the freestream Reynolds number based on the model diameter does not exceed 3 • 10 5 , and no indication of a boundary layer transition is mentioned in the analysis of the experimental data. For these reasons, the flows are simulated as laminar flows.

Table 7.1: Flow conditions used for the simulations of the Hollis experiments [START_REF] Hollis | Experimental and computational aerothermodynamics of a Mars entry vehicle[END_REF] Condition A mesh independence study is performed with the flow conditions n • 307. An example of the 2D-axisymmetric mesh used for the hemisphere simulation is shown in Appendix V. The effects of several mesh parameters are investigated with three different meshes:

M ∞ Test gas u ∞ T ∞ P ∞ T 0 h 0 n • (-) (-) (m • s -1 ) (K ) (P a) (K ) (J • kg -1 )
• a coarse mesh: cell number along the wall: 100 cells, expansion ratio : r = 10%, y + = 1.

• a medium mesh: cell number along the wall: 500 cells, expansion ratio : r = 5%, y + = 0.1.

• a fine mesh: cell number along the wall: 100 cells, expansion ratio : r = 1%, y + = 0.01.

As shown in Figure 7.1, the mesh independence is observed even for the coarse mesh. Downstream of the position at 10 • , no result discrepancy is observed between the different mesh sizes. For each run, the experimental wall heat flux is measured with two models: a Macor model and a quartz model. A wall heat flux difference up to 20% is measured between the materials. The thermal state of the model is not known, for this reason the wall temperature is set constant to 293K in the simulation. Even with this simple boundary condition, the simulated wall heat flux profiles are in good agreement with the measurements for all positions located downstream of 10 • . They are within the measurement uncertainties mentioned above. Near the stagnation point region, the wall heat flux shows erratic behaviour and is unstable. Indeed upstream of 5 • -10 • , a carbuncle phenomenon is observed. The y + or the expansion ratio does not mitigate this effect. In contrast, a better parietal discretization along the stagnation point surface reduces but does not remove the amplitude of the carbuncle effect. The use of a 3D-mesh or of 1st-order schemes could improve the wall heat flux profile at the stagnation point. All the following simulations are performed with 2D-axisymmetric meshes. Despite the carbuncle effect, the variation of the wall heat flux between the stagnation point and the 10 • -position is small and well predictable. If the rest of the wall heat flux profile is correctly simulated, the use of a 3D-mesh is not considered mandatory.

The effect of the local mesh refinement on the shock is presented in Figures 7.2 and 7.3. During the simulation, three successive mesh refinements of the shock wave are applied, as shown in Appendix V. The cells to refine are selected based on their pressure gradient value. A pressure gradient threshold is set by the user for each level of the shock refinement.

The first mesh refinement is a general mesh refinement of the shock and the post-shock regions. The successive second and third mesh refinements refine only the bow shock. The boundary layer along the hemisphere wall is never affected by these refinements in order to avoid excessively small cells. The improved discretization of the shock enhances the quality of the wall heat flux profile near the stagnation point. The region affected by the carbuncle effect is reduced and is limited to a region up to 5 • from the symmetry axis. Downstream of this location, the shock discretization has no significant effect on the wall heat flux profile.

In conclusion, the mesh independence is demonstrated with the three simulations of the run n • 307. Successive shock refinements improve the quality of the wall heat flux profiles, especially near the stagnation point. To combine fast simulations obtained with a coarse mesh and accurate results obtained after several mesh refinements, the following procedure is used for the following simulations presented in this work:

• The fluid fields are initialized with analytical correlations.

• The simulation on the fluid domain starts with a coarse mesh.

• The simulation on the fluid domain continues with a fine mesh, which uses the coarse mesh solution as the initial state. Then three shock refinements according to the pressure gradient are applied

• The solid domain is finally added to the computational domain and the heat transfer and ablation simulation is performed. The simulations of the runs n • 758-759 are presented in Figure 7.4. Very good agreement is found for the 19.05mm-hemisphere (3/4 ) and 25.4mm-hemisphere [START_REF] Denis | Standard Atmosphere. National Oceanic and Atmospheric Administration[END_REF]. In contrast, the wall heat flux profile computed for the 12.7mm-hemisphere (1/2 ) is ambiguous. The result obtained with the medium mesh appears better than the one obtained with the fine mesh with shock refinements because the carbuncle is less pronounced for the medium mesh.

The stagnation point wall heat flux obtained with this mesh also agrees better with measurements, however this statement must be taken with caution since only a single experimental data point is available. The origin of this significant carbuncle in not clearly understood, nonetheless comparisons of the temperature field show a bubble of cold flow located at the stagnation point for the small hemisphere diameter. This pattern, shown in Figure 7.5 does not appear on the bigger models. The vorticity field also appears to be different for the small diameter hemisphere, as shown in Figure 7.6. In the simulations, vorticity spots are periodically generated along the shock caused. They are artificially caused by the mesh. Usually they do not have an effect on the wall heat flux, but for the 12.7mm-hemisphere (1/2 ), a vorticity spot is located close to the symmetry axis and its vorticity trail interacts with the boundary layer on the stagnation point. This vorticity pattern likely contributes to the simulated carbuncle phenomenon .

The simulation results for the high-enthalpy expansion tube runs n • 760-761 are presented in Figure 7.7. The simulated wall heat fluxes are always overestimated by up to 50%. These simulations illustrate the limitation of the modelling assuming ideal, non-reacting gas currently implemented in the solver. Indeed the gas flow for these experiments is an air flow heated up to 6000K. Hence the dissociation of the air molecules and the formation of new ones can no longer be neglected. Not taking the flow chemistry into account leads to an overestimation of the temperature in the shock layer and consequently to an overestimation of the the wall heat flux.

As a summary, the simulations of the experiments of Hollis demonstrate the validity of the local mesh refinement strategy. Interactions between the shock region and the boundary layer near the stagnation point can occur, but this is also readily detected. The comparison with the experiment featuring a high-enthalpy air flow highlights that the solver is accurate and validated only for flows without dissociation and chemical reactions.

ISL measurements

Wall heat flux measurements on a sphere have already been conducted previously at ISL, resulting in an experimental database of these measurements ( [START_REF] Berlemann | Wärmestrommessungen an einem generischen Stoßrohrmodell mit einer Mach-10-Duse[END_REF], [START_REF] Biegger | Wärmestrommessungen an Kugeln bei Machzahl 10[END_REF], [8], [START_REF] Mielke | Wärmestromdichtemessung an Flugkörpern im Hyperschall-Hochenthalpie-Stoßwindkanal. Bachelorarbeit, Institut[END_REF] & Chapter 3). In this section, the focus is on the Mach 10 experiments reproducing laminar flow conditions equivalent to earth-atmospheric altitudes of 40km, 50km, 60km and 70km. The flow conditions are given in the Table 3.1 (for the experiments performed during this thesis) and in Table 8.9 (for the previous ISL experiments). The test gas is nitrogen and an appropriate thermodynamic model has already been described in Section 5.2. The flow is considered laminar. Due to the low flow enthalpy, no dissociation of Nitrogen needs to be modelled. The simulated wall heat fluxes are shown in Figure 7.8 for the experiments performed during this thesis and are available in Appendix V for the others. The experimental and numerical wall heat fluxes are in good agreement. A strong carbuncle effect is nonetheless observed for the simulation of the 150mm-sphere at 40km altitude. The temperature and velocity fields do not show a numerical 'cold bubble' pattern such as the one presented in Figure 7.5 and 7.6. This carbuncle is likely caused by the coarse discretization of the mesh near the stagnation point.

For the flow conditions equivalent to an altitude of 70km, the simulated wall heat flux is overestimated near the stagnation region. The stagnation temperature stays below 4000K and this effect appears at low density only. Therefore the nitrogen dissociation does not explain the overestimation. On the contrary, the thermodynamic properties of the flow are likely to freeze in the nozzle during expansion, so the vibrational temperature is not at equilibrium with the translational-rotational temperature. The non-equilibrium nozzle flow may reduce the temperature of the flow before the shock and inside the shock layer, which decreases of the wall heat flux. Also, as shown by Nompelis and Candler [START_REF] Nompelis | Effect of vibrational nonequilibrium on hypersonic double-cone experiments[END_REF], simulations of shock tunnel nozzle flows neglecting thermal non-equilibrium tend to over-estimate wall heat fluxes.

The Camel effect is visible for almost all wall heat flux profiles between 40km and 60km. The position of the maximum wall heat flux is measured by the second sensor of the profile and located at 12.5 • from the stagnation point. In the simulations, the maximal wall heat flux is also shifted to the side of the stagnation point for flow conditions between 50km and 70km. The shapes of the simulated wall heat flux differs from those presented in the publication [8] and obtained with the TAU code. The discrepancies can partially be explained by differences of the inlet conditions, the thermodynamic modelling and the mesh domain. In the present work, the real shock tunnel conditions are taken, whereas the US Standard Atmosphere table was used by Adeli [8]. The differences between US Standard Atmosphere and real shock tunnel conditions exceed 25%, as shown in Appendix 8.4. The simulated wall heat fluxes show a stronger Camel effect for the lowest density conditions. This behaviour is contrary to the measurements. The simulated maximal wall heat flux is located between 5 • and 10 • , which is in agreement with the measurements.

Eventually, the origin of the Camel effect is still unknown. The present work, however, again unveils this effect in the ISL shock tunnel. The laminar simulations produce a offcenter similar shift of the location of the maximal wall heat flux. This effect seems to be dependent on the freestream density. The uncertainties in the experiments and in the simulations demonstrate the need to continue the investigation of the Camel effect. In particular, the resolution of the measurements needs to be increased by adding more thin-film sensors in the region affected by the Camel effect. It will allow more reliable comparisons with the simulations and thus more reliable conclusions concerning the Camel effect.

Turbulent flow

The Mach 4.5 experiments presented in Chapter 3 are simulated with the solver ablation-FOAM. The experimental wall heat flux suggests that a laminar-turbulent transition occurs in the vicinity of the stagnation point. The flow conditions used in the simulations are those corresponding to the Mach 4.5 flows in Table 3.1. As for the laminar flow simulations, only the fluid domain is simulated since the temperature variation on the model surface is still negligible compared to the recovery temperature. Accordingly, the wall surface temperature is set to 293K. The thermodynamic model remains unchanged from the laminar flow simulations.

The boundary layer is expected to become turbulent on the model surface, hence the turbulence is modelled in the simulation. For each flow condition, three turbulence models are tested:

• a laminar model. The turbulence is not modelled. It is described in [START_REF] Menter | Two-equation eddy-viscosity turbulence models for engineering applications[END_REF] and [START_REF] Menter | Ten years of industrial experience with the SST turbulence model[END_REF] and the OpenFOAM implementation is available in [START_REF] Opencfd | OpenFOAM API guide v2012[END_REF]. This turbulence model does not simulate the laminar-turbulent transition.

• the k -ω SST γ-Re θ model. It is a four-equation eddy-viscosity model. It is described in [START_REF] Langtry | A correlation-based transition model using local variables for unstructured parallelized CFD codes[END_REF] and [START_REF] Langtry | Correlation-Based Transition Modeling for Unstructured Parallelized Computational Fluid Dynamics Codes[END_REF] and the OpenFOAM implementation can be found in [START_REF] Opencfd | OpenFOAM API guide v2012[END_REF]. This turbulence model is able to simulate transitional flows. Based on the two-equation k-ω SST model, blending functions are implemented to switch between laminar and turbulent flow. The blending functions are dependent on two local variables: the intermittency γ and the transitional momentum thickness Reynolds number Re θ . These two parameters add two transport equations to solve.

Figure 7.9 compares the wall heat flux profiles depending on the turbulence modelling. For all conditions, the laminar flow simulations underestimate the wall heat flux up to a factor of 2-3. Only at the stagnation point and by neglecting the Camel effect, the difference between experimental and simulated stagnation point wall heat fluxes is reduced and below the measurement uncertainties.

On the contrary, the turbulent simulations with the k -ω SST model are in good agreement with the experimental data downstream of the location 12.5 • -25 • irrespective of the freestream turbulence intensity. This strongly suggests that the boundary layer becomes turbulent upstream of the maximal wall heat flux location. The stagnation point wall heat flux obtained with the k -ω SST model is slightly overestimated because it depends on the freestream turbulence in the shock tunnel, which is unknown. Also, the k -ω SST model is not suited for transitional flow simulations.

The freestream flow is not perfectly quiet as it is expected for a shock tunnel flow. The best results for the wall heat flux are obtained for turbulence intensities around 0.2-0.5%, which are low levels of freestream turbulence. This turbulence level is low, probably because the test models are small compared to the nozzle diameter and are placed in the centreline of the nozzle. The stagnation point wall heat flux is controlled by the freestream turbulence: the higher the turbulence level, the higher the stagnation wall heat flux and vice-versa. Because the freestream turbulence intensity is low, the laminar simulations provide stagnation wall heat fluxes close to the experimental values. Just downstream of the stagnation point, the boundary layer becomes turbulent near 10 • . The resulting turbulent wall heat flux is added to the total wall heat flux, which results in an increase of the the wall heat flux as shown in Section 2.2.1. The k -ω SST model is able to reproduce the location of the maximal wall heat flux correctly.

The k -ω SST γ -Re θ model is able to simulate transitional flows. It is not especially developed for highly compressible flows. However, it has been chosen because of its similarities with the k -ω SST model and because it is already implemented in OpenFOAM. The wall heat flux profiles obtained with this turbulence model are compared with the others in Figure 7 The transition phenomenon is correctly modelled. For the condition D100mm, h = 6km, the wall heat flux increase between 10 • and 25 • is accurately captured, however before this location the profile appears unreliable because of a non-physical local maximum. This artefact is probably due to a combination of the carbuncle effect and the decrease of the laminar wall heat flux. This leads to an uncertain wall heat flux in this region. The last heat flux pattern is based on the condition D100mm, h = 13km. In Figure 7.9c, the laminar-turbulent transition is triggered near 25 • , thus too far from the measured location of the transition. No carbuncle is generated for this case.

The three preceding simulations illustrate the high sensitivity of the k -ω SST γ -Re θ model to the freestream conditions. Figure 7.10 shows the results of the sensitivity analysis based on the turbulence intensity. The turbulence intensity Tu is defined as:

Tu = 100 ||u || ||u ∞ || = 100 2 3 k ∞ ||u ∞ || (7.1)
From the turbulence intensity, the freestream turbulence parameters are set as:

k ∞ = 3 2 (Tu • ||u ∞ ||) 2 (7.2) ω ∞ = k ∞ L (7.3)
With L defined as the nozzle diameter.

γ = 1 (7.4)
The freestream turbulence parameters for the transition modelling are defined as follow: Simulations with turbulence intensities from 0.05% to 0.5% are performed. For the comparison, the laminar, the k -ω SST and the experimental wall heat flux profiles are presented in the diagram. It is shown that small freestream turbulence variations completely change the wall heat flux profile. For a turbulence intensity below 0.05%, the transition is not triggered by the turbulence model and the wall heat flux remains very close to the laminar one. Contrary to that, for a turbulence intensity above 0.5% the flow is fully turbulent and the k -ω SST γ -Re θ model behaves like the k -ω SST model. Between the laminar and turbulent cases, the k -ω SST γ -Re θ model produces a boundary layer transition for the turbulence intensities 0.1%, 0.2% and 0.25%. For the cases 0.2% and 0.25%, the transition occurs before the position at 10 • in agreement with the measurements. For the turbulence intensity of 0.1%, the transition is delayed. The boundary layer remains laminar until the position at 45 • . Figure 7.10 highlights how sensitive the simulations with the k -ω SST γ -Re θ model are. A change of 0.1% in the turbulence intensity completely changes the transition location. In conclusion, the turbulence parameters in the shock tunnel are unknown and there is no efficient method to determine the proper freestream turbulence conditions other than empirically testing different values until the results match the experimental data. Moreover the transition is difficult to trigger in the case of small models, such as the gallium models used in this work. The large negative pressure gradient along the model forces the boundary layer to remain laminar or forces it to re-laminarize. For all these reasons, the k -ω SST γ -Re θ model is deemed not suitable for the transition modelling in the present case. The classical k -ω SST model is thus preferred, even though it lacks a transition model.

Re θ = 1173

Simulation of the ablation of the hemisphere-cylinder models

In the preceding parts, the verifications of the solver ablationFOAM and the validations of the wall heat flux predictions for both low and high density hypersonic flows were presented. The present part is dedicated to the simulations and analysis of the gallium ablation experiments part of this thesis. The first section presents the simulated flow fields and their impacts on the ablation process. The second part provides the simulation results related to the ablation of the gallium model. The actual limitations of the solver are also discussed.

Hypersonic flow fields around the 14mm-hemisphere-cylinder model

In this section the boundary layer around the gallium 14mm-hemisphere-cylinder model is investigated. For that purpose, the flow field considered is based on the flow field condition n • 1 of Table 4.1 (Mach 4.5, h = 3km). Four simulations are performed: one laminar and three with the turbulence model k -ω SST and freestream turbulence intensities of 0.1%, 0.2% and 0.5%, respectively. The freestream conditions for k and ω are computed by Equations 7.2 and 7.3.

Figure 8.1a presents the wall heat flux profiles in function of the turbulence intensity. No wall heat flux measurement is available for this model diameter. However the correlations presented in the Table 2.4 of the Section 2. The increase of the stagnation wall heat flux observed for the k -ω SST simulations is caused by the freestream turbulence which is convected to the stagnation point. As a consequence, the resulting turbulent wall heat flux is added to the laminar one and the total stagnation wall heat flux increases. Because the model surface temperature cannot exceed the gallium melting temperature, the radiative cooling on the model surface is estimated below 500W • m -2 with the hypothesis of black body radiation. Thus, the contribution of the radiative cooling is negligible in the thermal balance at the model surface. Although this con-clusion generally holds for low melting point materials used in shock tunnel experiments, radiative surface cooling may not be negligible during real hypersonic flights.

The wall heat flux profiles are compared to the dimensionless ablation profiles presented in Figure 4.12a. As observed in the experiments, the ablation velocity decreases slightly from the stagnation point to a position around 10 • and then rapidly increases until 30 • . This location corresponds to the point of maximal ablation. Downstream this point, the ablation intensity decreases monotonically. This ablation profile agrees well with the wall heat flux profiles for the turbulent simulations with turbulent intensities of 0.1% and 0.2%. Indeed the location for the maximal wall heat flux is consistent with the location of the maximal ablation. This demonstrates the strong coupling between the thermal heating and the ablation.

The wall shear stress profiles are shown in Figure 8.1b. As expected the wall shear stress is higher for the turbulent simulations than for the laminar ones. The profiles are nonetheless all similar and show a wall shear stress maximum located at 50 • , which is thus far from the ablation maximum at 30 • . In addition, the maximal wall shear stress value is at about 0.006M P a whereas the elastic limit of gallium is estimated between 8M P a and 25M P a and the tensile strength between 15M P a and 40M P a [4]. In conclusion, the wall shear stress profile does not match with the ablation profile and the wall shear stress value is several orders of magnitude below the mechanical limits of solid gallium. Therefore the ablation by erosion, i.e. the mechanical ablation, does not play a significant role in the ablation of the gallium model. As gallium is a rather soft metal, this conclusion is still valid for other hard metals and for flow density similar or lower than the densities used in the simulation. It is demonstrated that the thermal heating is the main contributor to the ablation process in the hypersonic dense flow. An analysis of the boundary layer is proposed for a better understanding of the the wall heat flux profile and thus of the ablation profile. Seven boundary layer profiles are extracted on lines located at 5 .3b present the velocity and temperature profiles in the boundary layer, from the wall surface to a wall distance of 100µm. The increase of the boundary layer thickness along the model surface is clearly visible in both diagrams. The velocity at the outer edge of the boundary layer increases rapidly downstream of the stagnation point. In parallel to that, the temperature decreases, however the temperature profiles from 5 • to 20 • are still close to each other.

The high velocity gradient near the wall produces turbulence. This is shown in Figures 8.3c and 8.3d, with the turbulence kinetic energy profile and the specific turbulence dissipation profile, respectively. The turbulence dissipation decreases with the wall distance from As a consequence, the profiles for the turbulent viscosity and the turbulent thermal diffusivity are altered by the turbulence increase. Below 20 • , the boundary layer profiles for these variables have only one maximum corresponding to a laminar boundary layer, while the profiles for downstream positions show a second maximum developing very close to the wall. Thus, the turbulence enhances the momentum and energy diffusion near the wall. Finally, the temperature gradient at the wall, which is proportional to the wall heat flux, also presents a local maximum near 20 • .

It is demonstrated that the transition of the boundary layer occurs between 10 • and 20 • and that the turbulent boundary layer is the cause of the maximum heat flux located near 30 • . Furthermore, the wall heat flux is strongly coupled to the ablation intensity, so the transition to the turbulence is also the cause of the characteristic ablation shape. The coupling along the wall of the turbulence increase with the temperature decrease induces that the position 30 • is the optimum position for the maximal wall heat flux.

Ablation of the hemisphere-cylinder models

Simulations performed with the ablationFOAM solver aim to predict if ablation occurs during a hypersonic flight, and if yes, how the object shape changes due to the melting. The present part presents the results obtained with the solver and compares them with the experimental and theoretical data. The simulations are performed with the 14mm-and 10mmhemisphere-cone models. The flow conditions used are the same as those presented in Table 4.1. The turbulence model k -ω SST is used in a combination with a turbulence intensity of 0.2%. Table 8.1 presents the ablation start time obtained with the simulations. The location of the ablation start is located at the same position of the maximum wall heat flux. The theoretical ablation start time is based on the semi-infinite heat transfer theory of Section 4. The theoretical, experimental and numerical ablation start times are of the same order of magnitude. For the condition h = 13km no ablation is computed by the simulation before 2ms as it is observed during the shock tunnel experiment. An interval for the ablation start time is given rather than an exact number. The first interval bound corresponds to the time of the first variation of the phase fraction field and the second bound is the time when the mesh motion starts. In theory, for a pure mixture this interval is reduced to a single time, but because the first cell height at the interface is not infinitely small, the single value becomes a time interval in the simulation. Regarding the uncertainties of each method, none of them provides results superior over the others. Table 8.3 summarizes the conditions and assumptions used for the determination of the ablation start time. In the case of gallium with its Indeed the difference between melting temperature and initial temperature is of only around 10K , thus a variation of 1K leads to an ablation start time variation of around 20%. If the temperature variation is of 4K , the initial time varies by 60%. For a hypersonic vehicle, the difference between the initial wall temperature and the melting temperature is much higher, thus the sensitivity is strongly decreased. It is important to note that the theory of heat transfer in a semi-infinite medium is valid only for short durations, i.e. small Fourier numbers. This analytical correlation is hence valid for short duration flight, i.e. up to several seconds or for low thermal diffusivity materials. only, up to the melting temperature -and the computational mesh for the fluid domain does not move at all. The Figure 8.5b shows the ablation profiles obtained with the numerical simulations. First, the overall shape is similar to the experimental ones. As already explained, no ablation is detected at t = 2ms for the flow condition h = 13km. Numerical instabilities caused by the surface recession are clearly visible on the ablation profiles. These instabilities are proportional to the surface recession speed, such that their magnitudes are higher for the conditions h = 1km and h = 3km than for the condition h = 6km. Moreover, these instabilities increase when the solid mesh is finely discretized along the moving interface. The actual boundary condition used for the surface recession calculation does not include a smoothing algorithm of the surface velocity at the interface. In other words, the surface velocity of a boundary cell does not depend on the surface velocities of its neighbouring cells. As a consequence, the smooth initial boundary surface becomes more and more rugged over time. Small surface displacements can thus be simulated with the actual solver, however large surface displacements lead to overlapping cells and the simulations to crash. In absence of a smoothing algorithm, only simulations performed on meshes with a coarse surface discretization are possible. In spite of this limitation, the comparison of dimensionless ablation profiles with their dependence on the flow conditions, shown in Figure 8.5c, are in agreement with the expected profile shapes. The ablation maximum is located at the same position for all conditions. Moreover the region affected by the ablation is larger for the case with the higher wall heat flux, as also observed experimentally in 4.7. For the condition h = 6km, the mesh displacement downstream of 45 • is caused by numerical errors, probably rounding errors, generated during the mesh motion algorithm. The numerical error is determined to a value of around 0.01mm and is observed in all simulations including a mesh motion.

For the 14mm-hemisphere, the comparisons with the experimental results are presented in Figures 8.5d and 8.5e. The computed ablation speed is of the same order of magnitude as in the experimental data. The maximum region affected by the ablation is in very good agreement between the simulation (63 • ) and the experiments (69 • ). The ablation magnitude is underestimated in the simulation compared to the experimental results. The delay in the ablation start time and in the beginning of the computed surface recession, presented in Table 8.1, are the cause of this underestimation. Hence the simulated ablation profile corresponds to an experimental profile prior the time 2ms. This discrepancy in the ablation profile between simulations and experiments could be solved with a smaller cell at the interface, however the cell overlapping problem will then arise more quickly. At the end of this study, several novel contributions have been made in the research field of liquid ablation by melting.

Firstly, a new experimental method for the measurement of liquid ablation in a hypersonic impulse facility has been demonstrated. For the first time, the ablation by melting was studied within only two milliseconds of a steady flow time generated by a shock tunnel. Gallium was used as low temperature ablator, which included an innovative manufacturing process of the model based on low temperature casting. In parallel, a measurement system based on a direct optical visualization setup and on an image post-processing algorithm has been developed. The correction of the model vibration was highlighted as a critical key point for an accurate ablation measurement. For that purpose, the use of two high-speed cameras with two magnifications has been necessary: one camera has been dedicated to the vibration correction and the other one recorded an accurate ablation profile centered on the stagnation point region. The study of ablation by melting in the ISL shock tunnel was therefore possible, however restricted by the steady flow duration, such that only gallium provided useful results and only the high density flow conditions generated ablation within the two milliseconds. However, the flow conditions studied, with their high densities and high mass flow rates, are difficult to reproduce in continuous or long duration hypersonic facilities. The uncertainties of the shock tunnel experiments remain high and require time and space averaging. However, the measured ablation profiles confirm the previous observations made in the hyperballistic tunnel. In a turbulent flow, the ablation maximum is not located at the stagnation point but rather at a position near 30 • . The laminar-turbulent transition is located near 10 • -20 • depending on the flow condition but always before the point of maximum ablation. The characteristic ablation profile for a turbulent flow is therefore similar to the wall heat flux profile with a local ablation maximum at the stagnation point and an ablation intensity decay until the laminar-turbulent transition point. Then, the ablation intensity rises quickly until the ablation maximum near 30 • . Downstream of this location, the ablation intensity decreases monotonically. No ablation is detected on the cylinder or on the cone after the hemisphere. The numerical simulations of the flow field around the model confirm that the contribution of the ablation by erosion is negligible in the overall ablation process in high density flows.

The ablation profiles are controlled by the wall heat flux profile. Therefore, the wall heat flux profiles are measured in the shock tunnel under different flow conditions. At low density, the so-called Camel effect has been observed again, as it had been in previous investigations. The sensor limitations and the Carbuncle effect in the simulations do not allow to explain this shift of the maximum wall heat flux away from the stagnation point at low density. Up to now, the Camel effect has only been observed in the ISL shock tunnel, so we cannot exclude that it is caused by the facility during the hypersonic flow generation. Non-equilibrium effects, thermochemical effects or turbulence contributions are likely candidates that could explain the Camel effect. Finally, the wall heat flux profiles for the low density flows were taken as reference profiles for the validation of the numerical simulations. The wall heat flux measurements in the high density flows demonstrated that the boundary layer becomes turbulent upstream of the point of maximal wall heat flux, as also concluded based on the analysis of the ablation experiments. Thus, a characteristic wall heat flux profile is defined which is very similar to the characteristic ablation profile. The results from the wall heat flux experiment strongly suggest that the knowledge of the turbulence state is a key factor for the design of a hypersonic vehicle at low altitude. In the numerical simulation, the laminarturbulent transition is modelled with the k -ω SST γ -Re θ model. However, this model is not specifically dedicated to hypersonic flows and results show an excessive sensitivity to the freestream turbulence. That is why the k -ω SST model has been preferred for the ablation simulation despite the lack of a transition model.

The final contribution of this work has been the programming and the validation of the ablationFOAM solver dedicated to modelling the ablation by melting. All the programmed sub-modules of the solver have passed the validation on different test cases successfully. The global behaviour of the solver has been validated with both laminar and turbulent experiments. The solver is now validated for the simulation of liquid ablation in a hypersonic flow, which can include molecular vibration and turbulence. The actual solver does not support chemical reactions, non-equilibrium effects and ionization. In conclusion, the solver is suitable for hypersonic flight at comparatively low speed in a dense atmosphere. However exactly these conditions are generally encountered by hypersonic projectiles or hypersonic cruise missiles. After all, this work opens new perspectives for the wall heat flux and the ablation studies at ISL.

For the wall heat flux measurements, the regions of interest for both the Camel effect and the laminar-turbulent transition are now known. The future experimental models will be upgraded with more sensors in the regions close to the stagnation point.

For the ablation experiment, the experimental setup using gallium is now mature and readily applicable and the gallium experiments can be reproduced with more complex geometries such as projectile fins, for example. The measurements uncertainties can be further reduced, especially by the measurement of the initial model surface temperature.

In addition to this work, the ISL hyperballistic tunnel and the light gas guns have been modernized in order to conduct ablation experiments in hypersonic free flight. All the preparatory work has been finished and the tunnel is now ready for the experiments to commence. Based on Luneau's work and also the present work, the ablation of different metals such as steel, aluminium, aluminium alloy, aluminium oxide, titanium and Inconel will be investigated at speed between 2 and 9 km • s -1 .

As far as numerical simulations are concerned, the solver can now be used for the simulation of ablation of real hypersonic objects. The next task will be the inclusion of the modules for the thermochemical reactions and for the thermal radiation.

Last but not least, the study of ablation by melting is again an important research topic within modern hypersonics. In this context, this work demonstrates that the combination of experiments conducted in hypersonic facilities with the numerical simulations is essential.

Part V

Appendices

Altitude estimation for the flow conditions

Each flow condition is defined by a Mach number and an altitude. This altitude is computed from the freestream pressure and according to the US Standard Atmosphere model of 1976 [START_REF] Denis | Standard Atmosphere. National Oceanic and Atmospheric Administration[END_REF]. The designation of an equivalent altitude per flow condition is chosen for convenience and for the sake of simplicity but the reader should keep in mind that the experimental flow conditions can still deviate from the theoretical ones given by the US Standard Atmosphere model. Table 8.4 summarizes the deviation between the experimental and theoretical flow conditions. The deviation is computed as:

∆Φ = φ exp -φ t

Notes on the uncertainty calculation

Unless otherwise stated, all measurement uncertainties shown in the Figures are statistical uncertainties computed with a 95% confidence interval centred around the mean value. Thus, we assume a normal distribution of the measurement error, such that the uncertainty is computed as follows:

± φ = St i • σ (8.2)
The Student coefficient St i is a function of the number of measurements i and σ is the standard deviation. It is given in Table 8.5 It is important to note that the Student coefficient increases rapidly if the number of measurements is small, thus the uncertainty increases at the same rate. As a consequence, in the case of missing data or of a low number of tests the uncertainty can be systematically overestimated.

Raw data for the wall heat flux measurement experiments of Chapter 3

Material ρ C p λ T f T i (-) (J • kg -1 • K -1 ) (W • m -1 • K -1 ) (K ) (K )
q vol = 1E 8W • m -3 T w = 273K λ = 400W • m -1 • K -1 C p = 370J • kg -1 • K -1 ρ = 5904kg • m -3
Test case n • 2 / Steady conduction / Tube / Fixed temperature and fixed gradient g d z T w q z d z = 0.16m g = 0.3m

T w = 1000K qz = -4E + 6W • m -2 λ = 400W • m -1 • K -1 C p = 370J • kg -1 • K -1 ρ = 5904kg • m -3
T w = 1000K T z = 1000K λ = 400W • m -1 • K -1 C p = 370J • kg -1 • K -1 ρ = 5904kg • m -3
Test case n • 4 /Steady conduction / 1D-beam / Fixed temperature and fixed gradient

L T z q z L = 0.1m qz = 4E + 6W • m -2 T z = 273K λ = 400W • m -1 • K -1 C p = 370J • kg -1 • K -1 ρ = 5904kg • m -3
Test case n • 5 / Steady conduction / 1D-beam / Fixed temperature

L T z T w L = 0.1m T w = 1000K T z = 273K λ = 400W • m -1 • K -1 C p = 370J • kg -1 • K -1 ρ = 5904kg • m -3
Test case n • 6 / Steady conduction / 1D-beam / Heat conduction into two identical domains

L domain 1 T z T w L domain 2 L = 0.1m L d omai n 2 = 1m T w = 1000K T z = 273K λ = 400W • m -1 • K -1 C p = 370J • kg -1 • K -1 ρ = 5904kg • m -3
Test case n

L = 0.1m L d omai n 2 = 1m T w = 1000K T z = 273K λ d omai n 1 = 400W • m -1 • K -1 λ d omai n 2 = 40W • m -1 • K -1 C p = 370J • kg -1 • K -1 ρ = 5904kg • m -3
Test case n • 8 / Transient conduction / Sphere / Quasi constant temperature heat conduction

d T 0 T L = 1m T 0 = 288K T ∞ = 333K h = 20W • m -2 • K -1 λ = 63.9W • m -1 • K -1 C p = 434J • kg -1 • K -1 ρ = 7832kg • m -3
Test case n • 9 / Transient conduction / Semi-infinite medium / Fixed convective heat transfer coefficient

L T 0 h T L = 1m T 0 = 273K T ∞ = 1000K h = 200W • m -2 • K -1 λ = 400W • m -1 • K -1 C p = 370J • kg -1 • K -1 ρ = 5904kg • m -3
Test case n • 10 / Transient conduction / Semi-infinite medium / Fixed temperature

L T 0 T w L = 1m T 0 = 273K T w = 1000K λ = 400W • m -1 • K -1 C p = 370J • kg -1 • K -1 ρ = 5904kg • m -3
Test case n • 11 / Transient conduction / Semi-infinite medium / Fixed gradient

L T 0 q z L = 1m T 0 = 273K qz = 4E + 6W • m -2 λ = 400W • m -1 • K -1 C p = 370J • kg -1 • K -1 ρ = 5904kg • m -3
Test case n • 5 / 2D-beam / Solid mesh motion

L = 1m T w,1 = 273K T w,2 = 300K T f = 273K H f = 1000J • kg -1 λ = 100W • m -1 • K -1 C p = 100J • kg -1 • K -1 ρ = 1000kg • m -3
Test case n • 2 / Ablation / 1D-beam / Fixed temperature, fixed gradient and surface L q w,2 T w,1

L = 1m q w,2 = 2E + 4W • m -3 T w,1 = 263K T f = 273K H f = 1000J • kg -1 λ = 100W • m -1 • K -1 C p = 100J • kg -1 • K -1 ρ = 1000kg • m -3
CFD validations with the wall heat flux data of Chapter 7

Table 8.9: Flow conditions for the simulations of the ISL low density experiments prior to this work ( [START_REF] Berlemann | Wärmestrommessungen an einem generischen Stoßrohrmodell mit einer Mach-10-Duse[END_REF], [START_REF] Biegger | Wärmestrommessungen an Kugeln bei Machzahl 10[END_REF], [8]) 2D-axisymmetric mesh around an hemisphere Figure 8.12: View of the meshes used for the hemisphere simulations. On the left: initial mesh, on the center : mesh with one shock refinement, on the right: mesh with two shock refinements

spher e M ∞ Re ∞ h u ∞ T ∞ P ∞ T 0 h 0 P ∞ (mm) (-) (-) (km) (m • s -1 ) (K ) (P a) (K ) (J • kg -1 ) (W • m -2 )

6. 3

 3 Velocity profiles at t = 8E -6s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4 Details of the pressure profiles at t = 8E -6s . . . . . . . . . . . . . . . . . . . . . 6.5 Details of the temperature profiles at t = 8E -6s . . . . . . . . . . . . . . . . . . . 6.6 Details of the velocity profiles at t = 8E -6s . . . . . . . . . . . . . . . . . . . . . 6.7 The meshes used for the 15 • -wedge . . . . . . . . . . . . . . . . . . . . . . . . . . 6.8 Mach number fields for the Mach 3 flow . . . . . . . . . . . . . . . . . . . . . . . 6.9 Numerical and analytical results for the steady state heat transfer problems . . 6.10 Numerical and analytical results for the multi-domain and transient problems 6.11 View of the mesh motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.12 Results for the melting case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.13 Results for the ablation case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.1 Wall heat flux profiles obtained for the condition n • 307 002-003 . . . . . . . . . 7.2 Wall heat flux profiles obtained for the condition n • 307 004-005 . . . . . . . . . 7.3 Wall heat flux profiles obtained for the condition n • 307 006-007 . . . . . . . . . 7.4 Wall heat flux profiles obtained for the condition n • 758-759 . . . . . . . . . . . 7.5 Temperature field in the stagnation region . . . . . . . . . . . . . . . . . . . . . . 7.6 Vorticity field in the stagnation region . . . . . . . . . . . . . . . . . . . . . . . . . 7.7 Wall heat flux profiles obtained for the condition n • 760-761 . . . . . . . . . . . 7.8 Wall heat flux profiles at Mach 10 for the 100mm-hemisphere . . . . . . . . . . 7.9 Effect of the turbulence models on the simulated wall heat flux at Mach 4.5 . . 7.10 Sensitivity analysis to the turbulence intensity, Mach 4.5, D100mm, h = 3km . 8.1 Profiles for the 14mm-hemisphere-cylinder and for the condition h = 3km . . . 8.2 Location of the boundary layer profiles . . . . . . . . . . . . . . . . . . . . . . . . 8.3 Boundary layer profiles for the 14mm-hemisphere-cylinder, h = 3km . . . . . . 8.4 Evolution of the wall temperature gradient along the model . . . . . . . . . . . . 8.5 Numerical results with the ablation module obtained at t = 2ms . . . . . . . . . 8.6 Wall heat flux profiles at Mach 10 for the 100mm-sphere . . . . . . . . . . . . . . 114 8.7 Wall heat flux profiles at Mach 4.5 for the 100mm-sphere . . . . . . . . . . . . . 115 8.8 Ablation profiles for the condition Mach 4.5, h = 3km . . . . . . . . . . . . . . . 121 8.9 Seiler experiment at Mach 10, D = 150mm . . . . . . . . . . . . . . . . . . . . . . 132 8.10 Seiler experiment at Mach 10, D = 100mm . . . . . . . . . . . . . . . . . . . . . . 133 8.11 Seiler experiment at Mach 10, D = 50mm . . . . . . . . . . . . . . . . . . . . . . . 134 8.12 View of the meshes used for the hemisphere simulations . . . . . . . . . . . . . 135

  u or u or v = Velocity [m • s -1 ]x or y

  = Stagnation value = Initial value δ = Condition at outer edge of the boundary layer ∞ Measurements & Simulations scientific group ASM P A = Missile Air-Sol Moyenne Portée-Amélioré (Improved Medium-Range Air to Surface missile) C E P = Circular Error Probable C F D = Computational Fluid Dynamics C F L = Courant-Friedrich-Lewy Number C N E S = Centre National d'Etudes Spatiales (National Centre for Space Studies) DSMC = Direct Simulation Monte-Carlo E S A = European Space Agency G AMG = Geometric Agglomerated Algebraic Multigrid solver GCC = GNU Compiler Collection GP L = GNU General Public License HC M = Hypersonic Cruise Missile HGV = Hypersonic Glide Vehicle HV P = Hypervelocity Projectile IC B M = Inter Continental Ballistic Missile I SL = Institut franco-allemand de recherches de Saint-Louis (French-German Research Institute of Saint-Louis) LEO = Low Earth Orbit LI F = Laser-Induced Fluorescence LT S = Long Term Support M RB M = Medium Range Ballistic Missile M SI I = Mathématiques, Sciences de l'Information et de l'Ingénieur (Mathematics, Sciences of the Information and of the Engineer) MU SIC = MUlti Simulator In Combination N AC A = National Advisory Committee for Aeronautics N AS A = National Aeronautics and Space Administration N SM B = Navier-Stokes Multi-Block ON E R A = Office National d'Etudes et de Recherches Aérospatiales (The French Aerospace Lab) OpenCV = Open Computer Vision OpenF O AM = Open Field Operation And Manipulation P AT O = Porous material Analysis Toolbox based on OpenFOAM P BiCG = Preconditioned Bi-Conjugate Gradient PE A = Programme d'Etude Amont (Prospective Reasearch Programme) P I M P LE = Pressure Implicit Method for Pressure-Linked Equations RW T H = Rheinisch-Westfälische Technische Hochschule Aachen (Aachen University) SC AR AB = Spacecraft Atmospheric Re-entry and Aerothermal Break-up SC L = Space Conservation Law SST = Shear Stress Transport ST A = ISL Shock Tunnel A ST B = ISL Shock Tunnel B T P S = Thermal Protection System U M R = Unité Mixte de Recherche (Joint Research Unit) U S AF = United States Air Force U S or U S A = United States of America

6 CHAPTER 1 .Figure 1 . 1 :Figure 1 . 2 :

 611112 Figure 1.1: Geometry comparison between a reentry capsule and a warhead

Figure 1 . 3 :

 13 Figure 1.3: Cruise domains of hypersonic weapon systems

Figure 2 . 1 :

 21 Figure 2.1: Ablation of a sphere in a hypersonic flow

Figure 2 . 2 :

 22 Figure 2.2: Heat transfer coefficient as a function of local Reynolds and Mach numbers for an insulated flat plate in free flight (from [100])

Figure 2 . 3 :

 23 Figure 2.3: Stagnation wall heat flux computed with ρ = 1.2kg • m -3 & R = 0.01m

Figure 2 . 4 :

 24 Figure 2.4: Evolution of the function f 1 (θ)

Figure 2 . 5 :

 25 Figure 2.5: Energy transfers at the wall surface

Figure 2 . 6 :

 26 Figure 2.6: Duralium spheres after a free flight in nitrogen (left) and in air (right) (u 0 = 5200m• s -1 & ρ ∞ = 1.7kg • m -3 )

Figure 2 . 7 :Figure 2 . 8 :

 2728 Figure 2.7: Schematic of the ISL shock tunnel

3. 1 . 1

 11 Sphere with heat flux sensorsThe model geometry is a sphere of diameter 100mm made of brass. It mounted on a sting and placed in the test section. 15 thin-film heat transfer sensors are installed along an equatorial plane on the model wall. All sensors are spaced 12.5 • apart and the sensor positions cover a range from -37.5 • to 137.5 • with respect to the stagnation point.

  Figure 3.1: View of the model and the sensor

Figure 3 . 2 :

 32 Figure 3.2: Wall heat flux profiles

Figure 3 . 3 :

 33 Figure 3.3: Position of the maximal wall heat flux

Figure 3 . 4 :

 34 Figure 3.4: Comparisons of wall heat flux profiles

Figure 3 . 5 :3

 35 Figure 3.5: Comparison of profiles for the 100mm-sphere at Mach 10

Figure 3 . 6 :

 36 Figure 3.6: Comparison of dimensionless wall heat flux profiles

Figure 4 . 1 .

 41 They are mounted on a rearwardfacing sting. At the flow start, the model and the sting are at the same temperature. Due to the short flow duration, no thermal conduction occurs between the gallium model and the sting.

  1: 8mm-, 10mm-and 14mm-hemisphere model geometries Due to the low melting temperature of gallium the model cannot be machined. Thus a custom model generation technique based on low temperature casting has been developed.

  Figure 4.2: Model generation tools

Figure 4 . 3 :Figure 4 . 5 :

 4345 Figure 4.3: Optical setup

  , 4.6, 4.7 and 4.8. The model sting is not ablated by the flow, so the pattern geometry does not change during the experiment. Thus, for each frame of the low-magnification video these patterns are cross-correlated with those of the reference image. The result provided by the crosscorrelation corresponds to the displacement of the model caused by sting vibration and is used to correct the model vibration. The correction is applied to the both video recordings. For the high magnification video, the displacement is multiplied by the image magnification ratio.

  The post-treatment uncertainty is one pixel when the surface visualization is undisturbed. During the test, ablation trails flow around the model and partially shroud the model surface, inducing additional errors in the model vibration correction and consequently in the calculation of the ablation profiles. To compensate for this uncertainty, the ablation profiles shown in the following sections are temporally averaged over seven consecutive frames. Additionally, they are spatially averaged, retaining a resolution of 0.05mm.

Figure 4 . 6 :

 46 Figure 4.6: Images from the post-processing

Figure 4 . 7 :Figure 4 . 8 :

 4748 Figure 4.7: 10mm-hemisphere-cone models before and after testing at Mach 4.5

Figure 4 . 9 :

 49 Figure 4.9: Comparison for the analytical and measured heating phase results

Figure 4 . 10 :

 410 Figure 4.10: Ablation profiles for the 8mm-hemisphere at the flow condition h = 3km

Figure 4 . 11 :

 411 Figure 4.11: Ablation profiles at t = 2ms for condition h = 3km with the uncertainties represented by the blue areas

Figure 4 . 12 :

 412 Figure 4.12: Comparison of ablation profiles at t = 2ms

Figure 4 . 13 :

 413 Figure 4.13: Comparison of ablation profiles obtained with various facilities and model materials

Figure 5 . 1 :

 51 Figure 5.1: Simplified flow chart of the ablationFOAM solver

Figure 5 . 2 :

 52 Figure 5.2: Enthalpy -Temperature diagram

14 )Figure 5 . 3 :

 1453 Figure 5.3: Enthalpy -Phase fraction diagram

  Figure 5.4: The gas-solid boundary

α e 2 = e 2 nbr e 2 nbr + e 2 P

 22

Figure 5 . 5 :

 55 Figure 5.5: Schematics of the first cell during the mesh motion

Figure 6 . 1 :Figure 6 . 2 :Figure 6 . 3 :

 616263 Figure 6.1: Pressure profiles at t = 8E -6s

Figure 6 . 5 :

 65 Figure 6.5: Details of the temperature profiles at t = 8E -6s

Figure 6 . 6 :

 66 Figure 6.6: Details of the velocity profiles at t = 8E -6s

Figure 6 . 7 :Figure 6 . 8 :

 6768 Figure 6.7: The meshes used for the 15 • -wedge (bottom-left: coarse mesh, upper-left: medium mesh, bottom-right: fine mesh, upper-right:: fine mesh with shock refinements)

Figure 6 .

 6 10 presents the results for the transient, multi-domain heat transfer cases. The communication between two adjacent domains is validated with the test cases n • 6 and n • 7. Two 1D-beams are connected at their ends. The contact temperature and the thermal gradients are shared as explained in Section 5.4. The data shown in Figures 6.10a and 6.10b show a perfect agreement with the theoretical temperature profile, which validates the solver capability for the multi-domain simulation. The test cases from n • 8 to n • 11 are used for testing the transient capability of the solver. Simulations with low Biot number (case n • 8) and with semi-infinite medium (case n • 9, n • 10 and n • 11) are presented. The transient temperature profiles are correctly predicted by the solver.

5 Figure 6 . 9 :

 569 Figure 6.9: Comparison of the numerical and analytical results for the steady state heat transfer problems

11 Figure 6 . 10 :Figure 6 . 11 :

 11610611 Figure 6.10: Comparison of the numerical and analytical results for the multi-domain and the transient heat transfer problems

Figure 6 . 12 :

 612 Figure 6.12: Results for the melting case

Figure 6 . 13 :

 613 Figure 6.13: Results for the ablation case

Figure 7 . 1 :Figure 7 . 2 :

 7172 Figure 7.1: Wall heat flux profiles obtained for the condition n • 307 002-003

Figure 7 . 3 :

 73 Figure 7.3: Wall heat flux profiles obtained for the condition n • 307 006-007

Figure 7 . 4 :Figure 7 . 5 :

 7475 Figure 7.4: Wall heat flux profiles obtained for the condition n • 758-759

Figure 7 . 6 :Figure 7 . 7 :

 7677 Figure 7.6: Vorticity field in the stagnation region. From the left to the right: 25.04mmhemisphere (1 ), 19.05mm-hemisphere (3/4 ), 12.7mm-hemisphere (1/2 )

Figure 7 . 8 :

 78 Figure 7.8: Wall heat flux profiles at Mach 10 for the 100mm-hemisphere

. 9 .

 9 Downstream of 30 • -40 • , the results from the k -ω SST γ -Re θ model almost coincide with those from the k -ω SST model. From his position downstream, the flow is fully turbulent and the k -ω SST γ -Re θ model behaves like the k -ω SST model. Upstream 30 • -40 • , the three different heat flux patterns are shown as results obtained by the k -ω SST γ-Re θ model. For the condition D100mm, h = 3km, the wall heat flux profiles are close to the measurements with the exception of a Carbuncle effect at the stagnation point.

Figure 7 . 9 :

 79 Figure 7.9: Effect of the turbulence models on the simulated wall heat flux at Mach 4.5

Figure 7 . 10 :

 710 Figure 7.10: Sensitivity analysis to the turbulence intensity, Mach 4.5, D100mm, h = 3km

  2.1 give a stagnation wall heat flux between 4 • 10 6 and 9 • 10 6 W • m -2 . Previous comparisons with experimental and numerical data indicate that the stagnation wall heat flux estimations with the best agreement were based on the Oertel and the Seiler correlations. For the present flow condition, they predict a stagnation wall heat flux between 6•10 6 and 7•10 6 W •m -2 . Thus the stagnation wall heat flux computed in the laminar simulation around 7.5 • 10 6 W • m -2 appears realistic. Despite the Fay & Riddel and the van Driest expressions being more elaborate, they frequently underestimate the stagnation wall heat flux for the shock tunnel experiment. These correlation are rather suited for quiescent flows but not necessarily for noisy shock tunnel flows. In the present case the Fay & Riddel and the van Driest expressions give a stagnation wall heat flux between 4 • 10 6 and 5 • 10 6 W • m -2 .

Figure 8 . 1 :

 81 Figure 8.1: Profiles for the 14mm-hemisphere-cylinder and for the condition h = 3km

Figure 8 . 2 :

 82 Figure 8.2: Location of the boundary layer profiles (white lines)

Figure 8 . 3 :

 83 Figure 8.3: Boundary layer profiles for the 14mm-hemisphere-cylinder, h = 3km, turbulence intensity 0.2%

Figure 8 . 4 :

 84 Figure 8.4: Evolution of the wall temperature gradient along the model

Figure 8 .

 8 Figure 8.5 summarizes the simulation results obtained with the ablation module at the time t = 2ms. Figure 8.5a presents the wall heat flux profiles along the model. For the conditions h = 1km, h = 3km and h = 6km the wall heat flux maximum is still located near 30 • . The wall heat flux profiles at 2ms do not change significantly compared to the steady state simulations presented earlier. Indeed, the surface temperature increases by several Kelvin

Figure 8 . 5 :

 85 Figure 8.5: Numerical results with the ablation module obtained at t = 2ms

Figure 8 . 6 :Figure 8 . 7 :

 8687 Figure 8.6: Wall heat flux profiles at Mach 10 for the 100mm-sphere

Figure 8 . 8 :. 3 Test case n • 1 /

 8831 Figure 8.8: Ablation profiles for the condition Mach 4.5, h = 3km

Test case n • 3 /

 3 Steady conduction / Tube / Fixed temperature

• 7 /

 7 Steady conduction / 1D-beam / Heat conduc-

L. 5 Test case n • 1 /

 51 = 1m × 0.1m × 0.1m T w = 273K Compression: u w,1,∥ = 0.5m• -1 Extension: u w,2,∥ = 0.5m• -1 t f = 10sInitial and boundary conditions for the CFD test cases "Stephan problem" in Section 6

Figure 8 . 9 :Figure 8 . 10 :Figure 8 . 11 :

 89810811 Figure 8.9: Seiler experiment at Mach 10, D = 150mm

Table 2 .

 2

	1: Coefficients for equation 2.11 for nitrogen (N 2 ) [60]
	N 2	a 0	a 1			a 2	a 3	a 4
	For 200K < T < 1000K	3.531 -1.237E-4 -5.030E-7 2.435E-9 -1.409E-12
	For 1000K < T < 6000K 2.953 1.397E-3 -4.926E-7 7.860E-11 -4.608E-15
		ln µ = A • ln T +	B T	+	C T 2 + D	(2.12)
		ln λ = E • ln T +	F T	+	G T 2 + H	(2.13)

Table 2 .

 2

	4: Coefficients for the correlation 2.28
	Coefficients	C 1	n	m	N	M
	Laminar wall heat flux 5.81E -5 1/2 0.82 0.495 3.09
	Turbulent wall heat flux 3.25E -4 1/6 1.2 0.812 3.27

Table 4 .

 4

	h	Model geometry	Number of runs	Re D
	(km)	(-)	(-)	(-)
	1	Hemisphere 10mm-cone	2	9.06 E+5
	3	Hemisphere 8mm-cone	6	3.96 E+5
	3	Hemisphere 10mm-cone	8	4.95 E+5
	3	Hemisphere 14mm-cylindre	6	6.92 E+5
	6	Hemisphere 10mm-cone	3	5.86 E+5
	13	Hemisphere 10mm-cone	3	1.73 E+5

3: Test summary

Table 4 .

 4

		4: Heating phase durations
	h	Model geometry	Heating phase duration
	(km)	(-)	(ms)
	1	Hemisphere 10mm-cone	0.42
	3	Hemisphere 8mm-cone	0.50
	3	Hemisphere 10mm-cone	0.47
	3	Hemisphere 14mm-cylindre	0.39
	6	Hemisphere 10mm-cone	0.89
	13	Hemisphere 10mm-cone	n/a

Table 4 .

 4

		5: Freestream conditions for the ablation experiments
	Experiment	Ablated model	Facility	M ∞	Re ∞
	(-)	(-)	(-)	(-)	(-)
	Present work	Gallium hemisphere	Shock tunnel	4.5	[5E + 5; 1E + 6]
	Luneau [59]	Aluminium sphere	Hyperballistic tunnel [9; 16] [3E + 6; 4E + 6]
	Silton & Goldstein [89] Water ice hemisphere	Wind tunnel	5	1E + 6

Table 6 .

 6

	1: Initial conditions for the Riemann problem
	Shock tube part	Tube length	P i ni t	T i ni t	U i ni t
	(-)	(m)	(P a)	(K ) (m • s -1 )
	High pressure part	0.2	1.16E + 8 5000	0
	Low pressure part	0.2	6.93E + 6 3000	0

  Three additional simula-

	p (Pa)	3.00 3.25 3.50 3.75 4.00	1e7						100 cells, Co 0.1 100 cells, Co 0.01 100 cells, Co 0.001 1000 cells, Co 0.1 1000 cells, Co 0.01 1000 cells, Co 0.001 2000 cells, Co 0.1 2000 cells, Co 0.01 2000 cells, Co 0.001 Exact
		2.75
		2.50
		2.25
		2.00	0.26	0.28	0.30	0.32 x (m)	0.34	0.36	0.38	0.40
			Figure 6.4: Details of the pressure profiles at t = 8E -6s
			0.26	0.28	0.30	0.32 x (m)	0.34	0.36	0.38	0.40

Table 6 . 3

 63

: Relative errors with respect to the analytical solution

Table 6 .

 6 4: Test cases for the module dedicated to the heat transfer in solid

	n •	Test case	Geometry	Heat transfer condition
	1	Steady conduction	Cylinder	Heat source and constant wall temperature
	2	Steady conduction	Tube	Fixed temperature and fixed gradient
	3	Steady conduction	Tube	Fixed temperature
	4	Steady conduction	1D-beam	Fixed temperature and fixed gradient
	5	Steady conduction	1D-beam	Fixed temperature
	6	Steady conduction	1D-beam	Heat conduction into two identical domains
	7	Steady conduction	1D-beam	Heat conduction into two different domains
	8 Transient conduction	Sphere	Quasi constant temperature heat conduction
	9 Transient conduction Semi-infinite medium	Fixed convective heat transfer coefficient
	10 Transient conduction Semi-infinite medium	Fixed temperature
	11 Transient conduction Semi-infinite medium	Fixed gradient

Table 6 .

 6

		6: Test cases for the Stephan problem
	n • Test case Geometry	Heat transfer condition
	1

  .51 -589.428Tu + 0.2196

		1e6
	Wall heat flux qw (W m 2 )	2 4 6 8					CFD laminar CFD k SST CFD k SST Re Exp
		0
		0	10	20	30	40 Angle ( ) 50	60	70	80	90
									331.5 (Tu-0.5658) 0 .671	Tu	if Tu ≤ 1.3 if Tu > 1.3	(7.5)

  • , 10 • , 20• , 30 • , 40 • , 50 • and 60 • from the symmetry axis, as shown in Figure 8.2. Figures 8.3a and 8

Table 8 .

 8

					1: Ablation start time
	Condition	h	D	Ablation	Location of	Experimental Theoretical
					start time ablation start	start time	start time
	n •		(km) (mm)	(ms)	( • )		(ms)	(ms)
	1		3	14	0.1-0.7	25		0.4	up to 0.3
	1		3	10	0.1-0.3	26		0.4	up to 0.3
	2		1	14	0.1-0.7	25		≈ 0.4	up to 0.2
	2		1	10	0.11	24		0.4	up to 0.2
	3		6	14	0.3-1.5	25		≈ 0.9	up to 0.6
	4		13	14	>2	-		n/a	up to 2.1
				Table 8.2: Ablation parameters after 2ms
	Condition	h	D		Position of	Maximum ablation Maximal downstream
				ablation maximum	recession	position with ablation
	n •	(km) (mm)	( • )		(mm)	( • )
	1	3	14		25		0.032	62
	2	1	14		24		0.029	60
	3	6	14		22		0.006	42

Table 8 . 3 :

 83 Methods for the ablation start time determination

	Conditions	CFD simulation	Experiment	Analytic
	Assumptions	Pure melting	None	Semi-infinite heat transfer
	Flow condition	Constant	Not constant	Constant
			Pressure variation of 15%
	Stationary flow	Precisely known	Not precisely known	Exactly known
	start time		Determined through test video
	Initial wall	Set to 293K	Not precisely known	Set to 293K
	temperature
	Ablation start	Highly sensitive to the	Not precisely known	Highly sensitive to the
	time	initial wall temperature Determined through test video	initial wall temperature

Table 8 .

 8 heo 4: Experimental and theoretical flow condition deviations for the Table3.1

			φ exp			(8.1)
	Condition M ∞	h	∆P ∆ρ ∆T
	n •	(-) (km) (%) (%) (%)
	1	9.43	70	0	19	27
	2	9.65	60	0	27	42
	3	9.68	50	0	27	42
	4	4.50	13	0	24	17
	5	4.51	6	0	7	4
	6	4.42	3	0	38	23

Table 8 .

 8 5: Student coefficients used for the uncertainty calculation Number of measurements i Student coefficient St

	(-)	(-)
	2	12.71
	3	4.30
	4	3.18
	5	2.78
	6	2.57
	7	2.45
	8	2.37
	9	2.31
	10	2.26
	. . .	. . .
	∞	1.96

Table 8 .

 8 6: Flow condition per run

	Run	M ∞	h	Re x	u ∞	T ∞	P ∞	T 0	h 0
	(-)	(-)	(km)	(m -1 )	(m • s -1 )	(K )	(P a)	(K ) (J • kg -1 )

Table 8 .

 8 7: Parameters and uncertainties used for the ablation start time calculation[4]

Acknowledgements

I would like to express my warmest thanks to my supervisors Prof. Yannick HOARAU and Prof. Robert MOSE for having accepted the supervision of my work at the ICube laboratory of the University of Strasbourg. I am also particularly grateful to my supervisor at the French-German research institute of Saint-Louis (ISL), Dr. Robert HRUSCHKA for his continuous support and for the freedom of experimentation he gave me. Finally, this work would not have been possible without the trust given to me by Dr. Friedrich LEOPOLD, head of the Aerodynamics, Measurements & Simulations group at the ISL and I would like to express my gratitude to him. I would like to express my deep thanks to the members of the jury who accepted to evaluate my work. My deepest gratitude goes to Prof. Oliver CHAZOT, head of the Aeronautics and Aerospace Department at the von Karman Institute for Fluid Dynamics and Prof. Christian MUNDT,

Part IV

Conclusion

Composition of the Gallium model of Section 4.1.3 Test case n • 2 / 1D-beam / Mesh deformation normal to the wall

Test case n • 3 / 2D-beam / Mesh deformation tangent to the wall

Test case n • 4 / 2D-beam / Mesh deformation normal-tangent to the wall

Publications and communications

Papers:

Flavien DENIS

Study of the ablation by melting in hypersonic flow

Résumé

La dernière décennie a vu l'essor de nouvelles applications dans le domaine du vol hypersonique. Ainsi l'extension du domaine de vol des véhicules hypersoniques nécessite de nouvelles études sur l'échauffement aérodynamique et sur l'ablation des matériaux.

La présente étude porte sur la détermination des profils de flux thermique et d'ablation par fusion pour des matériaux métalliques en vol hypersonique en atmosphère dense. Les profils de flux thermique et d'ablation sont mesurés en soufflerie à choc. Un protocole expérimental, combinant métal ablatif à bas point de fusion et visualisation optique, permet la mesure de l'ablation dans des écoulements hypersoniques en quelques millisecondes. L'étude expérimentale met en évidence un profil caractéristique de l'ablation par fusion en écoulement turbulent.

L'étude numérique présente et valide le solveur ablationFOAM dédié à la prédiction de l'échauffement et de l'ablation en écoulement hypersonique. Le rôle du développement de la couche limite turbulente dans le phénomène d'ablation est explicité. Enfin, l'utilisation de différents modèles de turbulence est également discutée.

Mots-clés : écoulement hypersonique, échauffement, ablation, tunnel à choc, CFD

Abstract

The last decade has seen the development of new applications in the field of hypersonic flight. Thus, the flight domain extension of hypersonic vehicles requires new studies on aerothermal heating and on the thermal ablation.

The present study focuses on the determination of wall heat flux and ablation profiles for metallic materials in dense hypersonic flow. The heat flux and ablation profiles are measured in a shock tunnel. An experimental protocol, combining low melting point ablative metal and optical visualisation, allows the ablation measurement with a hypersonic impulse facility within the milliseconds. The experimental study highlights a characteristic profile for the ablation by melting in turbulent flow.

The numerical study presents and validates the ablationFOAM solver dedicated to the prediction of the aero-thermal heating and of the ablation by melting in hypersonic flow. The contribution of the development of the turbulent boundary layer is described. The use of several turbulence models is also discussed.