Christina Ronan

Matthéo, Marin Léa

Jean-Gabriel Paul

Keywords: Instabilité, mélange, turbulence, Rayleigh-Taylor, Faraday, fluides miscibles, Fusion par Confinement Inertiel, DNS, expérimentation, analyse de stabilité linéaire Instabilities, mixing, turbulence, Rayleigh-Taylor, Faraday, miscible fluids, Inertial Confinement Fusion, DNS, experiment, linear stability analysis mean density gradient . . 2.5.2 Wavebreaking detection . . . . . . . . . . . . . .

Cette thèse est consacrée à l'étude de la production de mélange turbulent par des instabilités de Faraday et Rayleigh-Taylor. La première partie est dédiée à l'instabilité de Faraday pour des fluides miscibles. Après avoir mis en évidence le mécanisme sélectionnant le mode de l'onde primaire stationnaire, nous étudions sa transition sous-critique à la turbulence. Nous obtenons, à partir de critères théoriques validés expérimentalement et numériquement, que la transition finale à la turbulence est causée par le déferlement venant du développement d'une instabilité secondaire paramétrique au noeud de la vague. Dans un deuxième temps, l'instabilité de Rayleigh-Taylor est étudiée à l'occasion de la mise en oeuvre d'une nouvelle installation expérimentale. Ce dispositif original permet de stabiliser le fluide lourd grâce à une grille jusqu'à ce qu'il soit mis en contact avec le fluide léger, permettant un meilleur contrôle des conditions initiales. Le but de l'étude est de déterminer, grâce aux expériences et à des simulations numériques, si la dynamique classique de l'instabilité de Rayleigh-Taylor est préservée malgré la présence de la grille. Une majeure partie de l'étude a porté sur la compréhension des mécanismes responsables de l'apparition d'un jet ascendant au moment du contact entre les deux fluides. Le modèle proposé confirme que la courbure des ménisques au niveau des trous de la grille, créés par un saut de pression entre le fluide lourd et l'air, est à l'origine de ce jet. Ce résultat nous oriente donc vers l'élaboration d'un nouveau dispositif permettant de contrôler finement les niveaux de pression.

" "

Long has paled that sunny sky: Echoes fade and memories die: Autumn frosts have slain July.

Still she haunts me, phantomwise, Alice moving under skies Never seen by waking eyes.

LEWIS CAROLL Through the Looking-Glass: A Boat, Beneath a Sunny Sky

List of Symbols

Introduction

Inertial Confinement Fusion (ICF) is a promising method [START_REF] Nuckolls | Laser Compression of Matter to Super-High Densities: Thermonuclear (CTR) Applications[END_REF] and [START_REF] Atzeni | The physics of inertial fusion: beam plasma interaction, hydrodynamics, hot dense matter[END_REF]) to study matter in extreme conditions of temperature and pressure that is currently under development. Recently, a milestone towards fusion has been reached at NIF (National Ignition Facility) with record values of fusion yield achieved (Zylstra et al. [2021]).

As such, it is of great interest to entities such as the Comissariat à l'Energie Atomique et aux Energies Alternatives (CEA).

This method is based on the compression of a spherical target using high intensity laser beams. The target is a capsule usually filled with a mixture of two hydrogen isotopes, the Deuterium D and the Tritium T, and the goal is to achieve fusion between these two atoms by increasing the pressure and temperature.

Two approaches exist: the 'direct drive' approach (upper right of Fig. 1, [START_REF] Betti | Inertial-confinement fusion with lasers[END_REF]) in which the lasers are focused directly on the target and the 'indirect drive' approach (upper left of Fig. 1, [START_REF] Lindl | Development of the indirect-drive approach to inertial confinement fusion and the target physics basis for ignition and gain[END_REF] and [START_REF] Lindl | The physics basis for ignition using indirect-drive targets on the National Ignition Facility[END_REF]) in which the lasers heat the inner walls of a gold cavity containing the target.

In the direct drive approach, the surface of the capsule is ablated due to the rapid heating. This causes an implosion compressing the fuel, and creates a shock wave, increasing further the pressure and temperature of the DT mixture and resulting in a self-sustaining burn.

The other approach, the indirect drive, is the one currently employed at the Laser Megajoule (LMJ) and the National Ignition Facility (NIF). The goal of this approach is to irradiate the target uniformly with X-rays. The X-rays are produced by the gold cavity called Hohlraum when it is heated by the laser beams. This irradiation of the capsule heats its outer surface, causing a high-speed ablation leading to the implosion of the target. This method should produce a much more homogeneous heating of the target than the direct drive approach possessing more energy losses.

The implosion process is the same in both approaches and it can be divided in three stages (Fig. 1): absorption of the energy coming from the laser beams by the capsule surface and ablation of the outer materials; implosion of the capsule, the extreme thermodynamic conditions required for the fusion are reached; creation of a hotspot through ignition of the fuel surrounded by colder and denser DT mixture. The four bottom images describe the implosion process. First, a spherical capsule is prepared at t = 0 with a layer of DT fuel on its inside surface. Second, the capsule surface absorbs energy which ablates the outer materials. Third, the capsule imploses. Lastly, a hotspot of DT is formed surrounded by colder and denser DT fuel. This last image on the bottom right evidences the effects of the hydrodynamics instability. Image taken from [START_REF] Betti | Inertial-confinement fusion with lasers[END_REF].

However, perturbation to the idealised configuration during the compression can completely alter the dynamics of the implosion and prevent the fusion. The ICF is a complicated process involving several very different physical phenomena such as shocks, magnetic field dynamics, hydrodynamic instabilities and complex equations of state. Thus, elementary configurations need to be studied to better understand the importance and impact of these perturbations and prevent them from hindering the fusion process. In hydrodynamics, these perturbations can happen at the interface between the ablator and the fuel. Thus, the elementary configuration of interest is that of an interface between two fluids of different densities under varying acceleration. As the acceleration in the compression process is unknown, several simplified systems can be considered. Among them: a stable configuration under periodic acceleration, the Faraday instability and an unstable configuration under constant acceleration, the Rayleigh-Taylor instability. Indeed, we know that when the perturbations develop at the ablator/fuel interface, they can be amplified by the Rayleigh-Taylor instability which can lead to the mixing of cold (dense) fuel in the hotspot plasma, decreasing the possible fusion yield.

Thus, this thesis is divided in two parts: the first two chapters focus on the transition to turbulence in the Faraday instability while the emphasis is put on the Rayleigh-Taylor instability in the last three chapters.

The first chapter is a brief state of the art on the Faraday instability. After introducing para-Chapter 0. Introduction metric instabilities, a short summary of some studies dealing with the onset of the Faraday instability with immiscible fluids is proposed. Lastly, the time evolution of the instability when miscible fluids are involved is addressed with the notions of turbulent mixing and saturated state.

The second chapter aims at characterising the mechanism leading to turbulence in the case of Faraday waves in miscible fluids. First, it was found that the mode selected by the primary wave at onset is a subcritical one as a result of a mode competition process. Then, two models for the wavebreaking mechanism are developed. They provide a criterion on the steepness of the wave after which the wavebreaking should occur. These models are then confronted to experimental data and a rather good agreement was found between theory and experiments. This shows that the mechanism involved in the wave breaking which then leads to turbulence is a secondary parametric instability coming from the movement of the primary standing wave.

The third chapter is a concise description of the Rayleigh-Taylor instability and the different studies made on the subject. A summary of the various most known Rayleigh-Taylor experiments are given along with a few numerical effects observed such as the difference between two dimensional and three dimensional simulations.

The fourth chapter presents a new type of Rayleigh-Taylor experiment involving the presence of a grid at the interface between the two fluids. The grid maintains the heavier fluid inside a reversed container until the contact with a bottom container filled with lighter fluid. The problematic is: does the grid's presence throughout the experiment impact the dynamics of the instability? In order to answer this question, experimental and numerical data are collected. It is found that the grid geometry greatly influences the shapes of the developing structures; but it is difficult to conclude on whether the long term dynamics is really influenced as classical values for the growth rate are recovered in some experiments.

The last chapter investigates the cause for the apparition of an ascending central jet observed in the majority of the experiments described in the fourth chapter. Numerical simulations are performed in order to test out several theories about the initial conditions that could reproduce this effect. However, as none of these theories give the expected results, the cause for the ascending central jet has to be found elsewhere. This is why we established an aspiration model translating a pressure jump at the interface between the two fluids, evidenced by the curvature of the meniscii inside the holes of the grid, into an imposed velocity field. This model is implemented inside the simulations and the different parameters are studied to reproduce as faithfully as possible the effect observed in the experiments.

Introduction

This introduction aims at describing the Faraday instability and presenting some of the main studies that were made since its discovery, both in the immiscible and miscible frameworks. In the following Chapter 2, the emphasis is put on the transition to turbulence of standing Faraday waves. This subject is inscribed in the continuity of the observations of the destabilisation of standing waves under horizontal vibrations made by [START_REF] Thorpe | On standing internal gravity waves of finite amplitude[END_REF], and more recently of the works of [START_REF] Gréa | What is the final size of turbulent mixing zones driven by the Faraday instability?[END_REF], [START_REF] Briard | Harmonic to subharmonic transition of the Faraday instability in miscible fluids[END_REF] and [START_REF] Briard | The turbulent Faraday instability in miscible fluids[END_REF] regarding the final saturated state of the miscible Faraday instability.

Parametric Instability

The general mechanism causing the excitation of an oscillator by modulating its natural frequency is called 'Parametric Instability'. Indeed, for specific values of control parameters involving a forcing frequency, some initially stable systems can be destabilized. The example of the pendulum is explained hereafter.

Considering a simple pendulum with natural frequency Ω whose suspension point is under periodic acceleration G = G 0 (1 + F cos(ωt)) with a forcing parameter F = a ω ω 2 /G 0 , G 0 the gravitational acceleration, a ω the oscillations amplitude and ω the oscillations frequency, we find that the oscillation of that pendulum is globally amplified when ω/Ω = 2/n, n being any integer [START_REF] Benielli | Excitation and breaking of internal gravity waves by parametric instability[END_REF]). This is shown in Fig. 1.1 where ω = 2Ω meaning that the apparent gravity oscillates by one period when the pendulum oscillates by half a period. Indeed, when the pendulum moves downward, the apparent gravity is higher than average, thus the pendulum is pushed downwards more strongly; when the pendulum moves upward, the apparent gravity is lower than average, the pendulum can thus move upwards more easily. The modulation of the apparent gravity δg is indicated on the figure. The mechanism of parametric instability illustrated by a sketch of a pendulum at successive phases with the corresponding modulation δg of the apparent gravity (the extremal positions of the pendulum are indicated by dotted lines and its velocity by the blue arrow). The apparent gravity is stronger than average when the pendulum moves downward and is weaker than average when the pendulum moves upward. This is possible when the excitation period is half T, the natural period of the pendulum. Figure greatly inspired by the figure 2 of [START_REF] Benielli | Excitation and breaking of internal gravity waves by parametric instability[END_REF].

Chapter 1. State of the art : The Faraday Instability

When the vertical periodic acceleration is strong enough, the effective potential of the system can be changed by the Kapitza effect and a new stable state of equilibrium can arise. This phenomenon is known as the Kapitza's pendulum (Kapitza [1951]) in which a rigid pendulum subjected to vertical vibrations can balance stably in an inverted position (with the dangling mass above the suspension point). In hydrodynamics, the same phenomenon can be found in the formation of frozen waves that happens when a system containing a lighter fluid above a heavier fluid is strongly oscillated horizontally [START_REF] Wolf | Dynamic stabilization of the interchange instability of a liquid-gas interface[END_REF], [START_REF] Lyubimov | Viscosity effect on the longwave instability of a fluid interface subjected to horizontal vibrations[END_REF] and Gréa & Briard [2019]).

These parametric instabilities have often been studied in internal and interfacial waves (Benielli & Sommeria [1998]). Several other examples of parametric instability can be found in physics as it is a common wave-wave interaction, such as the 'electron decay' instability arising from the decay of a strong electromagnetic wave into two electrostatic waves (an electron plasma wave and an ion acoustic wave). Another example of a parametric decay instability involves the sheath-plasma resonance as it can be parametrically excited through a large amplitude RF pump signal at twice the pump frequency [START_REF] Stenzel | Parametric instability of the sheath-plasma resonance[END_REF]). This type of instability can also be associated with other scientific domains such as medicine. Indeed, the association between the Korotkoff sounds (the sounds detected by a stethoscope in the auscultatory blood pressure measurement) and the parametric instability of a fluid filled elastic tube was evidenced.

The Faraday instability

The Faraday instability was first described by Faraday in 1831 [START_REF] Faraday | On a peculiar class of acoustical figures; and on certain forms assumed by groups of particles upon vibrating elastic surfaces[END_REF]) in the context of acoustical figures. It is a parametric instability involving the resonance between the forcing frequency and the natural frequency of an interface (review in [START_REF] Miles | Parametrically forced surface waves[END_REF]).

This instability consists of the superposition of two fluids of different densities in an initial stable configuration with the lighter fluid, ρ l , above the heavier one, ρ h . The interface between the two fluids can be destabilized under periodic vertical accelerations G(t) which changes the effective gravity. A sketch is shown in Fig. 1.2 to illustrate this principle with the destabilized interface as a dotted curve. A specificity of this instability is that it is subharmonic, meaning that the wave formed at the interface oscillates at half the frequency of the periodic vertical oscillations. The interface between the two fluids follows the dynamics of a Mathieu equation:

Immiscible fluids and onset

Mathieu equation (Benjamin & Ursell [1954])

ξ + 2γ(k) ξ + Ω(k) 2 (1 + F cos(ωt))ξ = 0, (

where ξ is the amplitude of the interface deformation, γ is the viscous damping term, F is the forcing parameter and ω is the frequency of the forcing. The frequency Ω is the natural resonant frequency of the interface in the case of immiscible fluids and the inviscid frequency of the diffuse interface in the case of miscible fluids. It is thus important to consider both cases. Both parameters γ and Ω depend on the horizontal wavenumber k.

It is to be noted that equation (1.1), coming from the decoupling of each inviscid interfacial mode, can only be written in the limit of small damping. Otherwise, the problem can be fully analysed using the method described in [START_REF] Kumar | Parametric instability of the interface between two fluids[END_REF].

Immiscible fluids and onset

Numerous experimental and theoretical works focused on the onset of the immiscible Faraday instability.

It was observed that organised structures appeared at the interface between the two fluids and their pattern formation was quite extensively studied [START_REF] Douady | Experimental study of the Faraday instability[END_REF], [START_REF] Convergence Study Fauve | Parametric instability of a liquidvapor interface close to the critical point[END_REF] and [START_REF] Edwards | Patterns and quasi-patterns in the Faraday experiment[END_REF]). For example, [START_REF] Edwards | Patterns and quasi-patterns in the Faraday experiment[END_REF] observed the emergence of a pattern of parallel lines in the case of a single-frequency forcing at high viscosity (see Fig. 1.3a). They also observed a hexagonal pattern in the case of two-frequency forcing (see Fig. 1.3b). They used a mixture of glycerol and water to obtain these results. These hexagonal patterns (see Fig. 1.3b) along with some square patterns were observed experimentally at large frequency in more recent studies [START_REF] Kityk | Spatiotemporal characterization of interfacial Faraday waves by means of a light absorption technique[END_REF]) and compared to Direct Numerical Simulations (DNS) by [START_REF] Périnet | Numerical simulation of Faraday waves[END_REF]. A review of the different patterns found is proposed in [START_REF] Skeldon | Can weakly nonlinear theory explain Faraday wave patterns near onset?[END_REF] along with a comparison between experiments and weakly non-linear analysis, underlying that three-wave interactions play a fundamental role in the stabilisation of certain patterns.

Chapter 1. State of the art : The Faraday Instability

Recent studies evidenced a secondary instability arising close to the transition towards chaos and leading to vibrations of the 2D Faraday pattern [START_REF] Domino | Faraday wave lattice as an elastic metamaterial[END_REF] and Fig. 1.3c). These vibrations, previously observed by [START_REF] Goldman | Lattice Dynamics and Melting of a Nonequilibrium Pattern[END_REF] and [START_REF] Shani | Localized Instability on the Route to Disorder in Faraday Waves[END_REF], are similar to a 2D transverse elastic wave. The authors fully characterized the dispersion relation of the transverse elastic waves and showed that a Faraday pattern presents an effective shear elasticity. The pattern is observed at the primary transition from the flat surface. The mixture used in both cases is composed of 88% glycerol and 12% water with a kinematic viscosity of 1.00 cm 2 .s -1 . Pictures taken from [START_REF] Edwards | Patterns and quasi-patterns in the Faraday experiment[END_REF]. (c) Top view of the oscillating Faraday pattern. Picture taken from [START_REF] Domino | Faraday wave lattice as an elastic metamaterial[END_REF].

Stability diagrams of the Mathieu equation (1.1) for the instability onset are determined using the Floquet theory. It was found that the unstable tongues lift up under the effects of viscosity and surface tension (see [START_REF] Kumar | Parametric instability of the interface between two fluids[END_REF]). Indeed, Fig. 1.4a presents the stability diagram of ideal fluids whereas Fig. 1.4b shows the stability diagram when taking into account the viscosity. Note the displacement of the minima of the tongues towards higher values of k and a when considering the viscosity: higher forcing are needed to destabilize the interface. Taking into account the viscosity also smooths the bottom of the tongues and widens the band of excited wavenumbers.

The stability analysis of a three-fluid system was very recently performed by [START_REF] Ward | Faraday instability in double-interface fluid layers[END_REF] as they investigated, both theoretically and experimentally, the mechanically forced Faraday instability in immiscible two-interface fluid layers systems. It was found that the coupling of the fluid interfaces gave double-tongued stability curves (see Fig. 1.5a) and that adding a third fluid could either enhance or delay the instability. Indeed, a two-fluid system can be destabilized when a third fluid is added and local stabilisation can be achieved by adjusting the fluid properties and layer heights. Moreover, the behavior of the tongues with frequency and viscosity (see Fig. 1.5b) was found to agree with the one seen in two-fluid cases. Thus, it was shown that the theory of [START_REF] Kumar | Parametric instability of the interface between two fluids[END_REF] can be extended to a three-fluid system and that a good estimate of the three-fluid stability curve can be obtained by superimposing the calculations for two two-fluid systems (top-middle and middle-bottom fluids) and then taking the lowest of the thresholds.

This study on the stability of a double interface system was of particular interest to us as we also performed Faraday experiments with a double interface: one miscible interface between Here, a c is the critical excitation amplitude, g is the gravitational acceleration and k is the wavenumber. This figure was taken from [START_REF] Kumar | Parametric instability of the interface between two fluids[END_REF].

(a) (b)

Figure 1.5: (a) Theoretical predictions for critical amplitudes. Double-tongued stability curves are observed for three-fluid Faraday instability systems, where each section of a tongue is dominated by one of the interfaces. The testing frequency is 9 Hz and the properties of the fluids considered are given in Table I of [START_REF] Ward | Faraday instability in double-interface fluid layers[END_REF]. (b) Effect of viscosity on three-fluid instability thresholds. The black lines are obtained using the same fluid properties than (a), while the red lines are obtained by increasing the viscosity of the middle fluid by a factor of 20. Increasing the middle fluid viscosity causes a shift of the threshold towards higher amplitude and higher wavenumber. The testing frequency is 10 Hz. Figures taken from [START_REF] Ward | Faraday instability in double-interface fluid layers[END_REF].

fresh and salt water and a free surface between fresh water and air 1 . The goal of these experiments was to understand the influence of the free surface on the miscible interface by carefully choosing the forcing frequency. However, the presence of jets forming on the free surface and penetrating the stratification below complicated the study and the preliminary works will be Chapter 1. State of the art : The Faraday Instability reported elsewhere than in this manuscript.

Miscible fluids and Turbulence

The miscible Faraday instability, less studied than its immiscible counterpart, was firstly investigated by Zoueshtiagh et al. [2009]. Of particular interest is the possible development of a turbulent mixing layer at the interface between the two fluids. For these vertically excited waves, the vertical periodic acceleration is G(t) = G 0 (1 + F cos(ωt)), thus the control parameters of the instability are the Atwood number A:

Atwood number A = (ρ h -ρ l )/(ρ h + ρ l ), (1.2)
and the acceleration ratio F:

Forcing parameter

F = a ω ω 2 G 0 , (1.3)
where a ω is the amplitude of displacement of the tank containing the fluids and G 0 is the gravitational acceleration.

The development of the instability is depicted in Fig. 1.6. It shows the time evolution of a Faraday experiment between salt and fresh water [START_REF] Briard | The turbulent Faraday instability in miscible fluids[END_REF]). Starting from the initial condition with a clear interface (Fig. 1.6a), the two fluids are destabilized via the periodic vertical acceleration G(t) with a time growing wave amplitude (Fig. 1.6b and Fig. 1.6c). This leads to the breaking of the wave (Fig. 1.6d), turbulent mixing (Fig. 1.6e) and the final saturated state (Fig. 1.6f). The specific mechanism leading to the wavebreaking will be adressed in Chap. 2.

The final saturated state has been investigated in recent studies. It was found in [START_REF] Gréa | What is the final size of turbulent mixing zones driven by the Faraday instability?[END_REF] that the mixing zone width L(t), which grows in size as the instability develops, reaches an asymptotic size L sat :

L sat = 2AG 0 ω 2 (2F + 4). (1.4)
Indeed, parametric resonances excite internal gravity waves, and as long as these waves keep being excited, turbulent kinetic and potential energies are produced by the periodic forcing. However, the mean density gradient decreases when L(t) increases, meaning that this mechanism cannot be sustained and a final state is reached with a saturated mixing layer size given in Eq. (1.4).

This was validated experimentally and numerically in [START_REF] Briard | Harmonic to subharmonic transition of the Faraday instability in miscible fluids[END_REF] and [START_REF] Briard | The turbulent Faraday instability in miscible fluids[END_REF]. In the latter, several experiments and numerical simulations were made to test the dependence of L sat in Atwood number A and forcing frequency ω. A good agreement was found between the experimental data and the theoretical prediction as shown in Fig. 1.7.

Introduction

This chapter is dedicated to the transition to turbulence of the miscible Faraday instability. As already explained in Chap. 1, the Faraday instability is a parametric instability involving a two-fluid system of different densities (the heavier below the lighter) subjected to a vertical periodic acceleration. In the case of this study, the two fluids are miscible ones with a small density contrast as the experiments were made in the hopes of studying the mixing zone generated by this instability. Moreover, this type of Faraday instability was less studied than the non miscible counterpart.

During the experiments, it was noticed that after some oscillations of the interfacial wave, vortices appeared at the node and the wave would then break leading in the end to the mixing of the two fluids. This phenomenon can be seen in Fig. 2.1, where (a) is an experimental observation of this wavebreaking, while (b) comes from Direct Numerical Simulations. These observations were very similar to the ones made by [START_REF] Thorpe | On standing internal gravity waves of finite amplitude[END_REF]. Lateral plungers were used to destabilize standing waves. These waves also consists of miscible fluids with a small density variation, the lighter one being placed above the heavier one. With this experimental setup, he 0.368 m 0.275 m identified several stages for the instability that we were also able to observe:

1. The wave amplitude increases gradually with a symmetrical and sharply defined profile (see Fig. 2.2a).

2. Then, when the wave's slope reaches around 0.2, blurred regions become apparent near the wave nodes while the rest of the profile remains sharp. This stage is called blurring (see Fig. 2.2b).

3. Irregularities, also called the interface roll-up, then arise at the node, coming from the formation of vortices, for wave slopes of about 0.4. Resulting from the mixing in the region of the node, the interface becomes gradually less distinct (see Fig. 2.2c).

4. At last, the two-dimensional motion breaks down and cross waves with the same frequency and in phase with the plungers occurs which results in mixing.

This breaking process has since then been observed for a Faraday wave between miscible fluids by [START_REF] Kalinichenko | Development of a shear instability in nodal zones of a standing internal wave[END_REF]. He agreed with [START_REF] Thorpe | On standing internal gravity waves of finite amplitude[END_REF] on the fact that the mechanism leading to this wavebreaking is a sort of Kelvin-Helmholtz instability. However, it has been found that a subharmonic secondary parametric instability can be involved in the breaking process of internal gravity waves [START_REF] Convergence Study Mcewan | Parametric instability of internal gravity waves[END_REF], [START_REF] Bouruet-Aubertot | Breaking of standing internal gravity waves through two-dimensional instabilities[END_REF], [START_REF] Benielli | Excitation and breaking of internal gravity waves by parametric instability[END_REF], [START_REF] Staquet | Internal gravity waves: From instabilities to turbulence[END_REF], [START_REF] Sutherland | Internal Gravity Waves[END_REF] and [START_REF] Yalim | Parametrically forced stably stratified flow in a threedimensional rectangular container[END_REF]). This leads to the question at the basis of this work: what is the mechanism at the origin of the breaking of the Faraday waves?

This chapter is organized as follows: we give in Sec. 2.1 and 2.2 a description of the experiments and numerical simulations used for this study. In Sec. 2.3, we analyze the characteristics of the primary Faraday wave, emphasizing in particular the mode selection mechanism. Section 2.4 is dedicated to the wavebreaking process with two theoretical approaches proposed and shedding light on the importance of a subharmonic secondary instability. We then detail our methodology in order to measure the wavebreaking amplitudes in Sec. 2.5. Finally, the analysis and discussion of the results in view of the theoretical predictions are provided in Sec. 2.5.3.

Description of the experiments

Description of the experiments

The experiments were carried out at the GazTransport and Technigaz ® (GTT) Motion Analysis and Testing Laboratory. The setup was initially designed to study free-surface flows and was used in order to study the turbulent mixing driven by vertical vibrations in [START_REF] Briard | The turbulent Faraday instability in miscible fluids[END_REF].

The instability studied here is between brine and fresh water, miscible fluids of small density contrast. To contain the fluids, an acrylic tank of inner dimensions 94 x 67 x 11 cm 3 was used.

The tank was first half-filled with fresh water coloured with food dye (Patent Blue V). The heavier salt water was made using commercial Axal ® softening agent consisting of sodium chloride and was injected very slowly into the container through a diffuser located at the bottom. This allowed to keep a sharp interface between the two fluids. A photograph of the tank and both fluids is visible in Fig. 2.3. In order to generate the vertical acceleration needed for this study, the tank was mounted on the platform of an Aquilon model Symétrie ® hexapod owned by GTT. Only vertical sinusoidal oscillations were generated with this hexapod even though it can induce motions along six axes (translation and rotation). The acceleration was limited to 6.87 m.s -2 (corresponding to a forcing parameter, defined in Chap. 1, F = 0.7), the vertical velocity to 1.6 m.s -1 and the vertical amplitude a ω to ± 65 cm.

The platform had an initial vertical position of Z(t = 0) = b, where Z refers to the vertical position of the hexapod, not to be mistaken with z the coordinate in the moving frame. After 20 s of rest, the hexapod moved in a series of n vertical oscillations Z(t) = b cos(ωt) starting at t = 0 s and ending at 2nπ/ω. Thus, at the end of the cycle, the platform was back at its initial position Z = b. It stood there for 20 s before being taken back slowly to Z = 0. Combining a cosine function, a start and end at Z = b ensures that no additional acceleration exists at the Chapter 2. The subcritical transition to turbulence of Faraday waves in miscible fluids beginning and end of the forcing.

Density profiles were measured before and after experiments with a MicroScale Conductivity and Temperature Instrument (MSCTI) from Precision Measurement Engineering, Inc. The probe was fixed above the tank and the hexapod moved the tank in a slow vertical motion in order to get the vertical conductivity and temperature profiles. A temperature correction on the conductivity measurements had to be applied since the salt water was always cooler than the fresh water due to the lights heating the latter during the filling.

A monochromatic Imaging Development Systems video camera was attached to the oscillating platform by a 1.5 m long Elcom ® arm that was mechanically reinforced in order to minimize the vibrations of the camera. That camera was used to capture, at 45 frames per second, images of LED back-lightened white screen as seen through the tank with the use of general purpose red gelatin filter that was put in front of the lens. This red gelatin filter allowed [START_REF] Briard | The turbulent Faraday instability in miscible fluids[END_REF] to operate in the absorption band of the dye and any effect of small persistent vibrations was removed during the post-processing.

In order to convert the image intensity into normalized concentration, the Beer-Lambert absorption law for the blue dye was used. The intensity values covered all 14 bits of the dynamical range of the camera as the result of the optimization of the aperture and exposure time during the experiments. The initial concentration field was used to perform a point by point calibration of the intensity taking into account the background lighting inhomogeneities. Hence, each image was converted into a depth-averaged and normalized concentration field ⟨C⟩ Y . This concentration field was then compared before and after each experiments to the measurements of the density probe and it was found that the change in optical index in salt water induces an under-estimation of the instantaneous mixing zone size L(t) when calculated from the camera images.

Each experiment lasted typically around 2 to 3 hours with the longest part being the filling of the tank. Moreover, the tank was washed between each experiments. [START_REF] Briard | The turbulent Faraday instability in miscible fluids[END_REF] had access to the GTT facility for 4 weeks and made 47 experiments with various parameters investigated and some configurations done several times in order to check the reproducibility of an experiment. However, due to the hexapod limitations, a whole region of the parameter space could not be explored (see Fig. 2.4). The initial mixing zone length of these experiments varied between 0.5 cm and 2 cm, depending on the quality of the tank filling. Of the 47 experiments performed, only those that had a very sharp initial interface were investigated in this work. This corresponds to 18 experiments shown in Tab. A.

Description of the experiments

Series Number

A F ω [rad s -1 ] k n,m [m -1 ] Mode (n, m) δ [cm] EXPa 0.

Direct Numerical Simulation

To complete the experiments, Direct Numerical Simulations (DNS) were performed (see Tab. B) using a code named 'Stratospec' [START_REF] Viciconte | Self-similar regimes of turbulence in weakly coupled plasmas under compression[END_REF] and [START_REF] Viciconte | Sudden diffusion of turbulent mixing layers in weakly coupled plasmas under compression[END_REF]). These simulations solve the Navier-Stokes-Boussinesq equations (2.1) on a triply periodic cubic box of length L box = 2π. The length of the box was converted into the physical dimension of the tank which is W = 94.6 cm.

Equations solved by the numerical simulations

∂U ∂t + (U • ∇)U = -∇P -2ACG(t)x 3 + ν∇ 2 U, ∂C ∂t + (U • ∇)C = D∇ 2 C, ∇ • U = 0, (2.1a) (2.1b) (2.1c)
where x 3 is the upward vertical direction, U is the total velocity field, C is the total concentration field, G(t) = G 0 (1 + F cos(ωt)) is the vertical periodic acceleration, P is the reduced pressure, ν is the kinematic viscosity and D the scalar diffusivity.

These DNS use a classical spectral Fourier collocation method with two-third rule dealiasing. Moreover, they have 1024 3 points and are performed using a pencil-decomposition with 1024 cores. A third-order low-storage strong-stability-preserving Runge-Kutta scheme with implicit viscous terms is used to calculate the time increment.

Periodicity in the vertical inhomogeneous direction needs to be ensured. To do so, an initial condition in which the concentration field changes sharply from its light to heavy value at the vertical boundaries was chosen (see Eqs. (2.2) and (2.3)). In addition to this initial condition on the concentration field, the velocity field is chosen null initially.

Initial condition

• Initial condition on the concentration field:

C(x, t = 0) = 1 2 1 + tanh z -ζ(x) σ + C inh (x), where ζ(x) = ϵ sin(k n,0 x), (2.2)
with

C inh (x) = - 1 2 tanh z + W/2 σ + tanh z -(W/2 -W/N p ) σ , (2.3)
where N p is the number of points in one direction.

• Initial condition on the velocity field:

U = 0. (2.4)
The 'true' initial concentration profile is given by the first term of Eq. (2.2), where the parameter σ is used to define the initial interface thickness δ = 3σ. This parameter is given for each of the simulations carried out in this work in Tab. B. The function ζ(x) indicates the position of (3,0) where the selected mode appears underlined. The parameter r expresses the initial amplitude ratio (r = 0 corresponding to a pure (2, 0) mode). The initial amplitudes ϵ of the interface perturbation, the y-spanwise perturbation amplitude ϵ 1 for DNSf together with the interface thicknesses δ are also detailed. The computation domain is of cubic size with length W = 94.6 cm or 2W for the series labelled with * . All the DNSs have a 1024 3 grid resolution.

Direct Numerical Simulation

Series Number

A F ω [rad s -1 ] Mode (n, m) r ϵ/ϵ 1 [cm] δ = 3σ [
the initial perturbed interface of sinusoidal shape with wavelength k n,0 and amplitude ϵ. In the Chapter 2. The subcritical transition to turbulence of Faraday waves in miscible fluids majority of our simulations ϵ was set at 1.5 cm. This gives a more diffused initial interface with a larger amplitude in the simulations than in the experiments, ensuring grid convergence in the simulations by having at least 20 grid points across the interface layer. The wavenumber k n,0 can be found from the modes given in Tab. B. Indeed, we have k n,0 = π n L box = n 2 , so, for DNSa1 for exemple, k 6,0 = 3.

DNSb * and DNSc * were used to study the effect of mode selection. For these simulations we have:

ζ(x) = ϵ[r cos(k 3,0 x) + (1 -r) cos(k 2,0 x)],
( 2.5) where the parameter r expresses the initial amplitude ratio between the modes (3,0) and (2,0). These simulations were conducted in a computational domain twice the size of the tank 2W in order to allow the development of odd modes that are otherwise forbidden due to the periodic boundary conditions. This means that the interface thickness δ is doubled and the viscosity and diffusion coefficients are multiplied by 4 to ensure grid convergence.

Since the flow remains 2D even after the wavebreaking for DNSa-e, we performed DNSf to study the full transition to turbulence where the interface position is slightly perturbed in the spanwise y-direction:

ζ(x, y) = ϵ sin(k n,0 x) + ϵ 1 f (y), (2.6)
with f (y) being the normalized white noise function and ϵ 1 being the amplitude of the y disturbance amplitude set at 10 -2 ϵ. Other values of ϵ 1 were also tested with a similar result to DNSf. The breaking of the spanwise symmetry invariance can also be produced by the lateral walls in the experiments. Hence, simulations using penalization layers to reproduce these walls were conducted (but are not presented here) and yield similar results to DNSf. The term C inh is added to ensure the vertical periodicity at the boundaries. However, this term gives two unstable regions at the top and the bottom of our domain where a Rayleigh-Taylor instability can be triggered. Hence, a penalization method, described in [START_REF] Briard | The turbulent Faraday instability in miscible fluids[END_REF], is implemented to avoid this phenomenon.

The goal of the penalization method is to mimic solid walls in our pseudo-spectral code, so that the velocity and concentration fields are frozen inside the penalized layers of height L η . Thus, the height of the unpenalized box is H = W -2L η and its depth is still D = W = 94.6 cm as penalization was only applied on layers at the top and bottom of the computational domain to prevent any triggering of a Rayleigh-Taylor instability. We did not try to reproduce the walls of the tank by adding vertical penalization layers as it does not play a role in the wavebreaking process. The height of the penalization layers L η was varied to verify that it did not affect the dynamics of the interface and in all the simulations presented, the Faraday wave's amplitude was always less than half the height of the unpenalized domain.

The penalization method for the velocity field is inspired from Jause-Labert et al.

[2012] and consists in adding a term in the right hand side of Eq. (2.1a) to imitate a porous medium:

Penalised momentum equation (see Jause-Labert et al. [2012])

∂U ∂t + (U • ∇)U = -∇P -2ACG 0 x 3 (1 + F cos(ωt)) + ν∇ 2 U - 1 η u χU, (2.7)
where η u is the velocity penalization parameter and χ is the mask function.

Mode selection mechanism of the Faraday wave

The velocity penalization parameter is set to 1 × 10 -3 s and the penalization mask is defined as a door function D(x), equal to 1 inside the penalization layer and 0 inside the fluid. The mask is filtered in the spectral space using χ(k) = D(k) exp(-C s (k/N p ) 2 ) to avoid Gibbs oscillations with C s the filter parameter. Here, we set N p = 1024 and C s = 256 to preserve the shape and intensity of the mask.

The concentration field is also penalized as in [START_REF] Kadoch | A volume penalization method for incompressible flows and scalar advection-diffusion with moving obstacles[END_REF]. The mask is the same for the two fields U and C but the penalization parameter is different. Thus, a scalar penalization parameter η c is defined as an effective diffusivity inside the solid domain. The advection equation (2.1b) is modified:

Penalised advection equation (see [START_REF] Kadoch | A volume penalization method for incompressible flows and scalar advection-diffusion with moving obstacles[END_REF])

∂C ∂t + ((1 -χ)U • ∇)C = ∇ • ((D(1 -χ) + η c χ)∇C) . (2.8)
In our simulations we chose η c = 2.3 × 10 -9 m 2 .s -1 as it does not impose any constraint on the time step unlike η u and can be chosen very small.

The main differences between the simulations and the experiments are: i) the computational domain in the simulations is of the same width as the width of the tank in the experiments but not of the same height nor depth, the simulations are thus two dimensional in a threedimensional cubic box whereas the experiments were two dimensional because the tank was of very small depth; ii) the simulations are periodic along x so the rigid walls of the experiments in the main direction are not taken into account; iii) the kinematic viscosity in the simulations is set as ν = 2.26 × 10 -6 m 2 .s -1 which is quite close to the viscosity of the fresh water at 20 • C (ν = 1.31 × 10 -6 m 2 .s -1 ), however, the limited spatial resolution forced us to choose D = ν for the diffusivity whereas in the experiments the Schmidt number, defined as Sc = ν/D, was around 700. The effect of a too large D can mainly be seen on the early time evolution of the mixing zone size L(t) and hindered us from exploring small F configurations as the instability would take a very long time to be triggered, time during which the interface would diffuse considerably numerically but would remain flat in the experiments.

Mode selection mechanism of the Faraday wave

Before studying any secondary instability at the origin of the wavebreaking mechanism, we first focus on the primary wave characteristics and more particularly on figuring out which mechanism selects its dominant wavelength.

Linear theory

As explained in Chapter 1, the amplitudes ξ k for the interface modes of wavenumber k are governed by a Mathieu equation:

ξk + 2γ(k) ξk + Ω 2 (k)(1 + F cos(ωt))ξ k = 0, (2.9)
where Ω is the inviscid frequency of the diffuse interface. The damping term γ will be described in Sec. 2.3.1.1.

Chapter 2. The subcritical transition to turbulence of Faraday waves in miscible fluids

In the deep water approximation, the inviscid frequency Ω(k) is a growing function of k and can be determined for a given vertical density profile. In this work, we considered a piecewise linear profile (see for instance [START_REF] Briard | The turbulent Faraday instability in miscible fluids[END_REF]):

Inviscid frequency

Ω(k) = AG 0 k 1 + kδ/2 1/2 .
(2.10)

The classical dispersion relation for an interface within the deep water approximation

Ω(k) =
√ AG 0 k is recovered for small wavenumbers (kδ ≪ 1) whereas the interface mode reduces to the local buoyancy or Brunt-Väisälä frequency at the interface, Ω = 2AG 0 /δ, in the large wavenumber limit (kδ ≫ 1).

The viscous damping term

The viscous dissipation term γ(k), corresponding to small interfacial mode damping, can have different origins:

1. it can come from the bulk flow for a sharp interface giving it the form:

Bulk damping (see [START_REF] Lamb | Hydrodynamics[END_REF] and [START_REF] Landau | [END_REF])

γ b (k) = 2νk 2 ;
(2.11)

2. due to velocity gradients, it can also come from the thin layer separating the two fluids.

Considering a piecewise linear profile, [START_REF] Briard | The turbulent Faraday instability in miscible fluids[END_REF] found the following expression for the damping:

Interfacial layer damping (see Briard et al. [2020])

γ δ (k) = (AG 0 νk 2 )/(Ω 2 δ) ∼ νk/δ for kδ ≪ 1;

(2.12)

3. the last contribution to the damping to take into consideration is the one coming from the boundary layers at the various walls in the experiments. Keulegan [1959], [START_REF] Miles | Surface-wave damping in closed basins[END_REF] derived an expression for the damping of a free surface wave in a rectangular basin due to the laminar boundary layers which can be generalized to our problem as [START_REF] Thorpe | On standing internal gravity waves of finite amplitude[END_REF] did. The calculations to derive the expression of that contribution are detailed below.

First, the boundary layer width can be expressed as δ w = (2ν/Ω) 1/2 giving δ w ∼ 1 -2 mm when considering the parameters of the experiments (ν = 1.31 × 10 -6 m 2 .s -1 , Ω = √ AG 0 k s -1 , A = 0.03 and k = π/W m -1 ). Note that this boundary layer width is much smaller than the characteristic wavelengths of the instability and the size of the tank. Then, taking a rectangular basin of height H, width W and depth D (see Fig. 2.5) and considering a wave of small amplitude a, horizontal wavenumber k and oscillating at the frequency Ω, we have derived the contribution of each wall to the damping. The rest interface position is located at the middle of the tank z = 0. We will consider two cases: one where the top boundary is a wall and one where it is a free surface.

Indeed, from [START_REF] Landau | [END_REF], we know that the damping coefficient is written γ = ⟨ Ėmech ⟩/2⟨E mech ⟩ where E mech is the mechanical energy. We consider that the amplitude of the wave a is very small (a ≪ 1) so that the wave moves linearly. The velocity potential ϕ for the primary mode can be expressed as (Keulegan [1959]):

   ϕ = - aG 0 Ω cosh kH/2 cosh k(z + H/2) cos kx cos Ωt, for z < 0, ϕ = - aG 0 Ω cosh kH/2 cosh k(z -H/2) cos kx cos Ωt, for z > 0.
(2.13) Our system being symmetric, all of the computations performed on only the lower part of the tank with E tot = 2E lower . From this velocity potential, we can calculate the particle velocities for the lower part as u = -∂ x ϕ and w = -∂ z ϕ.

   u = -G 0 ak Ω cosh k(z+H/2) cosh kH/2 sin kx cos Ωt, w = G 0 ak Ω sinh k(z+H/2)
cosh kH/2 cos kx cos Ωt.

(2.14)

We write the mean mechanical energy for the lower part of the tank as being twice the mean kinetic energy for that lower part which gives:

⟨E mech,low ⟩ = ρ V ⟨v 2 ⟩dV = ρ V u 2 + w 2 dV.
(2.15)

Chapter 2. The subcritical transition to turbulence of Faraday waves in miscible fluids

Taking the mean over a period of time of cos 2 (Ωt) brings a factor of a half, so we get:

⟨E mech,low ⟩ = 1 2 ρ G 0 ak Ω cosh kH/2 2 V cosh 2 k(z + H/2) sin 2 kx + sinh 2 k(z + H/2) cos 2 kx dV = 1 2 ρ G 0 ak Ω cosh kH/2 2 W 0 D 0 0 -H/2 sin 2 kx + sinh 2 k(z + H/2) dzdydx = 1 2 ρDW G 0 ak Ω cosh kH/2
2 sinh kH 4k .

(2.16)

In order to get the total mechanical energy, we have to sum the mechanical energy of the lower part of the tank and the one of the upper part of the tank. Thanks to the symmetry in our system, we have:

⟨E mech,tot ⟩ = 2⟨E mech,low ⟩ = ρDW G 0 ak Ω cosh kH/2 2 sinh kH 4k .
(2.17)

Now, we need to evaluate the damping coefficients for each wall of the tank.

Contributions of the vertical walls y = 0 and y = D: Considering the wall at y = 0 for the lower part of the tank (z ≤ 0), we have the following expression for the fluid velocities in the viscous layer:

   u = G 0 ak Ω cosh k(z+H/2)
cosh kH/2 sin kx e -y/δ w cos(Ωty/δ w ) -cos Ωt , w = -G 0 ak Ω sinh k(z+H/2) cosh kH/2 cos kx e -y/δ w cos(Ωty/δ w ) -cos Ωt , (2.18) giving a null velocity at the wall and the potential solution u 0 = ∇ϕ away from it.

Taking the square of the y-derivative, we get:

     (∂ y u) 2 = G 0 ak Ω 2 cosh 2 k(z+H)/2 cosh 2 kH/2 sin 2 kx 1 δ 2 w e -2y/δ w [1 + sin 2(y/δ w -Ωt)] , (∂ y w) 2 = G 0 ak Ω 2 sinh 2 k(z+H) cosh 2 kH/2 cos 2 kx 1 δ 2 w e -2y/δ w [1 + sin 2(y/δ w -Ωt)] . (2.19)
The mean kinetic energy dissipated on this wall can be written: .20) Adding (∂ y u) 2 and (∂ y w) 2 and integrating over a complete cycle of period 2π/Ω, we have:

⟨ Ė1 ⟩ = ρν Ω 2π 0 -H/2 W 0 ∞ 0 2π/Ω 0 (∂ y u) 2 + (∂ y w) 2 dtdydxdz. ( 2
2π/Ω 0 (∂ y u) 2 + (∂ y w) 2 dt = 2π Ω (G 0 ak) 2 Ω 2 e -2y/δ w δ 2 w cosh 2 kH/2 sin 2 kx + sinh 2 k(z + H/2) . (2.21)

Mode selection mechanism of the Faraday wave

Hence, at leading order in 1/δ w ,

⟨ Ė1 ⟩ = Ω 2π 2πρν(G 0 ak) 2 Ω 3 1 δ 2 w cosh 2 kH/2 ∞ 0 e -2y/δ w dy W 0 0 -H/2 sin 2 kx + sinh 2 k(z + H/2) dzdx = Ω 2π (G 0 ak) 2 Ω 2 ρνδ w 2δ 2 w Ω 2πW (H/4 + sinh kH/4k -H/4) cosh 2 kH/2 = ρνW 2δ w aG 0 k Ω cosh kH/2
2 sinh kH 4k .

(2.22)

The contribution from the wall at y = D is the same and considering the upper part of the rectangular basin, we get ⟨ Ė1,tot ⟩ = 4⟨ Ė1 ⟩.

From this, we can evaluate the contributions from these walls to the damping:

γ w 1 = ⟨ Ė1,tot ⟩ 2⟨E mech,tot ⟩ = ν δ w D . (2.23)
Contributions of the vertical walls x = 0 and x = W: Now considering the vertical wall x = 0 in the bottom part of the tank, the fluid velocities in the viscous layer are:

   u = 0, w = -G 0 ak Ω sinh k(z+H)
cosh kH/2 e -x/δ w cos(Ωtx/δ w ) -cos Ωt .

(2.24)

Hence,

(∂ x w) 2 = G 0 ak Ω 2 sinh 2 k(z + H/2) cosh 2 kH/2 1 δ 2 w e -2x/δ w [1 + sin 2(x/δ w -Ωt)] . (2.25)
The amount of energy dissipated on this wall is:

⟨ Ė2 ⟩ = Ω 2π 2πρν(G 0 ak) 2 Ω 3 1 δ 2 w cosh 2 kH/2 ∞ 0 e -2x/δ w dx D 0 dy 0 -H/2 sinh 2 k(z + H/2)dz = (G 0 ak) 2 ρνδ w 2δ 2 w Ω 2 D (sinh kH/2k -H/2) 2 cosh 2 kH/2 = ρνD 2δ w aG 0 k Ω cosh kH/2 2 sinh kH 4k - H 4 .
(2.26)

Again, to have the full contributions of the walls at x = 0 and x = W, we have ⟨ Ė2,tot ⟩ = 4⟨ Ė2 ⟩. Thus:

γ w 2 = ⟨ Ė2,tot ⟩ 2⟨E mech,tot ⟩ = ν δ w D 1 - kH sinh kH . ( 2
   u = -G 0 ak Ω 1 cosh kH/2
sin kx e -(z+H/2)/δ w cos(Ωt -(z + H/2)/δ w ) -cos Ωt , w = 0.

(2.28) So,

(∂ z u) 2 = G 0 ak Ω 2 sin 2 kx cosh 2 kH/2 1 δ 2 w e -2(z+H/2)/δ w [1 + sin(2((z + H/2)/δ w -Ωt))] . (2.29)
The corresponding loss of energy is given by:

⟨ Ė3 ⟩ = Ω 2π 2πρν(G 0 ak) 2 Ω 3 1 δ 2 w cosh 2 kH ∞ -H/2 e -2(z+H/2)/δ w dz D 0 dy W 0 sin 2 kxdx = (G 0 ak) 2 Ω 2 ρνδ w 2δ 2 w W L 2 cosh 2 kH/2 = ρνDW 4δ w aG 0 k Ω cosh kH/2 2 .
(2.30)

In the case where there is no free surface, we get: .31) whereas in the case of a free surface, we have: .32) Total contribution : The total contribution of the walls to the damping is γ w = γ w 1 + γ w 2 + γ w 3 : .33) for the case without free surface and .34) for the case with a free surface. However, we will only consider the dominant contribution of the lateral walls to the damping, which can be expressed as: .35) where D = 11 cm.

γ w 3 = 2⟨ Ė3 ⟩ 2⟨E mech,tot ⟩ = νk δ w sinh kH , ( 2
γ w 3 = ⟨ Ė3 ⟩ 2⟨E mech,tot ⟩ = νk 2δ w sinh kH . ( 2
γ w = ν δ w 2 D + k sinh(kH) 1 - H D , ( 2
γ w = ν δ w 2 D + k sinh(kH) 1 2 - H D , ( 2

Mode selection mechanism of the Faraday wave

Boundary layers damping

γ w ∼ ν Dδ w = √ νΩ √ 2D , ( 2
The numerical values of the damping coefficients coming from the bulk, the interfacial layer between the two fluids and the boundary layers at the walls are given in Tab. C. Even if we account for the viscosity contrast between fresh and salt water, the values do not vary more than a few percent. This indicates clearly that the dissipation occurs essentially in the viscous layers at the walls since γ w is the largest of all the contributions. Indeed, γ w scales like ν 1/2 whereas the bulk dissipation γ b and the interfacial dissipation γ δ are linear in ν. However, for larger wavenumbers, the contribution from the interface layer becomes the most important contribution to the damping as can be seen in Fig. 2.6.

These aspects have been extensively studied by [START_REF] Bechhoefer | An experimental study of the onset of parametrically pumped surface waves in viscous fluids[END_REF]. The authors suggested using fluids with high viscosities to better control the dissipation in experiments performed in order to study the instability threshold. As our work was focused on the wavebreak-

Mode (2,0) (3,0) (4,0) (5,0) (6,0) (7,0) (8,0) γ w = 9.75 ×10 -3 1.08 ×10 -2 1.15 ×10 -2 1.22 ×10 -2 1.28 ×10 -2 1.32 ×10 -2 1.36 ×10 -2 γ b = 1.16 ×10 -4 2.60 ×10 -4 4.62 ×10 -4 7.22 ×10 -4 1.04 ×10 -3 1.42 ×10 -3 1.85 ×10 -3 γ δ =
1.74 ×10 -3 2.61 ×10 -3 3.48 ×10 -3 4.35 ×10 -3 5.22 ×10 -3 6.09 ×10 -3 6.96 ×10 -3

Table C: Values for the damping coefficients γ w , γ b and γ δ in s -1 and evaluated for the largest wavelengths developing in the experiments. We assume here that the Atwood number is A = 0.03 and the thickness of the interfacial layer is δ = 0.5 cm. Here, the top boundary is taken as a wall to evaluate γ w (the values would be nearly the same for a free surface). Chapter 2. The subcritical transition to turbulence of Faraday waves in miscible fluids ing mechanism, we favored the use of low viscosity fluids.

Stability diagram

The stability diagram representing the first subharmonic tongue is shown in Fig. 2.7. This diagram is plotted for different large scale modes of the tank using the damping coefficients from the interfacial layer γ δ and from the boundary layers at the walls γ w . The neutral curves of the Mathieu equation (2.9) are computed using the method proposed by [START_REF] Kumar | Parametric instability of the interface between two fluids[END_REF], and previously used in [START_REF] Briard | The turbulent Faraday instability in miscible fluids[END_REF], assuming different Atwood number A values and an initial interface thickness δ = 1 cm. For a given wavenumber k, the minimum forcing F th needed to destabilize the interface occurs at a frequency corresponding to the first subharmonic resonance Ω(k) = ω/2. The classical asymptotic theory of the Mathieu equation gives F th = 8γ/ω, in the limit of small damping (see for instance [START_REF] Rajchenbach | Faraday waves: their dispersion relation, nature of bifurcation and wavenumber selection revisited[END_REF]). The minimum forcing varies very weakly for the different modes in the figure. Indeed, the contribution due to the damping from the viscous layers at the walls γ w scales like ω 1/2 , leading to a decrease of F th at larger ω which is compensated by the contribution from the damping at the interface γ δ scaling like ω 2 . The parameters taken in the experiments, shown in Fig. 2.7, are located in the unstable regions well above the viscous thresholds determined by the linear Floquet theory. Consequently, for most of the experiments, two or more modes can be simultaneously subharmonically unstable.

Linear or nonlinear mode selection?

In this section, we investigate the mechanisms leading to the mode selection of the primary wave. As discussed in the previous section, several modes can be linearly unstable, and play a role in the interface dynamics, due to the large forcing acceleration parameter F. Surprisingly, a single mode, corresponding nearly always to the smallest unstable wavelength, emerges from this process. There is a clear tendency to favor the modes pertaining to the right unstable tongues in Fig. 2.7, the mode reported in Tab. A being also indicated by the color of the symbol in the stability diagram. Moreover, the selection mechanism does not apparently discriminate between the even or the odd modes of the tank as both can be observed in the experiments.

We only considered the modes in the first subharmonic bands because we expected that the modes with the largest growth rates would be the ones selected and the higher resonance regions have much smaller amplification rates. For a given mode, the Floquet theory shows that the maximum amplification occurs for parameters close to the subharmonic resonance frequency located at the center of the instability tongue. However, the results in Fig. 2.7 reveal that in many cases, the selected mode does not have the largest linear growth rate. Moreover, some of the observed modes are hardly unstable and have very small growth rate from the linear theory (see for instance the EXPa1) and that is even if we take into account some experimental uncertainties in terms of Atwood number A (±0.001) or initial interface thickness δ (±0.5 cm) which would only slightly modify instability tongues of Fig. 2.7. Indeed, a larger interface thickness would decrease the natural frequencies Ω leading to a slight left-shift of the instability tongues. The movements of the hexapod were well controlled and remained sinusoidal so it is very unlikely that spurious forcing frequencies would appear in the system and change the linear stability of the problem.

The initial perturbation of the interface can also play a role in the mode selection mechanism. Indeed, a large initial amplitude on a given mode would explain why it appears even if it has not the largest growth rate during the linear phase. This would suggest that an initial condition at small scales is at work in the experiments although we have not observed such disturbance or could not identify a source able to generate it. However, this would not explain the appearance of linearly stable modes.

However, nonlinearities have the ability to select modes and generate transient chaotic regimes for Faraday experiments with immiscible fluids (see for instance [START_REF] Ciliberto | Chaotic mode competition in parametrically forced surface waves[END_REF]). Using weakly nonlinear approaches, [START_REF] Meron | Theory of chaos in surface waves: The reduction from hydrodynamics to few-dimensional dynamics[END_REF] and [START_REF] Meron | Parametric excitation of multimode dissipative systems[END_REF] have already detailed how the mode suppression phenomenon can occur. Considering two modes close to the first subharmonic resonance, the sign of the nonlinear cubic coupling terms in each amplitude equation is determined by the detuning parameter ∆ = (Ω 2 /ω 2 ) -(1/4). Here, the detuning parameter for a given mode expresses the departure of the forcing frequency from the first subharmonic resonance. Thus, if ∆ = 0 there is a perfect resonance between the forcing frequency ω and the subharmonic resonance frequency Ω. Therefore, when two modes compete, their respective detuning parameters ∆ generally have opposite signs because the forcing frequency ω lies between the two subharmonic resonance frequencies (see Fig. 2.7). One mode can develop by pumping energy from the other one, even if it is linearly stable, which explains the mode suppression. Indeed, as the wave amplitude grows, the frequency of the wave Ω tends to diminish [START_REF] Thorpe | On standing internal gravity waves of finite amplitude[END_REF]), hence the detuning parameter gets smaller [START_REF] Fauve | Pattern forming instabilities[END_REF]), left-shifting the instability bands and making the modes become more or less resonant depending on the value of their initial detuning parameter. This mode suppression process is shown in Fig. 2.8 with the shifting of the tongues either towards resonance (Fig. 2.8a) for ∆ ≥ 0 initially, or further away from resonance (Fig. 2.8a) for ∆ < 0 initially. The theory shows that the vanishing mode is supercritical (∆ < 0), as the nonlinear coupling damps the instability, while the dominant one is subcritical (∆ ≥ 0), as it is reinforced by the nonlinearities. Thus this process favors the subcritical modes at smaller wavelength.

Numerical analysis of the mode competition

At this stage, the mode suppression due to a nonlinear coupling between competing modes can explain the mode selection demonstrated in Fig. 2.7. Furthermore, the mode amplitude is no longer negligible compared to its wavelength when the selection process is at work, suggesting that the nonlinear effects are an important ingredient to account for. To further assess this hypothesis, direct numerical simulations were performed with well-characterized initial conditions.

Two series of DNS have been performed using 1024 3 grid points and an Atwood number of A = 0.03 (Series DNSb and c in Tab. B). The frequency ω and the forcing F taken in the simulations are represented in the phase diagram of Fig. 2.9. It is important to stress here that the phase diagram does not take into account the wall damping as simulations are performed in a triply periodic box. These two series of DNS start from the same location in the phase diagram (point A) with F = 0.8 and ω = 3.07 (or equivalently ω/ AG 0 k 2,0 = 2.2). This corresponds to parameters with the two unstable modes having nearly the same exponential growth rate as ∼ e µωt . Indeed, the Floquet exponent µ (Sec. 2.4.1.2) takes the value 0.09 for mode (3,0) and 0.07 for mode (2,0). We fix the forcing parameters in the DNSc series and vary the amplitude ratio r. In the other series, the relative amplitude is set at r = 0.5 and we decrease the forcing frequency ω in order to explore more deeply the (2,0) subharmonic instability tongue. The simulations are stopped when the wavebreaking occurs and we report which of the (2,0) or (3,0) mode prevails at this moment in Fig. 2.9. This procedure is performed both visually and by computing the Fourier modes of the interface.

The simulations clearly demonstrate mode suppression to the benefit of the modes with the smallest wavelength. The results, reported in Fig. 2.9 show the dominance of mode (3,0) The colored areas correspond to the subharmonic instability tongues for the modes (2,0) and (3,0) accounting for viscosity and a diffused interface δ = 1.8 cm. The symbols' color indicates which mode emerges from the simulation, the mixed colors express that both modes can be observed. The two series DNSb and DNSc (see Tab. B) are presented here starting from point A. In the DNSb group, the initial amplitude ratio r between modes (2,0) and (3,0) is set at r = 0.5 and the forcing frequency is decreased. In the DNSc group, the frequency is fixed and r is varied. The point corresponding to DNSe is also placed.

if it starts with a small initial amplitude (the transition occurs at r = 0.1) or if it is in a region where it is linearly stable (for small ω). The phenomenon can be observed in more detail on the snapshots extracted from the two series in Fig. 2.10, where the mode (3,0) emerges from cases with initial r = 0.25, 0.5 or with the frequencies ω = 2.6, 2.8. Indeed, in the last row of Fig. 2.10(b), one can observe that mode (2,0) prevails only for ω = 2.4 (ω/ AG 0 k 2,0 ) and for larger ω, say ω = 2.6 (ω/ AG 0 k 2,0 ) mode (3,0) is visible despite being linearly stable, as can be seen in Fig. 2.11, where the Fourier mode of the interface are shown for the last row of Fig. 2.10(b). As a consequence, the mode competition greatly enhances the sensitivity to initial conditions in the experiment. As importantly, this process breaks the symmetry of the primary wave as can be seen in both the experiments (Fig. 2.1) and the simulations (Fig. 2.10).

We now address some specifities of the mode suppression in the case of our Faraday experiments with miscible fluids. We have not observed oscillations between two specific modes as in [START_REF] Ciliberto | Pattern competition leads to chaos[END_REF] and [START_REF] Ciliberto | Chaotic mode competition in parametrically forced surface waves[END_REF] or similarly in [START_REF] Yalim | Parametrically forced stably stratified cavity flow: complicated nonlinear dynamics near the onset of instability[END_REF] in the context of a stable stratification. This is notably because, at large forcing parameter, the interface irreversibly grows allowing continuously new modes to be destabilized. The modes can change in our experiment as already reported in [START_REF] Briard | The turbulent Faraday instability in miscible fluids[END_REF]. But it always corresponds to a one way transition from large to small wavelength for interface modes. The more complex transitions evidenced in figure 14 of [START_REF] Briard | The turbulent Faraday instability in miscible fluids[END_REF], for instance, refer to modes pertaining We put the slices of width W of the concentration field at four instants starting from the initial interface at t = 0 and ending when wavebreaking occurs; pure fluids remain in white while mixed fluid with C = 0.5 ± 0.15 is in black. not to the same instability band or being of different nature. The irreversible mixing produced by the rapid breaking of the primary waves also explains this aspect.

Another difference with past Faraday immiscible experiments conducted in a shallow basin is that in our case the dominant waves correspond to those with the smallest wavelength as already discussed. This point clearly agrees with the nonlinear theory of mode suppression. It can be shown that within the deep water approximation, the subcritical modes are indeed those with small wavelength (see [START_REF] Rajchenbach | Faraday waves: their dispersion relation, nature of bifurcation and wavenumber selection revisited[END_REF] for details).

In this part, it has been evidenced that the mode selection of the primary wave may result from a complex nonlinear mode competition process. When this is the case, the subcritical mode is eventually selected. In the following, we explain the breakdown of the Faraday waves.

Modelling the breakdown of Faraday waves

Modelling the breakdown of Faraday waves

The goal of this study is to evaluate the critical wave steepness at which the wavebreaking occurs. To do so, two heuristic models are presented in this section (see Fig. 2.12). These models are dedicated to the breakdown of Faraday waves initiating the transition to turbulent mixing and allow us to explore various frameworks for the breakdown and to disentangle the inner mechanisms responsible for the instability.

Both models, although relying on different assumptions, suggest that the breakdown results from a subharmonic secondary instability at small scales. Therefore, one key ingredient in these approaches is to account for the unsteadiness of the primary wave. This aspect differs from secondary instability analysis relying on a frozen base flow used for instance in the context of the Kelvin-Helmholtz instability [START_REF] Salehipour | Turbulent diapycnal mixing in stratified shear flows: the influence of Prandtl number on mixing efficiency and transition at high Reynolds number[END_REF]).

For the first approach, hereafter referred to as global, we seek to identify the conditions for which a small disturbance can develop around the mean profile which characterizes a Faraday wave of given amplitude (see Fig. 2.12(b)). Indeed, this theory is based on the fact that monotonicity of the horizontally averaged density profiles is changed due to the breakdown of the Faraday waves. In contrast, the second model, later called local, relies on the local analysis of small perturbations at the node of the primary wave (see Fig. 2.12(c)). .12: Frameworks applied to model the wavebreaking of the primary wave (a) and detailed in Sec. 2.4. For the global approach (b), we consider the stability of the horizontally averaged concentration profiles. For the local approach (c), we study the development of small perturbations of wavenumber k wb at the node of the primary wave.

Chapter 2. The subcritical transition to turbulence of Faraday waves in miscible fluids

The global approach 2.4.1.1 A simple model equation

In order to derive a simple model for the breakdown of Faraday waves, the system of equations (2.1) is averaged. This gives equations for the mean flow and its fluctuations by decomposing the concentration and velocity fields as C = C + c and U = Ū + u. The mean quantities * are obtained by averaging along the horizontal x and y directions and c and u are the fluctuating part. Due to the symmetries and the incompressibility condition, it can be shown that the mean velocity is zero ( Ū = 0). Meanwhile, the mean vertical concentration profile C(z, t) evolves principally due to the vertical buoyancy flux wc as: .36) In this global approach, the primary Faraday wave is thus embedded in the mean vertical concentration profile C(z, t) while having fluctuation components satisfying (2.36).

∂ t C + ∂ z wc = κ∂ 2 zz C. ( 2
Further simplifications are now made in order to obtain a simpler model. We seek a 'rapid' secondary instability occurring at small scales and located at z = 0. In this context, the rapid acceleration theory, initially developed by [START_REF] Hunt | Rapid distortion theory and the problems of turbulence[END_REF] and applied to turbulent mixing layer in [START_REF] Gréa | The rapid acceleration model and the growth rate of a turbulent mixing zone induced by Rayleigh-Taylor instability[END_REF], provides a convenient framework for expressing the dynamics of small scale disturbances. We thus discard all the dissipative and the non linear terms, except for the coupling between the fluctuations and the mean concentration field. The small scale disturbance only sees a uniform mean concentration gradient at z = 0 determined by the mixing width:

Inverse mean concentration gradient

L = -1/∂ z C(0, t).
(2.37)

Within this quasi-homogeneous approximation, the small scale fluctuating quantities and pressure p are determined by:

∂ t u = -∇p -2AG(t)cx 3 , ∂ t c = w L(t) , ∇ • u = 0. (2.38a) (2.38b) (2.38c)
The classical equations for an internal gravity wave driven, here, by a time varying acceleration and mean density gradient can be recognized in the system of equations (2.38). This mean gradient can be evaluated from the mixing zone width (see Eq. (2.39)).

Modelling the breakdown of Faraday waves

Integral mixing zone width [START_REF] Andrews | A simple experiment to investigate two-dimensional mixing by Rayleigh-Taylor instability[END_REF])

L int = 6 +∞ -∞ C(1 -C)dz.
(2.39)

It was introduced by [START_REF] Andrews | A simple experiment to investigate two-dimensional mixing by Rayleigh-Taylor instability[END_REF] and previously used in [START_REF] Briard | The turbulent Faraday instability in miscible fluids[END_REF] for the fully turbulent regime in the case where this mean density gradient is uniform across the layer. However, this property is lost when the flow takes the form of a single laminar wave. In this case, the inverse concentration gradient is maximum at z = 0. Different ways to evaluate this gradient will be discussed in a next section.

These waves depend on their orientations but for this heuristic model we focus only on waves with a wavevector in the horizontal plane. Differently oriented modes are thought to be less relevant in the secondary instability, partly because they are less likely to modify the mean concentration profile. Indeed, the feedback of the fluctuations to the mean concentration profile is controlled by the vertical buoyancy flux term wc which is weaker for vertically oriented modes. Eliminating w in Eqs. (2.38), we obtain the following Mathieu-like equation (a detailed derivation can be found in [START_REF] Gréa | What is the final size of turbulent mixing zones driven by the Faraday instability?[END_REF]):

c + L L ċ + 2AG(t) L c = 0. (2.40)
The concentration fluctuations c are therefore fully determined by L expressing the amplitude of the primary Faraday wave. The condition on L, for which the perturbations can develop, needs to be determined. Indeed, the rise of the perturbation foreshadows the breaking process of the primary wave and the onset of turbulence. More precisely, Eq. (2.40) exhibits the buoyancy frequency defined as Ω B = (2AG 0 /L) 1/2 and a damping term L/L expressing the variations of Ω B as the mixing zone width L evolves. As will be seen in Sec. 2.4.2, the Ω B frequency is relevant for the secondary instability because the shear at the nodes of the primary wave is directly driven by the wave amplitude.

We now detail the implications of this model equation regarding the wavebreaking phenomenology as observed in Fig. 2.1.

The subcritical nature and the criterion for the wavebreaking

The stability of the model equation (2.40) allows to evaluate the saturation criterion of the instability. This gives when the subharmonic instability stops or, equivalently, when the inner frequencies of the layer and the forcing frequency are no longer in resonance. Hence, this has been extensively discussed in [START_REF] Gréa | What is the final size of turbulent mixing zones driven by the Faraday instability?[END_REF] and [START_REF] Briard | The turbulent Faraday instability in miscible fluids[END_REF] in order to predict the final widths of the turbulent mixing zones. But the authors did not seek to account for the unsteadiness of L. However, as we interpret the wavebreaking as the development of a secondary subharmonic instability at small scales, we aim at finding the onset of this instability. Indeed, it is our belief that a small disturbance, characterized by the buoyancy frequency Ω B , will become parametrically unstable as the result of both the forcing and the primary wave oscillations. Thus, the enlargement of the primary wave amplitude determines not only the instability threshold but also the later amplification of the secondary instability growth rate. This explains the subcritical nature of this secondary instability and why it develops rapidly at the Chapter 2. The subcritical transition to turbulence of Faraday waves in miscible fluids interface. This peculiarity is indeed well-known for non linear Mathieu systems such as (2.40) as detailed in [START_REF] Soliman | Indeterminate sub-critical bifurcations in parametric resonance[END_REF] or [START_REF] Fauve | Pattern forming instabilities[END_REF].

In order to derive an analytic criterion for the wavebreaking, we consider the inverse mean concentration gradient L having the following simple form L(t) = L 0 (1 + β cos(ωt)). The length L is thus expected to be proportional to the amplitude ξ p of the primary Faraday wave, L(t) ∼ |ξ p | in the laminar phase. Here, the proportionality coefficient depends on the shape of the nonlinear primary wave. Also, the parameter β expresses the relative amplitude of the Faraday wave oscillations while L 0 is the mean over one oscillation period. This expression does not account for the primary mode growth, which is assumed small over an oscillation period. It also expresses that for a subharmonic instability, L oscillates at the frequency ω while ξ p oscillates at the frequency ω/2. However, the higher temporal harmonics of L(t), or ξ p for the primary wave, are discarded. In addition, it is important to note that the primary Faraday mode amplitude is in phase with the acceleration G(t).

In this context, we use Floquet analysis to find the secondary subharmonic instability onset. Indeed, with L(t) = L 0 (1 + β cos(ωt)) injected in Eq. ( 2.40), we have: (2.41) where τ = ωt and Ω 2 0,B = 2AG 0 /L 0 . Applying the Floquet theorem, the solutions for a linear equation with time periodic coefficients (Eq. (2.41)) are known and take the form:

(1 + β cos(τ)) c -β sin(τ) ċ + Ω 2 0,B ω 2 (1 + F cos(τ))c = 0,
Floquet Theorem on c c = N=+∞ ∑ n=-∞
Y n e (µ+i(n+α))τ , (2.42) where µ is the Floquet exponent characterizing the growth rate of the instability. The most unstable modes generally correspond to α = 0 or α = 1/2 which here represent the harmonic or subharmonic modes respectively.

Taking µ = 0 allows to determine the neutral branches of the instability. Inserting the solution (2.42) in Eq. ( 2.41), one gets:

N=+∞ ∑ n=-∞ -(n + α) 2 1 + β 2 (e iτ + e -iτ ) Y n e i(n+α)τ - β 2 (e iτ -e -iτ )(n + α)Y n e i(n+α)τ + N ∑ n s 1 + F 2 (e iτ + e -iτ ) Y n e i(n+α)τ = 0, (2.43) with s = Ω 2 0,B /ω 2 , which can be simplified into Y n (n + α) 2 + β 2 [Y n-1 (n -1 + α)(n + α) + Y n+1 (n + 1 + α)(n + α)] = s Y n + F 2 (Y n-1 + Y n+1 ) . (2.44)

Modelling the breakdown of Faraday waves

This constitutes a generalized eigenvalue problem of the form AX = sBX, where X is constructed from the real and imaginary parts of the vector Y. In order to restrict the problem to 0 ≤ n ≤ N, we follow [START_REF] Kumar | Parametric instability of the interface between two fluids[END_REF] who expressed the condition for c being real and truncated the solution. This leads to (2N+2)×(2N+2) matrix sizes for A and B. Lastly, we only focus on the subharmonic instability so we take α = 1/2 and we have the realisability condition Y -1 = Y * 0 , where * indicates the complex conjugate. The matrices A and B are thus:

A =                        d + 0 0 β 2 b 0 0 0 d - 0 0 β 2 b 0 . . . β 2 a 1 0 d 1 0 . . . . . . 0 0 β 2 a 1 0 d 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . d N-1 0 β 2 b N-1 0 . . . . . . 0 d N-1 0 β 2 b N-1 0 . . . β 2 a N 0 d N 0 0 β 2 a N 0 d N                        (2.45) and B =                        1 + F 2 0 F 2 0 0 1 -F 2 0 F 2 . . .
F 2 0 . . . . . . 0 1 0 F 2 0 . . . F 2 0 1 0 0 F 2 0 1                        (2.46) with a n = (n -1 + α)(n + α), b n = (n + 1 + α)(n + α).
For the diagonal terms in matrix A, we have d n = (n + α) 2 for all n > 0 and for n = 0, we have

d ± 0 = α 2 ± β 2 a 0 .
The results of this generalized eigenvalues problem are represented in the stability diagram of Fig. 2.13 exhibiting the instability tongues for β = 0 and β = 0.7. Note that the case β = 0 corresponds to the classical stability diagram of a Mathieu equation. In this representation, the right-hand sides of the neutral branches (solid black lines of Fig. 2.13) determine the critical amplitude of the primary wave and the beginning of the secondary instability as L 0 grows. This gives a critical threshold L crit which should be close to the one characterizing the wavebreaking L wb if the instability develops quickly (we still have L crit ≤ L wb ).

An analytical approximation for the critical threshold L crit can be derived by taking N = 0 and solving (for F and β ≪ 1):

(Ω 2 0B /ω 2 , F) with Ω 0B = (2AG 0 /L 0 ) 1/2 .
The two colored areas indicate the first subharmonic tongues with β = 0 and β = 0.7 respectively. The dashed blue lines correspond to the approximation (2.48) while the red dashed dotted ones correspond to (2.49).

1/4 -β/8 -s(1 + F/2) 0 0 1/4 + β/8 -s(1 -F/2) Y r 0 Y i 0 = 0 (2.47)
Indeed, one of the solutions of Eq. ( 2.47) is given by: First approximation of L crit from the global theory

L crit ∼ 2AG 0 (4 -2F) ω 2 (1 + β/2) , for F ≪ 1 and β ≪ 1. (2.48)
As shown in Fig. 2.13 in blue dashed lines, the criterion (2.48) slightly underestimates L crit at small F and large β while being reasonably correct for the parameters taken in the experiment. However, it becomes very bad at large F, even leading to negative values for F ≥ 2. Despite being more complicated, a better approximation can be derived, taking N = 1, corresponding to Eq. ( 2.49) below and shown by the red dashed dotted lines in Fig. 2.13.

Second approximation of L crit from the global theory

L crit = 2AG 0 ω 2 4F 2 + 8F -16 F(9 + 3β) -(β + 20) + M , (2.49) with M = 9(13 + 8β)F 2 -6(48 + β(17 + 4β))F + 256.
In equation ( 2.48), we see that the forcing parameter F together with the movement of the primary wave characterized by β contribute to the destabilization of the primary wave. In par-

Modelling the breakdown of Faraday waves

ticular, even without acceleration forcing (F = 0), the secondary instability can be triggered by the primary wave oscillations. This also leads to striking differences in terms of growth rates. For instance, for the subharmonic resonance at F = 0.7, we find a Floquet exponent µ = 0.09 for β = 0 (also corresponding to the growth of the primary wave), while µ = 0.2 for the case with β = 0.7. Therefore, the acceleration induced by the primary wave increases noticeably (but not drastically) the growth rate of the secondary modes. This effect explains why the breakdown occurs in a time scale shorter than the growth of the primary wave.

The local approach

The main advantage of the previous approach is of being relatively simple since it relies on the horizontal averaging process. However, the drawback of this theory is the loss of track of the physical mechanism responsible for the secondary instability. Moreover, it works under the assumption that the secondary instability mode only results from the interaction with the mean component of the primary wave, which can be excessive. This was our motivation to come up with a complementary method.

For the local approach, we perform a stability analysis of the flow generated at the node of the primary wave. To do so, we assume we have an inviscid interfacial wave of small amplitude. The equations describing the field generated by this type of wave have already been written in many classical textbooks (see for instance [START_REF] Sutherland | Internal Gravity Waves[END_REF]). The procedure leading to these equations is recalled below.

We first express a 2D incompressible velocity disturbance of the primary wave by its stream function as (u p , w p ) = (-∂ z ψ p , ∂ x ψ p ) and we seek modal solutions of the form ψ p (x, z, t) = A p (t) ψp (z)e ikx , assuming that the interface deformation is ζ p (x, t) = ξ p (t)e ikx . Neglecting the nonlinear terms, the momentum equations are written:

   ∂ t u p = -1 ρ ∂ x P ⇐⇒ ∂ t ∂ z ψ p = -ik ρ p, ∂ t w p = -1 ρ ∂ z P -G(t) ⇐⇒ ∂ t ikψ p = -1 ρ ∂ z p, (2.50)
with P = p + P Hydro = pρg(t)z with p of the form p = pe ikx . Deriving the first equation ( 2.50) along z, multiplying the second by ik and substracting the first equation to the second, we obtain the degenerate Rayleigh equation: (2.51) leading to ψp = e ±kz on either side of the interface. Then, discarding the second order terms, we write the kinematic condition expressing the interface dynamics: .52) giving ikA p = ξp and the continuity of ψ p across the interface. Moreover, the equation for the horizontal momentum, not taking into account the second order terms, writes: -

∂ 2 zz ψp -k 2 ψp = 0,
w p = ikψ p = Dζ p Dt = ζp , ( 2
∂ t u p = -∂ 2 tz ψ p = - 1 ρ i ikp, ( 2
ρ 1 G(t)ζ p + ρ 1 ∂ 2 tz ψ p ik = -ρ 2 G(t)ζ p + ρ 2 ∂ 2 tz ψ p ik . (2.54)
Leading to, in z = 0, ξp + AG(t)kξ p = 0.

(2.55)

Thus, for a standing wave, which is the sum of two traveling waves propagating in opposite directions, we have: kx), and u = ± ξp e ∓kz cos(kx).

ζ p = ξ p sin(kx), w = ξp e ∓kz sin(kx), ψ p = - ξp k e ∓kz cos(
(2.56)

The goal is to perform the stability analysis of this flow at the node of the primary wave where we have x = 0 and z = 0. Therefore, we rescale the dimensions with the typical wavenumber k wb of the secondary instability. Thus, the flow induced by the primary wave reduces to an horizontal interface subjected to a discontinuous tangential velocity, in the limit of small perturbation wavelength κ = k wb /k ≫ 1:

ζ p = 0, u p = ± ξp = ±U and w p = 0.
(2.57)

Now, we perform the linear stability analysis of the base flow coming from the primary wave, defined by Eq. ( 2.57). Using the same steps as previously, we consider a small velocity perturbation (u, w) = (-∂ z ψ, ∂ x ψ), we seek modal solutions of the form ψ ± (x, z, t) = A ± (t) ψ± e ik wb x with an interface deformation ζ(x, t) = ξ(t)e ik wb x and we obtain the degenerate Rayleigh equation:

∂ 2 zz ψ± -k 2 wb ψ± = 0, (2.58)
which gives, using boundary conditions, ψ+ = e -k wb z and ψ-= e k wb z .

At z = 0, we thus have, neglecting the nonlinear terms:

w = ik wb ψ = Dζ Dt = ζ ± Uik wb ζ, ( 2.59)
and so

A + = ξ ik wb + Uξ, A -= ξ ik wb -Uξ.
(2.60)

Moreover, still at the interface, the continuity of pressure gives:

p + -ρ 1 G(t)ξ = p -G(t)ξ, (2.61)
with p = p ± e ikx and the momentum equation for the horizontal velocity writes:

   ∂ t u + U∂ x u = -1 ρ 1 p, ∂ t u -U∂ x u = -1 ρ 2 p.
(2.62)

Modelling the breakdown of Faraday waves

Thus, we have:

-ik wb p + = ρ 2 k wb Ȧ+ + ik 2 wb U A + , -ik wb p -= ρ 1 -k wb Ȧ-+ ik 2 wb U A -. (2.63)
Combining all these conditions, we obtain:

ρ 1 ξ i + 2k wb U ξ + ξ[k wb U + ik 2 wb U 2 ] + ik wb ρ 1 G(t)ξ = ρ 2 - ξ i + 2k wb U ξ + ξ[k wb U -ik 2 wb U 2 ] + ik wb ρ 2 G(t)ξ. (2.64)
Thus, the equations driving the evolution of an interfacial perturbation amplitude ξ at the node of the primary wave are given in the inviscid limit by

ξ -2iAk wb U ξ + AG(t)k wb -k 2 wb U 2 -iAk wb U ξ = 0, (2.65)
already derived by Kelly [1965], which expresses the dynamics of an oscillating sheared interface, and within the Boussinesq limit here.

In this theory, supposing the instability amplitude and wavelength to be small, respectively kξ ≪ 1 and κ = k wb /k ≫ 1, allows us to neglect the interface tilting at the node and to assume the quasi-homogeneity of the perturbation. Note that often at least two vortices may appear in the experiments at the node of the primary wave during the breaking process suggesting the validity of the homogeneity assumption.

Many works have been dedicated to the stability analysis of (2.65) in the case of G constant and U oscillating at a single frequency ω. In addition to the parametric resonant modes Kelly [1965] identified, [START_REF] Lyubimov | Development of a steady relief at the interface of fluids in a vibrational field[END_REF] and [START_REF] Khenner | Stability of plane-parallel vibrational flow in two-layer system[END_REF] have shown and derived a criterion for the existence of Kelvin-Helmholtz type modes at the interface. These latter modes generate a longwave instability observed in most experiments and which may exhibit frozen wave patterns at high forcing frequency [START_REF] Wolf | Dynamic stabilization of the interchange instability of a liquid-gas interface[END_REF], [START_REF] Wunenburger | Frozen wave induced by high frequency horizontal vibrations on a CO 2 liquid-gas interface near the critical point[END_REF], [START_REF] Yoshikawa | Oscillatory Kelvin-Helmholtz instability. Part 2: an experiment in fluids with a large viscosity contrast[END_REF], [START_REF] Gaponenko | Dynamics of the interface between miscible liquids subjected to horizontal vibration[END_REF], [START_REF] Lyubimov | Viscosity effect on the longwave instability of a fluid interface subjected to horizontal vibrations[END_REF] and Gréa & Briard [2019]). Frozen waves are completely steady structures analogous to the inclined equilibrium positions of a strongly and horizontally oscillated pendulum.

As is emphasized from the local analysis, the shear is an important part of the breakdown process of the primary wave which was early recognized in the works of [START_REF] Thorpe | On standing internal gravity waves of finite amplitude[END_REF] and [START_REF] Kalinichenko | Development of a shear instability in nodal zones of a standing internal wave[END_REF]. However, the exact nature of the instability, being a Kelvin-Helmholtz (KH) or parametric resonance (PR) type, has yet to be identified. Indeed, the results of Kelly [1965] and [START_REF] Khenner | Stability of plane-parallel vibrational flow in two-layer system[END_REF] do not apply to our specific configuration where the acceleration oscillates at frequency ω while the shear velocity is subharmonic, it oscillates at frequency ω/2. Thus, we reconsider the problem of Eq. ( 2.65) by taking a primary wave of the form ξ p (t) = a cos(ωt/2) leading to: .66) where in the last expression, the subharmonic resonance condition ω 2 = 4AG 0 k/(1 + 4∆) for an inviscid interface is used. In order to study the stability of (2.65), it is further convenient to Chapter 2. The subcritical transition to turbulence of Faraday waves in miscible fluids

U(t) = - aω 2 sin(ωt/2) and U 2 = a 2 ω 2 8 (1 -cos(ωt)) = 1 2 AG 0 k 1 + 4∆ a 2 (1 -cos(ωt)), ( 2
introduce the new variable Y defined from ξ = Y exp[ t 0 iAk wb U(t ′ )dt ′ ] giving, at leading order in A: Ÿ + (AG(t)k wb -k 2 wb U 2 )Y = 0. (2.67)
Due to the change of variable expression and that U oscillates at ω/2, the response in ξ will also be subharmonic. However, for small Atwood number, ξ ∼ Y and the response can also be nearly synchronous. Indeed, by substituting the expression for U(t) into (2.67), we obtain a simple Mathieu equation of the form: In Fig. 2.14, we show the stability diagram of Eq. ( 2.67) in a κ -ka plane (the subharmonic resonance condition for the primary wave is again used). The instability tongue corresponding to the KH type modes appears for P ≤ 0 which stands as the classical criterion for the inviscid KH instability, k wb Ū2 ≥ AG 0 [START_REF] Chandrasekhar | Hydrodynamic and Hydromagnetic Stability[END_REF]). The parametric resonance zones start for P ≥ 0 but have also a continuation in the opposite half-plane. Remarkably, the instability zones exhibit a very weak dependence on κ for κ ≫ 1. Therefore, at given κ and as the primary wave amplitude grows, the perturbation passes through the successive instability zones, first the PR types then lately the KH one. The growth rates can be computed with the Floquet exponent µ and show a maximum approximately in the middle of each zone. The KH and PR1 growth rates are larger compared to the other instability zones. Therefore, the breakdown of the Faraday wave is expected to occur when the wave steepness ka lies in the instability KH or PR1 zones.

Ÿ + (P + Q cos(ωt))Y = 0, (2.68) with P = AG 0 k κ - 1 2(1+4∆) κ 2 (ka) 2 and Q = AG 0 k κF + 1 2(1+4∆) κ 2 (ka) 2 .

Data analysis of the experiments and simulations

This local theory is inviscid which explains that the growth rates are higher at large κ. Of course, the viscosity and the thickness δ of the interface should play a role and moderate this aspect.

Similarly to the global approach, we can propose an approximation for the critical wave steepness corresponding to the onset of the PR1 instability. Using the asymptotic expression for the neutral curves of the Mathieu equation in the limit P → -∞, Q → +∞ detailed in [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF], we obtain Eq. ( 2.69).

Critical steepness from the local theory

ka crit ∼ 1 3 (1 + 4∆)(1 + F). (2.69)
This simple expression (2.69) corresponds indeed to the plateau (it does not depend on κ) separating the PR1 and PR2 bands in the small perturbation wavelength limit κ ≫ 1 as shown on Fig. 2.14 (black dotted line).

Two theoretical approaches have been proposed to study the breaking process of Faraday waves. In order to assess their validity, it is necessary to detect the wavebreaking in both the experiments and simulations. In the next section, the detection method will be detailed.

Data analysis of the experiments and simulations

In this section, we detail the method used to detect the instant at which the wavebreaking occurs and how the primary wave amplitudes and inverse mean density gradients were measured. Then, we discuss the validity of the global and local approaches, developed in the previous section, against the data from the experiments and the simulations.

Primary wave amplitude and inverse mean density gradient

As was seen in Sec. 2.4, the primary wave is characterized, in the global theory, through the inverse mean concentration gradient L(t) at the interface rest position z = 0 whereas, in the case of the local theory, it is characterized through the wave amplitude ξ p (t) (see Fig. 2.12). Thus, to compare these two theories, we have to find the link between those two quantities. Assuming that the shape of the primary wave is frozen, a proportionality relation is expected between L(t) and ξ p whose proportionality coefficient is determined by the shape of the wave. Considering a sharp sinusoidal interface with wave amplitude larger than the interface thickness, we can write the interface position as ζ p (x, t) = ξ p (t) sin(kx) and thus show that the proportionality coefficient is π. Hence, we have

L(t) = π|ξ p |.
This choice was first motivated by the fact that a sinusoidal profile is a good fit for the interface position (see Fig. 2.15(a)). Though if the amplitude is too small, the fit is less satisfactory (probably because the interface thickness should be accounted for in the evaluation of the mean concentration profile). And secondly by the expression for the finite amplitude standing wave profile given by [START_REF] Thorpe | On standing internal gravity waves of finite amplitude[END_REF] in deep water approximation, knowing that the higher order corrections are negligible even at moderate wave steepness.

However, the inverse local gradient L is a difficult quantity to obtain from the mean concentration profile which is the only quantity to which we have access in the experiments. Indeed, these profiles can be noisy, especially when the secondary instabilities start, as they are obtained from the post-processing of the camera images as detailed in Sec. 2.1 and in [START_REF] Briard | The turbulent Faraday instability in miscible fluids[END_REF].

As for the amplitude of the wave |ξ p |, it can be obtained either from the crest-to-crest amplitude of the wave measured on the raw images or from taking the height between two arbitrary thresholds on the mean concentration (here we take the height where 0.1 ≤ C ≤ 0.9), which are two methods very sensitive to small variations. Therefore, we wanted to use a less sensitive method to measure indirectly both the inverse local gradient and the wave amplitude. Indeed, using the integral length L int introduced in Eq. ( 2.39) and, again, assuming the shape of the concentration profiles frozen, we can deduce all the characteristic lengths. For a sinusoidal interface, we have:

L int = 2.4|ξ p | = 0.76L.
(2.70)

Data analysis of the experiments and simulations

To ensure that these relations apply in our problem, we have measured the inverse mean concentration gradient L, the integral length L int and the amplitude |ξ p | at a local maximum just before the wavebreaking using the threshold method on C. In Fig. 2.15(b), we see that the sinusoidal profile provides a good evaluation of the mean concentration profile. The method shows that the inverse mean gradient is maximum at the center of the layer z = 0, in accordance with the fact that the wavebreaking starts at the nodes of the primary wave. Besides, we recover the expected relation between L and the wave amplitude |ξ p |. Therefore, the different lengths L and |ξ p | can be correctly evaluated from L int . In Fig. 2.15(c), we further check that the relation between L int and |ξ p | holds on multiple experimental data just before the wavebreaking. The correlation is thus satisfactory and gives us confidence in our method to extract the amplitudes and the inverse concentration gradients necessary to explore the wavebreaking phenomenon.

Moreover, the global theory assumes an inverse density gradient oscillating as

L(t) = L 0 (1 + β cos(ωt))
, while in the local theory the primary wave amplitude evolves as ξ p = a cos(ωt/2). By expanding |ξ p (t)| in Fourier series and truncating at leading order, we simply get L 0 = 2κa/π and β = 2/3. This allows to express the results from the global theory in terms of critical steepness as for the local theory, to better compare the two. In particular, the critical wave steepness for the global theory expressed by (2.48) in the limit of small F, using the resonance condition κ = π, gives:

Critical steepness from the global theory

ka crit = 3 8 (2 -F)(1 + 4∆).
(2.71)

Wavebreaking detection

Now that we are capable of measuring the amplitude of the primary wave, we have to find out at which instant the wavebreaking occurs. This part is used to detail the method for this wavebreaking detection. When the wavebreaking occurs, vortices appear at the node of the wave giving birth to a mushroom-like structure, as we can see in Fig. 2.16. Thus, one method of detection would be a simple visual criterion. However, this is somewhat subjective so we came up with an automated detection procedure using a simple algorithm based on the Thorpe displacement, first introduced by [START_REF] Thorpe | Turbulence and mixing in a Scottish Loch[END_REF]. Indeed, the wavebreaking phenomenon can be defined as a local overturning of the interface at the node and this displacement characterizes the distance a parcel of fluid has to move vertically in order to be in a stable equilibrium with the surrounding water. Thus, the existence of some local overturns would produce a non zero Thorpe displacement while a completely stable profile would have no displacements. Hence, for every image in an experiment, we compute the displacements δ T = z * -z for each 1D vertical transect, where z is the position of a fluid parcel on the instantaneous concentration field and z * is the position of the same parcel in the vertical sorted concentration field (see Fig. 2.16). Then, when this displacement exceeds a certain value (here |δ T | > 5 -12 px = 5.5 -13.2 mm), chosen a bit smaller than the initial interface widths, it gives us the image, and thus the instant, at which the wavebreaking occurs. For a given experiment, this algorithm provides almost the same image, and thus the same time and value for L, as the one chosen by eyes only (see Fig. it is not possible to detect from this method the beginning of the 'blurred' region as defined in [START_REF] Thorpe | On standing internal gravity waves of finite amplitude[END_REF] and indicating the secondary instability onset leading to the wavebreaking at the nodes. Now that we have the value of L for which there is wavebreaking, the last thing we need is to find the value of the ratio κ, and thus k wb , as it is necessary for the local theory. For this, we once again use the Thorpe displacement to measure the vortex size at wavebreaking. To determine k wb we take the maximum displacement δ T evaluated at the image given by the Thorpe displacement method. By construction, this measurement cannot be smaller than the arbitrary threshold chosen for the wavebreaking detection. In practice, δ T exceeds by more than three times this value.

As a summary, we illustrate in Fig. 2.17 the whole procedure allowing for the wavebreaking detection and the measurements of the wave amplitudes and inverse concentration gradients. The time evolution of the integral length L int in a simulation revealing the subharmonic oscillations (L int oscillates at frequency ω) and the growth of the primary wave are presented in Fig. 2.17(a). The instant corresponding to the wavebreaking is determined by evaluating the Thorpe displacement. This corresponds to the appearance of the vortices at the nodes. In addition, we show the evolution of the integral length L int,s computed from the sorted concentration profiles as in [START_REF] Briard | Harmonic to subharmonic transition of the Faraday instability in miscible fluids[END_REF]. This quantity expresses the evolution of irreversible mixing by distinguishing the available and the background potential energies (see also for instance [START_REF] Winters | Available potential energy and mixing in density-stratified fluids[END_REF], [START_REF] Peltier | Mixing efficiency in stratified shear flows[END_REF], Davies Wykes & Dalziel [2014]). We observe that the growth of the irreversible mixing starts just after the wavebreaking time. In Fig. 2.17, we verify further that this process is due to the wavebreaking as the curve of L int,s gets detached from the interface width L int,d expressing the thickening of the interface by pure diffusion. However, the slow evolution of L int,s would make it difficult to build a wavebreaking detection criterion from it. In Fig. 2.17(b), we also present the different mean concentration profiles before and after the wavebreaking time and renormalized by L int which confirms that the mean concentration profiles can be considered as frozen and well represented by a sinusoidal interface before the wavebreaking. At and after the wavebreaking, the mean concentration profiles are drastically distorted, and exhibit inversions of the stratification due to the roll-up of the interface. The predictions from the local and the global theories are also plotted in Fig. 2.17(a) showing a good agreement with the measured amplitude at the wavebreaking, and as expected we have L wb > L crit for both theories. We discuss more thoroughly this point in the next section.

Data analysis of the experiments and simulations

Chapter 2. The subcritical transition to turbulence of Faraday waves in miscible fluids

Results and discussion

Critical steepness values

The goal now is to compare the measurements of the wavebreaking in the experiments and in the simulations against the local and global approaches presented in Sec. 2.4.

Using the procedure described in the Sec. 2.5.2, we detect the moments corresponding to the wavebreaking. These moments are usually close to a maximum of the primary wave and give us the L int of the wavebreaking. However, it is the amplitude a that is needed to compare theories and experiments. In order to find this amplitude at wavebreaking, we perform a linear interpolation between two successive maxima. Thus, the wave steepnesses ka and wave ratios κ are reported and superimposed to the stability curves obtained from the local and the global approaches for both the experiments, in Fig. 2.18, and the DNS, in Fig. 2.19.

From these figures, we see that the wave steepness measured at wavebreaking is around ka ∼ 0.75, both for the experiments and the DNS. The critical steepness values predicted by the theories, for the moderate forcing parameters F investigated in this work, are both around ka = 0.5. Therefore, we notice that the experimental and numerical values are located above the thresholds provided by the theories and inside the parametric resonance instability band. Some points happen to be slightly below the threshold given by the global theory at ∆ = 0.04 in Fig. 2.18(d) and 2.19(d). These points correspond to negative detuning cases which have lower critical steepness value.

This result gives a strong credit to a wavebreaking process due to a secondary subharmonic instability appearing when the primary wave reaches a critical amplitude. Therefore, it is a subcritical instability. Moreover, we performed a numerical simulation, DNSe, whose parameters puts it close to the neutral curve (see Tab. B and Fig. 2.9). In this case, the primary wave reaches a saturation amplitude below the critical steepness and does not experience wavebreaking.

Another feature to be noticed from Fig. 2.18 is that the ka measured in the experiments tends to decrease when the wavenumber ratio κ increases. This stands despite the difficulty to evaluate κ but is less obvious in the DNS (see Fig. 2.19. Although the neutral curves depend very weakly on κ (they do not depend on it with the global theory), the instability growth rates evaluated from the Floquet theory increase when the small perturbation wavelength κ becomes large in the local theory. This is a characteristic result of classical inviscid theories such as for the Kelvin-Helmholtz instability. It is therefore not surprising to see that the wavebreaking occurs earlier, i.e. at smaller ka, when the secondary instability is initiated at smaller wavelengths. Moreover, the unstable modes from the higher parametric bands (PR2 for instance) may play a role in the destabilization process of the primary wave. More precisely, it may explain the blurring of the interface already observed by [START_REF] Kalinichenko | Development of a shear instability in nodal zones of a standing internal wave[END_REF], well before the wavebreaking, at ka ∼ 0.4. These PR bands are, however, expected to be seriously damped by the viscosity in addition of having a lower growth rate.

We also found that the growth rates of the PR1 modes, as computed in the local theory, are very large compared to the primary Faraday wave ones. Indeed, we find Floquet exponents µ varying between 2 and 8, corresponding to an amplification around 500 over one oscillation period. By contrast, the primary wave has Floquet exponents around 0.8. This is due to the additional forcing generated by the primary wave and explains why the breaking process develops over one or two oscillation periods. It seems also justified to neglect in the theory the growth of the primary wave amplitude for these parameters. By contrast, the growth rates evaluated by the global approach are much more modest (∼ 0.2). This point is probably due to the assumption linking the local mean concentration gradient to the amplitude of the primary wave detailed in Sec. 2.5.1. Indeed, as soon as the secondary instability develops, the shape of the mean concentration profile is distorted due to the perturbation feedback. This aspect is well demonstrated in the simulation result of Fig. 2.17. It eventually leads to a local change of sign in the stratification which would in return considerably increase the instability growth rate. However, this effect is not accounted for in the present global approach which still remains effective at predicting the secondary instability onset.

A large range of κ values has been found in the experiments and simulations, which comes equally from the primary and secondary wavelength measurements. The simulations provide larger κ (up to 110 in Fig. 2.19) partly because the primary wavelength is large. Therefore, an important question is: how is the wavenumber k wb selected? Neither the global theory, which does not depend on k wb , nor the local theory, which relies on the analysis of an infinitely thin interface, can explain this aspect. When the interface thickness reaches comparable size with the secondary instability wavelength, the natural frequency becomes bounded and the instability growth rate is expected to be limited. By analogy, the pure KH instability (represented by the Taylor-Goldstein equation) has a maximum growth rate around δk wb ∼ 1 at low Richardson number. This also depends on the shear and density profiles considered for the analysis (see for instance [START_REF] Taylor | Effect of variation in density on the stability of superposed streams of fluid[END_REF], [START_REF] Hazel | Numerical studies of the stability of inviscid stratified shear flows[END_REF], [START_REF] Caulfield | Multiple linear instability of layered stratified shear flow[END_REF]). The same process for the PR mode seems at work in our experiments as suggested by the measurements of k wb . This has been also evidenced theoretically by [START_REF] Poulin | Parametric instability in oscillatory shear flows[END_REF] on oscillatory uniform shear layer configurations. However, the local analysis accounting for the interface thickness is rendered complex by the coupling between the internal layer modes induced by the time varying shear and acceleration.

Data analysis of the experiments and simulations

Exploring further the forcing parameter effect

From Eqs. (2.69) and ( 2.71), we see that the local and global approaches differ in their dependence on the forcing parameter F. Indeed, the critical wave steepness grows linearly in F in the case of the local theory (see dashed line on Fig. 2.20) while, for the global theory, the dependence on F is more complex, it even exhibit a diminution of the wave steepness at small F (see dotted line on Fig. 2.20). Several DNS were performed, the DNSd cases, to investigate the dependence of the wave steepness on the forcing parameter and the results are shown on Fig. 2.20). These simulations are conducted using the same initial conditions and frequency ω, varying only the forcing parameter F and the results show a critical wave steepness globally evolving linearly with F. This seems reasonably well predicted by the local approach although it ends up overestimating the critical wave steepness at very large F (above F = 4). By contrast, the global approach clearly underestimates the critical wave steepness. However, when F become large, the growth of the primary wave amplitude cannot be neglected over one oscillation period, so the global and local approaches become limited. In the DNSd cases of Fig. 2.20, the wavebreaking occurs during the first oscillation for F ≥ 3, when the sign of gravity is inverted. Therefore, at large F, the secondary instability is expected to change nature as it becomes triggered more by the growth of the primary wave than its oscillations. In that respect, the wavebreaking process is no longer parametric but becomes of KH type similar to the one appearing in Rayleigh-Taylor instability [START_REF] Birkhoff | Helmholtz and Taylor instability[END_REF], [START_REF] Daly | Numerical study of two fluid Rayleigh-Taylor instability[END_REF], [START_REF] Baker | Singularity formation during Rayleigh-Taylor instability[END_REF]).

The final transition to turbulence

So far, we have not demonstrated that the secondary instability developing at the node of the primary wave triggers the full transition to a turbulent regime. Indeed, the simulations with a 2D interface initialisation are unable to develop the cascading process leading to turbulence. We thus provide as an example the simulation DNSf where the 2D initialisation is perturbed by a small random white noise in the spanwise y direction. In Fig. 2.21, we compare at different times DNSf with the simulation DNSd3 using the same parameters (see Tab. B) but having a simple 2D initial condition.

Both simulations evolve similarly until wavebreaking (here at ωt = 14.5). This means also that the critical amplitude is not modified by the small spanwise y perturbation, hardly visible since its amplitude is two orders of magnitude smaller than that of the x 2D perturbation. This feature has been reproduced on other simulations at different ϵ 1 parameters or in simulations where the spanwise invariance is broken by the presence of lateral walls. After wavebreaking, the simulations differ sensitively as DNSf (or simulations with the spanwise invariance broken) exhibits a rapid transition to turbulent mixing. The process is so violent that turbulence does not remain confined to the node of the Faraday wave and spreads around the whole layer. The images of Fig. 2.21 at ωt = 25.7 reveal that the final mechanism leading to turbulent mixing is related to the merging and the stretching of the secondary vortices at the node forced by the oscillations of the Faraday wave.

In order to identify the specificity of this transition, we evaluate the Reynolds and bulk Richardson numbers as classically introduced for the Kelvin-Helmholtz instability [START_REF] Caulfield | Layering, instabilities, and mixing in turbulent stratified flows[END_REF]). Taking for the interface half width σ and for the half shear velocity U = aω/2, we find Re = Uσ/ν = 313 at the wavebreaking (ωt = 14.5). For the bulk Richardson number, we obtain

Conclusion

Ri = AG 0 σ/U 2 = 0.023. While the low value of the Richardson number is consistent with the existence of a strong shear instability, the Reynolds number is unexpectedly small for a mixing transition. This is despite being probably underestimated due to the blurring of the interface before the wavebreaking. Perhaps the fact that the transition does not result from a convective secondary instability inside the vortices as often observed for the Kelvin-Helmholtz instability [START_REF] Salehipour | Turbulent diapycnal mixing in stratified shear flows: the influence of Prandtl number on mixing efficiency and transition at high Reynolds number[END_REF]) may explain this aspect. By contrast, we can consider a global Reynolds number based on the width and velocity of the layer. This gives Re = a 2 ω/ν ∼ 10 4 and agrees with a mixing transition criteria of Re > 10 4 proposed by [START_REF] Dimotakis | Turbulent mixing[END_REF].

Conclusion

In this chapter, the wavebreaking mechanism leading to turbulence of growing Faraday waves at the interface between miscible fluids of small density contrast was investigated experimentally, numerically and theoretically along with the mode selection mechanism of the primary wave.

We have evidenced that the mode selection of the Faraday wave results not only from a linear process but also a nonlinear competition favoring the modes with smaller wavelength. Indeed, it was shown that, considering the large forcing accelerations and quantification induced by the geometry of the tank used in our experiments, several subharmonic modes can be simultaneously unstable by the Floquet linear theory but only one appears. The parameters taken in the experiments are located well above the instability threshold even if we take into account the viscous damping at the walls and from the bulk flow. Hence, the inviscid theory provides a good first approximation for this problem. Moreover, as the mode amplitude grows, the natural frequencies of the system decrease. Therefore, when the primary amplitude reaches critical value, the modes with positive detuning, called subcritical, are favored whereas the modes with negative detuning, called supercritical, become damped as they are less in parametric resonance with the forcing frequency. This mode competition phenomenology also explains the sensitivity to initial conditions and the symmetry breaking of the dominant mode which has been observed in our experiments.

Two theories were constructed to explain the breakdown of the primary wave as the appearance of a secondary subharmonic instability at small scales. The first theory, referred to as 'global', is based on a horizontal averaging process and reveals that the secondary instability is mostly due to the oscillations of the primary wave. This explains why this instability develops very rapidly compared to the growth rate of the primary wave. From this theory, we propose a criterion giving the critical steepness of the primary wave at which the wavebreaking is expected to occur. However, this theory is too simple and it cannot explain why the secondary instability appears at the node of the wave and how it depends on its wavelength. This is why we thought of a second theory, referred to as 'local'. In this theory, we study the flow in a local frame attached to the node of the wave. It reveals the importance of the shear in the development of the secondary instability. From a stability analysis, it was shown that the unstable modes can be either Kelvin-Helmholtz or of parametric resonance type, with the latter developing earlier during the growth of the primary wave.

Then, we compared these theories to the experimental and numerical data. Thus, we evaluated the primary wave amplitudes and the wavenumber associated with the wavebreaking Chapter 2. The subcritical transition to turbulence of Faraday waves in miscible fluids thanks to the Thorpe displacement method which indicates a local overturning. The results show that the wavebreaking is detected for wave steepnesses around 0.75 inside the parametric resonance band. Thus, it was deduced that the wavebreaking came from a parametric resonant subharmonic instability developing when the primary Faraday wave amplitude reaches a critical value: it is of subcritical nature. Therefore, the present mechanism for Faraday waves shows similarities with the breaking process of internal gravity waves. However, when the forcing parameter is increased, our approaches come to a limitation as the wavebreaking process changes in nature resulting more from the amplitude growth of the primary wave than its oscillations. In this case, the phenomenon becomes similar to the secondary vortices appearing in the classical Rayleigh-Taylor instability.

In addition, to show that this wavebreaking does lead to turbulence, numerical simulations with 2D initial conditions perturbed along the spanwise direction were performed. It was indeed found from these simulations that the final transition to turbulence originates from the secondary instability at the node of the wave. The stretching and merging of the secondary vortices driven by the oscillations of the primary wave is an efficient mechanism to produce mixing, differing substantially with the transition scenarios observed in the context of the Kelvin-Helmholtz instability.

Description of the instability

Introduction

The second part of this thesis is dedicated to the Rayleigh-Taylor instability through a grid. Thus, this chapter is aimed at describing the classical Rayleigh-Taylor instability and its main features. Moreover, a small review of the experimental and numerical studies and some recent findings are given.

Description of the instability

Unlike the Faraday instability, the Rayleigh-Taylor instability, also involving the superposition of two fluids of different densities, is not a parametric instability. Indeed, in this case the heavier fluid is on top of the lighter one in a vertical downwards gravitational field. Thus, it is naturally unstable to any sort of perturbations at its interface (in a case without viscosity or surface tension) and is solely driven by the gap in potential energy coming from the difference in density of the two fluids. The bigger this difference, the more rapidly the instability develops. However, it bodes well to remember that if there are no perturbations at the interface between the two fluids, the hydrostatic pressure will maintain the system.

The Rayleigh-Taylor instability (see Fig. 3.1b) has been extensively studied throughout history starting with the observations of [START_REF] Rayleigh | Investigation of the character of the equilibrium of an incompressible heavy fluid of variable density[END_REF] and [START_REF] Taylor | The instability of liquid surfaces when accelerated in adirection perpendicular to their planes. I[END_REF] describing the very characteristic development of this instability: we have 'bubbles', or mushroom-like shapes, of lighter fluid ascending inside the heavier fluid while 'spikes' of heavier fluid fall inside the lighter fluid, thus creating a mixing layer expanding away from the initial position of the interface between the two fluids (see Fig. 3.1a).

At small Atwood number (A < 0.1, arbitrary threshold), defined as A = (ρ hρ l )/(ρ h + ρ l ), the up-welling (bubble) and down-welling (spike) penetrations are almost the same, thus the disturbances form a symmetrical mixing region. Therefore in small Atwood incompressible cases, the Boussinesq assumption is used. This assumption states that the density can be taken as constant everywhere in the equations except in the buoyancy terms where the density difference effects need to be accounted for. However, at large Atwood number, the spike narrows, the drag on its head decreases and thus, its penetration is greater than that of the bubble. The tip of the structures are destabilised by the Kelvin-Helmholtz instability, involving the local shear, which leads to mixing and turbulence.

The reason behind the extensive study of this instability is that it can be found in several natural instances like the cirrus clouds where characteristic mushroom-like shapes can be seen. It also plays an important role in many scientific domains such as oceanography, climate dynamics, astrophysics (for example in the filamentary structure on the sun) and in many industrial applications such as the Inertial Confinement Fusion or painting (a more detailed list can be found in Zhou [2017a] and Zhou [2017b] where the instability is reviewed).

Self-similar regime and turbulence

The evolution of a Rayleigh-Taylor mixing layer can be divided into three stages: first the linear phase, then the nonlinear growth and finally, self-similar growth. The linear phase, in which the initial perturbations grow exponentially and independently until they reach a critical amplitude called the saturation amplitude, can be predicted by linear stability theory. The nonlinear growth results from the mode coupling which leads to the development of an inertial range of scales in the kinetic energy spectrum and ultimately the mixing transition. This mixing transition occurs for large Reynolds numbers, Re = Ul/ν > 1 -2 × 10 4 , where U is the velocity, l is the tranverse extent of the flow and ν is the kinematic viscosity [START_REF] Dimotakis | The mixing transition in turbulent flows[END_REF]). The self-similar growth occurs after the Rayleigh-Taylor mixing layer transitions to fully developed turbulence [START_REF] Cook | Transition stages of Rayleigh-Taylor instability between miscible fluids[END_REF] and [START_REF] Cook | The mixing transition in Rayleigh-Taylor instability[END_REF]).

In mathematics, a self-similar object is similar to itself at a different scale. This property can be found in everyday objects such as the Romanesco broccoli (Fig. 3.2). Indeed, the same pattern is found in the Romanesco broccoli as one zooms in on the picture. For the Rayleigh-Taylor instability, the temporal growth of the mixing layer is self-similar with a well known expression in the Boussinesq approximation given by [START_REF] Youngs | Numerical simulation of turbulent mixing by Rayleigh-Taylor instability[END_REF]:

Self-similar regime and turbulence

Mixing layer width

L(t) = 2αAG 0 t 2 , (3.1)
where α is the bubble growth rate. This expression was popularised by [START_REF] Youngs | Numerical simulation of turbulent mixing by Rayleigh-Taylor instability[END_REF]. However, [START_REF] Birkhoff | Taylor instability and laminar mixing[END_REF] first derived a self-similar growth for mixing layer that is consistent with Eq. (3.1). Indeed, he found that, in the long term, the penetration distance h b of the bubbles in the heavy fluid can be written h b = α b AG 0 t 2 , α b being the growth rate of the bubbles, while the penetration distance of the spikes is h s = α s AG 0 t 2 , where α s is the growth rate of the spikes. Thus, using the notation of Eq. ( 3.1), we have

L(t) = (α b + α s )AG 0 t 2 .
The height h b is given by the saturated linear mode with the biggest amplitude. Moreover, he supposed that the spikes fall with a free fall velocity such that their height is proportional to G 0 t 2 . Self-similarity in the Rayleigh-Taylor instability is often associated with the large-scale properties of the flow, dominated by production mechanisms, which were largely studied in CEA [START_REF] Gréa | The rapid acceleration model and the growth rate of a turbulent mixing zone induced by Rayleigh-Taylor instability[END_REF], Soulard et al.

[2015] and [START_REF] Soulard | Influence of the mixing parameter on the second order moments of velocity and concentration in Rayleigh-Taylor turbulence[END_REF]).

The growth rate of the mixing layer is a critical parameter, not only for engineering applications, but also because it characterises the self-similar dynamics of the Rayleigh-Taylor instability. It will be evaluated, using the Boussinesq assumption [START_REF] Glimm | A critical analysis of Rayleigh-Taylor growth rates[END_REF] and [START_REF] Jacobs | Rayleigh-Taylor instability in complex stratifications[END_REF]), as:

Growth rate α = L2 8AG 0 L . ( 3.2)
Other techniques exist to compute this growth rate and were compared in [START_REF] Cabot | Reynolds number effects on Rayleigh-Taylor instability with possible implications for type-Ia supernovae[END_REF] as shown in Fig. 3.3a: the similarity method (Eq. (3.2)) is more robust than other evaluations such as α = L/AG 0 t 2 , or fitting with polynomials, or using a criterion based on a concentration threshold.

The value of the growth rate is still a much discussed topic due to the large discrepancies found in the literature. However, the work of the α-group [START_REF] Dimonte | A comparative study of the turbulent Rayleigh-Taylor instability using high-resolution three-dimensional numerical simulations: The Alpha-Group collaboration[END_REF]) should be mentioned here. It was a collaboration between several teams working to find a value for the growth rate of the mixing layer when the self-similar regime is reached. It allowed the comparison between experimental results and identical simulations performed with different numerical codes. They observed a clear difference between the experimental and numerical values for the growth rate, namely around α ∼ 0.025 while in the experiments α ∼ 0.055. Moreover, the miscibility seems to strongly affect the value of the growth rate α [START_REF] Roberts | The effects of forced small-wavelength, finite-bandwidth initial perturbations and miscibility on the turbulent Rayleigh-Taylor instability[END_REF]) indicating that it plays a role in the late-time turbulent instability growth. Indeed, the growth rate for the miscible experiments of [START_REF] Roberts | The effects of forced small-wavelength, finite-bandwidth initial perturbations and miscibility on the turbulent Rayleigh-Taylor instability[END_REF] was found to always be significantly smaller than the growth rate of the immiscible experiments. LEM is from Linear Electric Motor [START_REF] Dimonte | Density ratio dependence of Rayleigh-Taylor mixing for sustained and impulsive acceleration histories[END_REF]), RR from 'rocket rig' [START_REF] Read | Experimental investigation of turbulent mixing by Rayleigh-Taylor instability[END_REF]), K from Kucherenko et al. [START_REF] Kucherenko | Experimental study of the gravitational turbulent mixing self-similar mode[END_REF]), AS is from Andrews and Spalding [START_REF] Andrews | A simple experiment to investigate two-dimensional mixing by Rayleigh-Taylor instability[END_REF]).

Experimental works and Numerical simulations

The Rayleigh-Taylor instability has been extensively studied through experimental works (for example [START_REF] Schneider | Large and small scale structure in Rayleigh-Taylor mixing[END_REF], Ramaprabhu & Andrews [2004] and [START_REF] Akula | Effect of shear on Rayleigh-Taylor mixing at small Atwood number[END_REF]) and numerical studies (for example [START_REF] Cook | Transition stages of Rayleigh-Taylor instability between miscible fluids[END_REF], [START_REF] Cook | The mixing transition in Rayleigh-Taylor instability[END_REF] and [START_REF] Glimm | A numerical study of bubble interactions in Rayleigh-Taylor instability for compressible fluids[END_REF]). However, it can be difficult to control and characterize the initial perturbation in these experiments which is the principal drawback of experimental investigations. This issue can be resolved using computational studies as the initial conditions are more easily controlled and characterized. Both methods have their advantages and drawbacks and one must complement the other. The second part of this thesis being dedicated to the study of the turbulent mixing, a sample of large-Reynolds-number Rayleigh-Taylor experimental and numerical studies are reviewed hereafter. For a more comprehensive list, see the reviews Zhou [2017a] and Zhou [2017b].

Experimental works

A classic experiment in the Rayleigh-Taylor research field is the 'Rocket Rig' (RR) of [START_REF] Read | Experimental investigation of turbulent mixing by Rayleigh-Taylor instability[END_REF], which consists of two-dimensional or three-dimensional boxes filled with the lighter fluid above the heavier one. These are then accelerated downwards with an acceleration up to 50 times that of the gravity which sets up an unstable Rayleigh-Taylor interface (see Fig. 3.4a for the experimental setup). However, only high Atwood numbers were considered and most experiments were made with immiscible fluids. Indeed, the interface in RR experiments with miscible fluids would smear due to diffusion leading to a delay in the development of the mixing layer. These Here, the facility allows for an acceleration up to a thousand times the earth's gravity and variable acceleration profiles unlike the Rocket Rig. Both of these experiments have the particularity to possess a non constant gravitational field. Small Atwood number experiments were not carried out with these setups as it would require for the box to move a much higher distance S to reach the late-time mixed state than with high Atwood number. This accelerating distance could easily become a difficulty as we have S > 300H, with H the box height, for Atwood numbers smaller than 0.01.

Hence, stationary Rayleigh-Taylor experiments were created to explore small Atwood numbers. These use miscible fluids and innovative methods to place the heavy fluid over the light one and three of them are listed hereafter.

The overturning tank

The first stationary experiment is the overturning tank engineered by [START_REF] Andrews | A simple experiment to investigate two-dimensional mixing by Rayleigh-Taylor instability[END_REF] (see Fig. 3.5). The experimental setup is shown in Fig. 3.5a and consists of a narrow tank of thickness 0.5 cm, height 36 cm and breadth 25 cm filled with fresh water on top of a heavy brine solution. It pivots around its center and when it is released, it spins 180 • thanks to springs and is caught at the top by a door catch thus placing the heavy fluid on top of the lighter one. The experiments could be performed either with a horizontal interface or with an initially tilted rig with a wedge placed on one side. The tilted rig generates an initial long wavelength of half a sawtooth. The result of this tilted-rig experiment for an initial angle of 0.9 • is shown in Fig. 3.5b at four different times with an Atwood number of A = 0.048. In the beginning, the center of the mixing region develops as if it were a flat interface with visible bubbles and spikes.

Then, an overturning motion appears with the creation of jets at the sides of the box. As this overturning motion develops, the central mixing zone is rotated and stretched out due to the longer wavelength initial perturbations. This leads to a thinning of the central mixing region.

Thus, this new type of experiment found a way to overcome the issue of mixing development time at small Atwood number by inverting the position of the two fluids and holding the tank steady until the mixing was complete. This was used by [START_REF] Voropayev | Experiments on the evolution of gravitational instability of an overturned, initially stably stratified fluid[END_REF] to explore the Rayleigh-Taylor instability in a linearly stratified fluid and by Dalziel et al. [2008] to study the Rayleigh-Taylor mixing in a high aspect geometry. However, the narrow depth of the tank imposed a two dimensional flow although it is necessary to prevent sloping during the overturn. Moreover, the initial condition is not well controlled as catching the tank at the top causes it to shake, which introduces small perturbations. These perturbations were found to have higher wavenumbers than the most unstable one and thus they were not significant. Other issues were found with these experiments and are listed in [START_REF] Andrews | Small Atwood number Rayleigh-Taylor experiments[END_REF].

Experimental works and Numerical simulations

The sliding barrier

The second experiment is the sliding barrier experiment (see Fig. 3.6a) which was developed by [START_REF] Linden | Molecular mixing in Rayleigh-Taylor instability[END_REF] in Cambridge University in order to provide a low Atwood number three dimensional setup. This experiment consists in separating initially the two fluids with a barrier that is later on removed. Thus, the lighter fluid (fresh water) is placed beneath the barrier while the heavier fluid (salt water) is put above it inside a cuboidal tank. Equal volumes are occupied by the fluids. Then, the barrier is removed by sliding it through a slot at one end of the tank (see Fig. 3.6b). The initial experimental setup introduced by Linden & Redondo [1991] had a simple metal barrier which was then replaced by a pair of horizontal plates each wrapped in nylon fabric (see Fig. 3.6c). The goal of this new type of barrier was to eliminate the viscous boundary layers that forms on the upper and lower surfaces of the barrier as it is withdrawn. Indeed, the plates are removed at speed U by pulling on the nylon fabric from the tank at speed 2U, confining the shear between the plates and the nylon.

The two main issues with this experiment are: the finite time of the barrier withdrawal and the asymmetrical filling of the void left after the removal of the plates by the upper fluid. The former introduces a delay in the contact between the two fluids. The latter initiates a twodimensional vortex sheet seeding a large-scale two-dimensional component persisting throughout the experiments. This can result in a jet that is accelerated as can be seen in Fig. 3.7 down the right-hand wall in an experiment made by [START_REF] Dalziel | Self-similarity and internal structure of turbulence induced by Rayleigh-Taylor instability[END_REF].

Despite these drawbacks, this experimental setup was quite recently used by Holford et al.

The water tunnel

The third experiment is that of the water tunnel (see Fig. 3.8) developed by Andrews and collaborators [START_REF] Snider | Rayleigh-Taylor and shear driven mixing with an unstable thermal stratification[END_REF], [START_REF] Wilson | Spectral nonequilibrium in a turbulent mixing layer[END_REF] and [START_REF] Wilson | Spectral measurements of Rayleigh-Taylor mixing at small Atwood number[END_REF]) at Texas A&M University.

In this experiment, the heavy fluid (cold water with typical temperature 20 • C) is put over a splitter plate and the lighter fluid (warm water with typical temperature 25 • C) is below as shown in Fig. 3.8. Hence, at the end of this plate, the lighter fluid is below the heavier one and the mixing layer brought by the Rayleigh-Taylor instability develops downstream. The densities of both fluids can be computed through an equation of state for water. Thus, in this case the Atwood number is 7.4×10 -4 . It results in a spatially evolving mixing layer as the fluids flow downstream as shown in Fig. 3.9.

Experimental works and Numerical simulations

Numerical simulations

As said previously, the numerical simulations' main advantage is the control over the initial conditions that is not easily achieved in experiments.

The first direct three-dimensional numerical simulations of the Rayleigh-Taylor instability have been performed by [START_REF] Youngs | Numerical simulation of turbulent mixing by Rayleigh-Taylor instability[END_REF], [START_REF] Youngs | Three-dimensional numerical simulation of turbulent mixing by Rayleigh-Taylor instability[END_REF] and [START_REF] Youngs | Numerical simulation of mixing by Rayleigh-Taylor and Richtmyer-Meshkov instabilities[END_REF] to investigate turbulent mixing due to this instability. It was notably found in [START_REF] Youngs | Three-dimensional numerical simulation of turbulent mixing by Rayleigh-Taylor instability[END_REF] that the value of α for miscible fluids may have been overestimated when evaluated from the immiscible fluids experiments of [START_REF] Read | Experimental investigation of turbulent mixing by Rayleigh-Taylor instability[END_REF] or from the miscible fluids experiments of [START_REF] Linden | Molecular mixing in Rayleigh-Taylor instability[END_REF]. Indeed, in those experiments they found α = 0.07 whereas [START_REF] Youngs | Three-dimensional numerical simulation of turbulent mixing by Rayleigh-Taylor instability[END_REF] found α = 0.04 (see Fig. 3.10).

As mentioned before, the value of the growth rate was studied by the α-group [START_REF] Dimonte | A comparative study of the turbulent Rayleigh-Taylor instability using high-resolution three-dimensional numerical simulations: The Alpha-Group collaboration[END_REF]) through an extensive numerical work using a variety of high-resolution, multimode, three dimensional numerical simulations. Indeed, seven different codes named TUR-MOIL3D, FLASH, WP/PPM, NAV/STK, RTI-3D, HYDRA and ALEGRA from five different institutes (AWE, U. Chicago, LLNL, Texas A&M and Sandia NL) were compared, each with their own methods (Eulerian, Piecewise Parabolic Method, NS and Arbitrary Lagrangian Eulerian approach) to study the turbulent Rayleigh-Taylor instability. They also used this important numerical alternating phases of unstable and stable stratification, initially produces a growing turbulent mixing layer. However, after a few periods of the acceleration, an asymptotic state is reached (see Fig. 3.12) in which the turbulence decays in time and becomes unable to develop further instabilities. The finite extension of the mixing layer in this asymptotic state was found to have an amplitude proportional to T 2 where T is the period of the acceleration.

This thesis focuses on the dynamics of the Rayleigh-Taylor instability through a grid. As such, the study of the Rayleigh-Taylor instability through a porous media was of particular interest and is described hereafter.

Rayleigh-Taylor instability in porous media

A recent study by [START_REF] Boffetta | Dimensional effects in Rayleigh-Taylor mixing[END_REF] addressed the effects of the dimensional confinement on the evolution of the incompressible Rayleigh-Taylor mixing. This study used numerical simulations both in a bulk flow and in a porous medium, either in a three dimensional case or in an ideal two dimensional case where the confinement is not caused by the presence of side walls. In the bulk flow, the dynamics is ruled by the Navier-Stokes equation whereas, in the porous medium, it follows the Darcy equation. The vertical sections of the density field are shown in Fig. 3.13 for these simulations.

It is evident that the structures are very different between a classical Rayleigh-Taylor instability (NS system) and a porous one (Darcy system). Indeed, in the NS configuration ((a) and (b) of Fig. 3.13), we have the expected bubbles and mushrooms, while in the Darcy configuration ((c) and (d) of Fig. 3.13), there are only thin elongated 'fingers'. Thus, the turbulent properties and the mixing parameter are very likely to be very different as well. However, it was found that both configurations show similar dependence on the dimensionality. Indeed, at fixed time, for both configurations, the structures are more elongated in two dimensions and less blurred than in three dimensions, meaning there is less mixing due to the faster growth of the large structures. It was found that, in two-dimensional geometry, the mass flux is more intense and the correlation between the vertical velocity field and the density field is stronger, giving a bigger Nusselt number Nu, representing the dimensionless mass flux, for the same Rayleigh number Ra. This increase correlation makes for a faster growth rate of the mixing layer in two dimensions. Indeed, it is demonstrated by the authors that the growth rate of the mixing zone width is larger by about 57% in the 2D simulations than in the 3D simulations in both the Navier-Stokes and Darcy simulations. The α-group also reported a higher growth rate in 2D than in 3D [START_REF] Dimonte | A comparative study of the turbulent Rayleigh-Taylor instability using high-resolution three-dimensional numerical simulations: The Alpha-Group collaboration[END_REF]). The downside of this larger correlation and faster growth rate is the loss of homogeneisation of the density fluctuations within the mixing layer in the two-dimensional geometry. This was shown by the much larger variance of the density fluctuations in 2D compared to 3D (around 100% larger).

Introduction

The Rayleigh-Taylor instability, described in Chapter 3, is known to be extremely hard to study experimentally. Indeed, briefly, it consists of the superposition of a heavier fluid on top of a lighter one in a gravitational field [START_REF] Taylor | The instability of liquid surfaces when accelerated in adirection perpendicular to their planes. I[END_REF]). For ideal fluids (no viscosity nor surface tension), this configuration is unstable to any perturbations. Thus getting the heavier fluid on top of the lighter one, without perturbing the interface between the two fluids, is not an easy task. Moreover, this instability dynamics is very sensitive to the initial conditions even in the late time dynamics (Zhou [2017b]). This means that it is quite crucial to control perfectly these initial conditions. However, most Rayleigh-Taylor experiments [START_REF] Linden | Molecular mixing in Rayleigh-Taylor instability[END_REF], [START_REF] Lawrie | Rayleigh-Taylor mixing: confinement by stratification and geometry[END_REF]) introduce a shear at the very beginning of the instability, thus the initial conditions are not perfectly controlled. A good control of the initial conditions is important as it impacts the growth rate in the long term.

In an attempt to better control the initial conditions, a new type of Rayleigh-Taylor experiment is presented in this chapter (see Fig 4.1). The idea is to use a grid to maintain the heavier fluid on top of air in a reversed container. Then the heavier and lighter fluids are put into contact allowing for the instability to develop. The grid remaining at the interface between the two fluids throughout the duration of the experiment, it is important to assess its impact on the overall dynamics. This chapter is organised as follows. First, the experimental setup, the linear stability analysis for the two interfaces of the system (water-air and brine-fresh water) and the experimental data are presented in Section 4.1. The code and different numerical parameters used are afterwards given in Section 4.2 and lastly, a numerical study to investigate the effect of the grid is detailed in Section 4.3.

Chapter 4. The Rayleigh-Taylor instability through a grid: Influence of the grid

Experimental work

The experiments were made at the Ecole Centrale de Lyon with the collaboration of Louis Gostiaux and Marc Michard of the LMFA (CNRS, UMR) and the idea of using a grid first came from Olivier Soulard (CEA, DAM, DIF).

Principle and Linear Stability Analysis

The entire experiment is based on the use of a grid and surface tension to maintain a body of brine inside a reversed container. This means that a stable configuration, involving the superposition of salt water on air, must be obtained. This can only occur if the holes of the grid are small enough that the surface tension becomes strong enough to stop the brine from flowing through them. Performing the linear stability analysis of the salt water-air interface, we can find the critical wavelength above which this configuration becomes stable. Considering a 2D incompressible and inviscid system containing two uniform inviscid and immiscible fluids of constant density separated by a horizontal boundary at z ∼ 0 (see Fig. 4.2), the linearised equations of motion for the perturbation fields of small amplitude in each fluid are [START_REF] Chandrasekhar | Hydrodynamic and Hydromagnetic Stability[END_REF]):

   ρ∂ t u = -∂ x p, ρ∂ t w = -∂ z p -G 0 ρ, (4.1)
where p is the perturbation around the hydrostatic state of the pressure field G 0 is the gravitational acceleration, u is the x-velocity component and w is the z-velocity component. The density ρ is considered constant in each fluid and different from the density of the other fluid.

The incompressibility condition is given by:

∂u ∂x + ∂w ∂z = 0. (4.2)
We seek modal solutions of the stream function ψ, defined as (u, w) = (-∂ z ψ, ∂ x ψ), of the form:

ψ = ψ(z) exp(ikx + Ωt), (4.3)

Experimental work

where Ω is the growth rate of the perturbation and k is the horizontal wavenumber. We also write p = p(z)e ikx . Rewriting equations (4.1) for the stream function, we obtain the degenerate Rayleigh equation:

∂ 2 zz -k 2 ψ = 0. (4.4)
The solution of this equation is of the form ψ+ = Ae -kz for z > 0 and ψ-= Ae +kz for z < 0. The constant A is chosen to ensure the continuity of ψ across the interface (kinematic condition written in Eq. ( 4.5)). Moreover, the normal velocity w = ∂ x ψ must also be continuous at the interface.

The interface deformation is written ξ = ξ(z) exp(ikx + Ωt) and its amplitude is assumed to be very small so we have ξ ≪ 1. Thus, neglecting the second order terms, the kinematic condition applied in z = 0 becomes:

D ξ Dt = ∂ ξ ∂t = w = ik ψ ⇐⇒ ξ = ik ψ Ω . (4.5)
We know that the surface tension must play an important role for the stability of a water-air interface, thus we need to take it into account in our analysis. For this, we write the Young-Laplace equation at z = 0 giving:

∆P = γ T (∂ 2 xx ξ) (4.6)
where ∆P = P + -P -is the pressure jump and γ T is the surface tension. The solution of Eq. ( 4.4) remains unchanged by the addition of the surface tension as it will only impact the boundary condition at the interface which is derived form the Young-Laplace equation. Moreover, from Eq. ( 4.1), we have ∂ 2 zz ψ = ikp/ρ. Combining these with the Young-Laplace equation, we find that the boundary condition that must be satisfied at the interface writes:

ρ h ∂ z ψ + 0 -ρ l ∂ z ψ - 0 = - k 2 Ω 2 [G 0 (ρ h -ρ l ) -k 2 γ T ]ψ 0 , (4.7)
with ρ h the density of the heavier fluid, ρ l the density of the lighter fluid and ψ 0 is the value of

ψ at z = 0.
Combining the solutions of Eq. ( 4.4) and Eq. ( 4.7), we find:

Dispersion relation with surface tension (see [START_REF] Chandrasekhar | Hydrodynamic and Hydromagnetic Stability[END_REF]) system. Moreover, taking into account the surface tension gives a maximum growth rate Ω 0 at a certain wavenumber k 0 . For this most unstable mode the amplitude of the disturbance grows the most rapidly. The wavenumber k 0 can be evaluated from the critical wavenumber as it is given by k 0 = k c / √ 3 and thus the maximum growth rate is given by:

Ω = G 0 k A - k 2 γ T G 0 (ρ h + ρ l ) 1/2 . ( 4
Ω 0 = 2 3 3/2 (ρ h -ρ l ) 3/2 G 3/2 0 (ρ h + ρ l )γ 1/2 T . ( 4
.9)

In our case, we have k 0 = 211.98 m -1 and Ω 0 = 37.20 s -1 . It gives a most unstable wavelength of

λ 0 ∼ 3 cm.
The critical wavenumber given in Fig. 4.3a is k c = 367.34 m -1 meaning that the smallest wavelength that can destabilize the system is λ c = 1.7 cm. Thus, by choosing a grid whose holes are smaller than this value, the system brine-air should be stabilised by the surface tension, as the wavelengths larger than the grid mesh size cannot develop.

As a conclusion, we have seen that salt water can indeed be maintained on top of air using a grid and the surface tension. However, the effect of the surface tension will no longer retain the brine once it is put into contact with the fresh water, thus letting a Rayleigh-Taylor like instability develop. This principle allows for a new type of experimental setup using the grid for the filling and maintaining of the heavier fluid. The experimental setup is explained hereafter.

Experimental setup

The experimental setup consists of two cylindrical tanks one on top of the other of height H cyl = 30 cm and diameter 2R cyl = 19.4 cm (see Fig. 4.4a). The one on the bottom is fixed and filled from below, using a pump, with a mix of fresh water, ethanol and sulforhodamine. The one on the top is movable and the grid is fixed on the lower part. It is filled with salt water which gives us a small density contrast between the two fluids inside the two tanks. The Atwood number A was varied in our experiments from 0.9 × 10 -3 to 1.34 × 10 -2 . Having very small Atwood numbers allows for slower dynamics of the instability and thus a better visualisation. Indeed, when the self similar regime is reached, the mixing zone width grows as a function of the Atwood number (Eq. ( 3.1) of Chap. 3). Thus, the smaller the Atwood number, the slower the growth of the mixing zone.

The upper tank is filled through aspiration created using pressurised air at the top of the tank. Thus the water is sucked inside the tank through the grid as shown in Fig. 4.5, which is a photo taken during the filling of the upper tank. We can see on it the upper tank, the grid fixed at its bottom and the tube allowing for the sucking of the salt water using pressurised air. As seen previously, the surface tension at the grid and the slight depression created at the top of the tank are sufficient to retain the body of salt water inside, even as we move it through air to put it into contact with the bottom reservoir.

Three grids were chosen for our experiments and are given in Tab. A. The grid mesh size d and thread diameter l are shown in Fig. 4.6 and are used to define the porosity ϕ of the grid as written in Eq. (4.10). The effects of the grid porosity, mesh size and thread diameter are studied and discussed in Section 4.3.

Grid n

• d [mm] l [mm] ϕ 1 1 0.

Porosity of the grid

ϕ = d 2 (d + l) 2 .
(4.10)

Note that the mesh size used is much smaller than the critical wavelength given by the linear stability analysis, thus these grids can indeed be used to stabilise the system. Once the contact between the fluids is made, the brine falls inside the fresh water in a Rayleigh-Taylor instability like manner (Fig. 4.7). The linear stability analysis is also performed (Sec. 4.1.3) on this brinefresh water interface to better understand how the grid may help in getting more reproducible experiments. The visualisation of this phenomenon is made using the Planar Laser-Induced Fluorescence (PLIF) method [START_REF] Van Cruyningen | Quantitative imaging of concentration by planar laser-induced fluorescence[END_REF]). As shown in Fig. 4.4b, this method requires the projection of a continuous laser sheet onto the flow field in the bottom cylinder in order to induce fluorescence. Indeed, as the fresh water is mixed with a fluorescent material, it emits light at a very specific wavelength when illuminated by the laser sheet. In the experiments, a 5 W continuous Spectra Physics-Millenia eV CW DPSS laser emitting at 532 nm (green) and sulforhodamine as the fluorescent material emitting at 552 nm (yellow) were used.

Experimental work 1 cm

The light emitted by the fluorescent material is recorded at 30 frames per second using a camera and a filter for the laser light. Only the bottom cylinder can be recorded through this process as it is where the laser sheet is created and the grid prevents it from being present in the upper cylinder. This process gives images as shown in Fig. 4.7 with the fresh water in light gray, the heavier fluid in black and the grid visible at the top of the image. From these images, the concentration fields inside the bottom cylinder can be computed at any given time. Indeed, the concentration in salt of a parcel of fluid can easily be deduced from the colour of the corresponding pixel: the darker the pixel, the more salted the fluid is. These concentration fields can the be used to compute other quantities like mean concentration profiles or mixing zone widths.

Ethanol is also added to the mixture of fresh water and sulforhodamine in order to match the optical index of the salt water and thus have a clear picture throughout the experiment.

Linear stability analysis of the brine-fresh water interface: experimental grids

When the two cylindrical tanks are put into contact, a salt water-fresh water interface appears. The linear stability analysis of this interface is made in order to compare the growth rate Ω(k) to the grid mesh size. The viscous and diffusive effects are taken into account in this study while the surface tension is disregarded. Two cases are considered: the case of a sharp interface and the case of a diffuse interface.

The linear stability theory for Rayleigh-Taylor instability onset states that the growth rate can be written as Eq. (4.11) (see [START_REF] Duff | Effects of diffusion on interface instability between gases[END_REF] and [START_REF] Morgan | Parametric investigation of the transition to turbulence in Rayleigh-Taylor mixing[END_REF]) when both the viscous and diffusive effects are considered. (4.11) where ν is the kinematic viscosity, D is the diffusion coefficient and ψ = 1 + 0.375kδ and δ is the diffusion thickness of the interface.

Ω = AG 0 k ψ + ν 2 k 4 1/2 -(ν + D)k 2 ,
This value of ψ was found by [START_REF] Duff | Effects of diffusion on interface instability between gases[END_REF] for A = 0.5, which is much bigger than the Atwood number used in our experiments. As shown in [START_REF] Duff | Effects of diffusion on interface instability between gases[END_REF], the smaller the Atwood number, the bigger the slope of ψ as a function of kδ. This would shift the maximum of Ω towards the smaller wavenumbers k slightly. However, ψ is only slightly dependent on the Atwood number so the shift would be negligible.

Sharp interface

In the case of a sharp interface, ψ = 1 and thus the dispersion relation is written:

Ω = AG 0 k + ν 2 k 4 1/2 -(ν + D)k 2 . (4.12)
This dispersion relation is plotted against the wavenumber k in Fig. 4.8 along with the fastest growing wavenumber k 0 in dashed blue line and the critical wavenumber k c in red dashed line. The parameters taken here are those of a typical experiment: A = 0.0037, G 0 = 9.81 m.s -2 , ν = 1 × 10 -6 m 2 .s -1 and D taken so that the Schmidt number (Sc = ν/D) is around 700. largest mesh size in dashed green line (k g ). The presence of a grid at the interface between the two fluids prevents the development of any perturbation having a wavelength bigger than the mesh size. Thus, looking at the experimental grids, the largest perturbation wavelength allowed by the grids is 2 mm giving a wavenumber k g of around 3141 m -1 . This wavenumber is much smaller than the critical wavenumber k c , which is around 23130 m -1 , above which the growth rate becomes negative and diffusion is expected to stabilize the growth. This means that the grid does not stabilize the instability between fresh water and salt water. The other two experimental grids also have mesh grids giving wavenumbers smaller than the critical wavenumber. However, we notice that k g is larger than k 0 : the most unstable wavelength cannot develop in our system due to the presence of the grid. We would need a grid of mesh size 3.8 mm at the smallest if we wanted to allow for the fastest growing wavelength to develop. Not having the most unstable wavelength available means that the Rayleigh-Taylor instability early dynamics will be impacted. This could then alter the long term dynamics. In the absence of the grid, the selected wavelength will be the most unstable one given by 2π/k 0 whereas when the grid is present, it is very likely that the selected wavelength will be very close to the one given by the mesh size (2π/(2 1/2 k g ) because it is not a 2D mode) as it is the one with the biggest growth rate. The advantage however, is that the grid imposes an initial condition at small scales.

For bigger Atwood numbers, the maximum growth rate would be shifted towards the higher wavenumbers. In order to verify if none of the Atwood numbers tested in our experiments allow the development of the fastest growing wavelength, we plot in Fig. 4.9a the growth rate as a dark gray area and corresponding k 0 and k c for the smallest and largest Atwood numbers used in our experiments. On this figure is also shown k g , the wavenumber of Grid3 which is the grid with the largest mesh size.

Notice that even the largest Atwood number, A = 0.0134, is not big enough to have k 0 larger than k g which would allow for the development of the most unstable wavelength even with the Chapter 4. The Rayleigh-Taylor instability through a grid: Influence of the grid grid present. Indeed, the Atwood number from which k 0 becomes larger than k g is around 0.026 for Grid3 and 0.2 for Grid1 (see Fig. 4.9b).

As a conclusion, it is clear than in every experiments made with these grids, the wavelength that developed was always the one selected by the holes of the grid. This can help to better control the initial conditions and thus to make more reproducible experiments.

Diffuse interface

The case of the diffuse interface is studied to understand the effect of the thickness of the interface on the growth rate of the instability. Figure 4.10 shows the values the growth rate takes when varying the interface thickness from the sharp interface (δ = 0 mm) to a diffusion thickness δ = 5 mm. We do not know the diffusion thickness of the interface in the experiments, it is a data that can only be speculated on. Here we consider that the interface thickness is not thicker than 5 mm.

We notice that having a thick diffuse interface compared to a sharp one tends to shift the growth rate towards smaller values and k 0 and k c towards the smaller wavenumbers. Thus, the most unstable wavenumber is shifted away from the smallest wavenumber allowed by the experimental grids k g . Meaning that it would take an even higher Atwood number to allow for the most unstable wavelength to develop despite the presence of the grid than the one found in the case of a sharp interface. Moreover, there are less wavenumbers available between the minimum wavenumber of the grids and the critical wavenumber, thus less wavelengths can perturb our system.

The experiments made and their observations are described in the next section and it was indeed observed that filaments of the size of the holes developed and then mixed very rapidly with each other. Evolution of the viscous growth rate Ω(k) with the increase of the diffusion thickness of the interface δ (dark area) from 0 to 5 mm. The evolution of the most unstable wavenumber k 0 (blue area), the critical wavenumber k c (red area) are also given and the grid wavenumber k g (dashed green line) are also given. The parameters are A = 0.0037, G 0 = 9.81 m.s -2 , ν = 1 × 10 -6 m 2 .s -1 , D = ν/700.

Experimental work

Observations

We made several experiments varying the grid and the Atwood number. These experiments and their parameters are given in Tab. B. The series of runs are organized depending on the grid used.

Series N In 74% of experiments, we observed the appearance of an ascending central jet at the very beginning of the instability, in which the lighter fluid goes up in the center of the cylinder, or at least in the center of the plane formed by the laser sheet, while the heavier fluid goes down on the sides of the cylinder (see Figure 4.11). In Tab. B, the velocity for the descending lateral jets Chapter 4. The Rayleigh-Taylor instability through a grid: Influence of the grid (V DLJ ) of heavier fluid is given in the last column. It seems that a larger grid porosity leads to a larger V DLJ . Indeed, for A ∼ 0.0037 and ϕ = 0.41, the mean V DLJ is around 1.29 cm.s -1 while for the same Atwood number but ϕ = 0.48, the mean descending lateral jet velocity is around 1.85 cm.s -1 . The cause of the apparition of this central jet is investigated in Chapter 5.

• d [mm] ϕ A ρ h [kg.m -3 ] ρ l [kg.m -3 ] ACJ V DLJ [cm.s -1 ]
In the other 26% experiments, a Rayleigh-Taylor like instability was observed as the one seen in Fig. 4.7. Indeed, even if thin filaments of the size of the holes develop, as the grid's threads are extremely small, the mixing happens very close to the grid and mushroom like structures form. In order to assess whether we recover the dynamics of the Rayleigh-Taylor instability, some of these experiments are further studied.

Experimental work

Experimental data

The homogeneous experiments are studied in this section with the purpose of focusing only on the influence of the grid on the late time dynamics of the Rayleigh-Taylor instability. From the recorded images, the mean horizontal concentration profile (Eq. (4.14)) can be computed. These profiles are shown in Fig. 4.12 at times t ∼ 2 s and t ∼ 4 s for various experiments (EXPb2, EXPb5, EXPb9, EXPb13 and EXPb14 of Tab. B) having the same grid but different Atwood numbers. Only the values of ⟨C⟩ XY below z = 0 can be accessed as only the bottom cylinder is recorded through the PLIF method.

The mean horizontal concentration is written:

Horizontal mean concentration ⟨C⟩ XY (Z, t) = 1 S S C(X, Y, Z, t) dXdY, (4.13)
where S is the surface of the XY plane.

As we only have access to the XZ plane in the experiments, the experimental mean horizontal concentration is written:

⟨C⟩ XY (z, t) = 1 D cyl R cyl -R cyl C(X, Y, z, t) dX, (4.14)
A small jump in the concentration can be seen at the initial interface position (z = 0 mm) which we attribute to the presence of the grid. In order to make this jump more visible, the mean horizontal concentration ⟨C⟩ XY is plotted against the height z for EXPb5 on Fig. Chapter 4. The Rayleigh-Taylor instability through a grid: Influence of the grid it tends to become smaller with time. Indeed, the amplitude of the jump on Fig. 4.13a is much larger at t ∼ 2 s than at t ∼ 4 s. This means that the jump in the mean concentration profile tends to disappear as the instability develops. The reason behind the presence of this jump is investigated further through direct numerical simulations (DNS) in Section 4.3.

Using the mean horizontal concentration profiles, the mixing zone width given in Eq. ( 4.15) along with the growth rate α of Eq. ( 4.16) are computed for the same set of homogeneous experiments (EXPb2, EXPb5, EXPb9, EXPb13 and EXPb14) and are plotted in Fig. 4.14. A moving average was applied to L(t) in order to smooth the fluctuations.

Mixing zone size and Growth rate

• Mixing zone size [START_REF] Andrews | A simple experiment to investigate two-dimensional mixing by Rayleigh-Taylor instability[END_REF]):

L(t) = 6 fluid ⟨C⟩ XY (z, t) (1 -⟨C⟩ XY (z, t)) dz, (4.15)
• Growth rate:

α(t) = L2 8AG 0 L . (4.16)
Note the presence of a small artifact at the very beginning of the experiment with the mixing zone width decreasing from an initial non zero value. This artifact probably comes from the lack of definition of the origin point. Indeed, z is defined as zero at the top of the image but not at the specific position where the instability begins.

We would expect that as the Atwood number increases, the mixing zone width becomes larger. Indeed, a larger Atwood number means a greater buoyancy force. This is what we see between experiments EXPb5, EXPb9 and experiments EXPb13, EXPb14. However, this statement is true when the initial conditions are strictly the same. Otherwise, the transient regimes between two experiments can be very different with L(t) growing very rapidly initially and not being caught up by the other L(t) with larger A. This is probably what happens for EXPb2 which has the smallest Atwood number but a larger mixing zone width than EXPb5 and 9.

Direct Numerical Simulations

Looking at the growth rate α(t), a plateau seems to be reached at least for the experiments EXPb2, 5 and 9. This indicates that the self-similar regime of the Rayleigh-Taylor instability is reached. This is less evident for the last two experiments considered which could be due to their obvious vertical confinement. Indeed, the height of the cylindrical tank is of 30 cm which is very close to the width of the mixing zone at the end of these experiments. Moreover, if the instability is vertically confined then it will no longer develop proportionally to t 2 , thus the desired plateau of α will no longer exist. The values of the growth rate are obtained by taking the mean value of α between t = 3 s and t = 7 s (gray area in Fig. 4.14b). Note that the values are around 0.034 which is an expected value for an experimental Rayleigh-Taylor instability [START_REF] Dimonte | A comparative study of the turbulent Rayleigh-Taylor instability using high-resolution three-dimensional numerical simulations: The Alpha-Group collaboration[END_REF]).

Experimental data is lacking to truly deduce whether the grid's presence influences the long term dynamics of the instability and thus whether the Rayleigh-Taylor instability is recovered. Hence, in order to gather more data, direct numerical simulations of a Rayleigh-Taylor-like instability through a grid are performed in Section 4.3. The next Section (4.2) is dedicated to the description of the direct numerical simulations.

Direct Numerical Simulations

As stated previously, several direct numerical simulations (DNS) are performed in this study in order to investigate the dynamics of the instability when a grid is present at its interface. This section is devoted to the description of the initial conditions and the penalisation method.

Description of numerical simulations

These simulations solve numerically the Navier-Stokes equations under the Boussinesq approximation and are performed in a triply periodic cubic box of length L box = 2π using the code Stratospec already described in Chapter 2 [START_REF] Viciconte | Sudden diffusion of turbulent mixing layers in weakly coupled plasmas under compression[END_REF], [START_REF] Viciconte | Sudden diffusion of turbulent mixing layers in weakly coupled plasmas under compression[END_REF], [START_REF] Briard | Harmonic to subharmonic transition of the Faraday instability in miscible fluids[END_REF] and [START_REF] Briard | The turbulent Faraday instability in miscible fluids[END_REF]), which is based on a pseudo-spectral collocation method with two-third rule dealiasing. The time advancement is realized through a third-order low-storage strong-stability-preserving Runge-Kutta scheme, with implicit viscous terms. The simulations use a N p X × N p Y × N p Z grid box with a pencil-decomposition on N p cores and an Atwood number of 0.003, N p i being the number of points in one direction.

Initial Conditions

The initial conditions on the velocity and concentration fields of the Stratospec code used for this study are the following:

Initial conditions: velocity field

• Velocity field: U(x, t = 0) = 0. (4.17)

• Concentration field (hyperbolic tangent concentration profile):

C(x, t = 0) = 1 2 1 + tanh z -S(x, y) σ + C inh (x), (4.18)
where σ defines the initial width of the interface δ = 3σ for δ = 6 ⟨C⟩ XY (1 -⟨C⟩ XY )dz.

The initial condition on the concentration field is given by Eq. 4.18: its steepness is controlled by σ (the smaller σ the steeper the tangent). This profile is shown in Fig. 4.16 in blue line.

The function S(x, y) gives the form of the interface and of the initial perturbations at this interface. In the case of a planar interface, S(x, y) is the perturbation surface which is defined from a 2D perturbation spectrum of the form:

E(k) = E 0 k k peak s i exp - s i 2 k k peak 2 = S k Ŝ(k) Ŝ * (k)d 2 k ′ , (4.19)
where E 0 is a renormalisation factor (linked to the amplitude of the deformation of the interface) and s i is the infrared slope of the spectrum, taken in the simulations as 1. In Eq. ( 4.19), k peak defines the wavenumber of the maximum of the spectrum E. The amplitude of the perturbations at the interface is thus given by S ′ 2 with S ′ 2 =

+∞ 0 E(k)dk and is set at 3×10 -4 m (0.15% of the height of the box) in our simulations. Moreover, in Eq. (4.19), S k is the sphere of radius k, where k = |k| and Ŝ is the Fourier transform of S ′ . In Fig. 4.15, the 2D perturbation spectrum is represented with k peak shown with a red dashed line.

In order to ensure the periodicity in the vertical inhomogeneous direction, we need to change the concentration field from its heavy value to its light value at the vertical boundaries. This is done with C inh (x) which is given by: C .20) This contribution to the initial concentration profile is shown as the red dashed line in Fig. 4.16.

E k k peak s i exp -s i 2 k k peak 2 ∼ k s i with s i = 1 k peak
z [m] C(z, t = 0) = 1 2 1 + tanh z σ C(z, t = 0) = 1 2 1 + tanh z σ + C inh (z)
inh (x) = - 1 2 tanh z + L box /2 σ + tanh z -L box /2 σ . ( 4

The penalization method

Penalised layers with no-slip boundary conditions are added for various reasons, explained below, inside the cubic box. The goal of this section is to describe the principle of the penalisation method, also called immersed boundary method, and the different penalised layers added.

Chapter 4. The Rayleigh-Taylor instability through a grid: Influence of the grid

Equations

The penalisation method is described in Chap.

Momentum equation with penalisation

∂U ∂t + (U • ∇)U = -∇P + ν∇ 2 U - 1 η u χ(U -U Pen ), (4.21)
where P is the reduced pressure p/ρ 0 , where ρ 0 ≡ (ρ h + ρ l )/2, η u is the velocity penalisation parameter, χ is the mask function and U Pen is the velocity field imposed inside the penalisation layer.

The mask function is defined as a door function which is equal to 1 inside the penalisation layer and 0 inside the fluid. As explained in [START_REF] Briard | The turbulent Faraday instability in miscible fluids[END_REF], this mask function is filtered in spectral space to avoid Gibbs oscillations with exp(-C s (k/N p ) 2 ), where N p is the number of points in one direction and C s is the filter parameter. Thus, the additional force only exists inside the penalisation layer. The velocity penalisation parameter imposes a constraint on the time step ∆t ≤ η u : the smaller η u , the smaller the time step has to be to achieve convergence.

Advection-diffusion equation with penalisation

∂C ∂t + ((1 -χ)U • ∇)C + χ (U Pen • ∇C) = ∇ • ((D(1 -χ) + η c χ)∇C) , (4.22)
with D the diffusion coefficient and η c the scalar penalisation parameter which can be interpreted as an effective diffusivity inside the penalised layers.

The scalar penalisation parameter decreases the diffusion of the concentration inside the penalised layer: as it does not impose any constraint on the time step, it can be chosen very small in order to completely freeze that diffusion. The term χ (U Pen • ∇C) is important if a nonzero velocity field is imposed inside the penalised layer as it transports the artificial scalar field created inside the layer.

In order for the created boundary to correspond to a solid wall, small values of the pair of parameters (η u , η c ) must be taken.

Cylindrical configuration

The Rayleigh-Taylor experiments through a grid were made inside cylindrical tanks whereas the code used for the simulations only simulates a cubic domain. Thus, in order to reproduce the experimental geometry in the simulations, penalisation layers are added on the sides of the cubic box to create a cylinder, thus breaking the periodicity of the simulations. Penalisation layers are also added at the top and bottom of the domain to ensure the periodicity in this direction and prevent any fluid exchange at this boundary. Two configurations are thus considered: the Large Cylinder (LC) configuration and the Thin Cylinder (TC) configuration. The TC reproduces the exact dimensions of the experimental tanks with a diameter of 2R cyl = 0.194 m

Direct Numerical Simulations

η u [s] η c [m 2 .s -1 ] U Pen LC L box /2 -R cyl = 0.01 L box /2 -H cyl /2 = 0.01 1 × 10 -2 1.2 × 10 -10 0 TC L box /2 -R cyl = 0.23 L box /2 -H cyl /2 = 0.03 1 × 10 -2 1.2 × 10 -10 0
Table C: Parameters of the penalisation layers giving the two cylindrical configurations, Large Cylinder (Fig. 4.17b) and Thin Cylinder (Fig. 4.17a): width of the penalised layer on the sides of the cylinder L Pen,S , width of the penalised layer on the top and bottom of the cylinder L Pen,TB , velocity penalisation parameter η u , scalar penalisation parameter η c and the velocity field in the penalisation layers U Pen .

and a height of H cyl = 0.3 m per cylinder. However, this configuration penalises a large part of the computational domain. The LC represents a zoom in on the interface with a height of only H cyl = R cyl = 0.097 m per tank. With this configuration, the structures are much better resolved at the interface, but it will be more easily vertically confined. Both of these configurations are shown in Fig. 4.17 where the penalised areas are presented in light gray.

The penalisation layers used to create these cylindrical tanks in the simulations need to mimick solid walls, meaning that the smallest possible pair (η u , η c ) should be chosen, without constraining too much the time step. Moreover, a null velocity field has to be imposed inside these layers (U Pen = 0). The parameters taken for the penalised layers of each configurations are given in Tab. C with L Pen,S and L Pen,TB being the width of the layers on both sides and on the top and bottom of the cylinder respectively.

Validation of the cylindrical configuration

The cylindrical configuration can be achieved inside our simulation's cubic box by adding penalisation layers on the sides (see Fig. 4.17). Thus, it is important to assess whether adding these walls has an impact on the instability. Two configurations are compared, the cubic one and the Large Cylinder one (see Tab. D), with the exact same initial conditions: A = 0.003, ν = 0.6 × 10 -6 m 2 .s -1 , D = ν, η u = 1 × 10 -2 s, η c = 1.2 × 10 -10 m 2 .s -1 and k peak = 585.90 m -1 . Simulation DNSCubVSCyl1 is made inside a triply periodic cubic box whereas simulation DNSCub-VSCyl2 is made inside the cylindrical geometry.

Series

N In Fig. 4.19, the horizontal mean concentration profiles, ⟨C⟩ XY (see Eq. ( 4.13)), are plotted at t ∼ 5 s (see Fig. 4.19a) and the time evolution of the mixing zone size, L(t) (replacing ⟨C⟩ θr by ⟨C⟩ XY in Eq. (4.15)), for both simulations DNSCubVSCyl (see Fig. 4.19b).

As expected from the instantaneous concentration fields, very little difference is visible either on the horizontal mean concentration or on the mixing zone width. The mixing zone width grows slightly less in the cylindrical configuration as we see in Fig. 4.19b when comparing the cubic configuration (blue) and cylindrical configurations (orange). This is probably due to the no-slip boundary condition imposed by the penalised layers in the cylindrical configuration which makes the fluid near the walls adhere to them a little.

We would expect greater differences between the two configurations in the case of a very thin cylinder as the important quantity is the horizontal integral scale. Indeed, as was investigated in [START_REF] Lawrie | Turbulent diffusion in tall tubes. I. models for Rayleigh-Taylor instability[END_REF], when the instability develops inside a laterally confined domain, the evolution of the mixing zone width is constrained by the lateral walls rendering the scaling in L(t) ∝ t 2 obsolete in the long term. The instability does not seem horizontally confined by the lateral walls in the configuration considered here. As a conclusion, using penalisation layers to recreate a cylinder tends to slightly slow down the dynamics of the instability. Since this effect is very small it should not prevent us from using the cylindrical geometry to conduct our study.

Vertical confinement

As said previously, two cylindrical configurations (Thin Cylinder and Large Cylinder) can be considered for the following numerical studies. The first is a better reproduction of the experiment's geometry whose main disadvantage is to lose a large part of the computational domain to penalised layers; while the second is more of a close up on the interface between the two fluids whose main disadvantage is to be vertically less elongated. Thus, it is important to assess whether the configuration impacts the instability's development: the goal of this section is to evaluate whether the vertical confinement issues appear in the LC configuration.

Two simulations, whose parameters are given in Tab. E, were made to compare these two configurations. The physical parameters are the same for the two simulations, the difference is the size of the computational domain that is changed from L box ∼ 0.21 m for the LC configuration, to L box ∼ 0.66 m for the TC configuration. However, the physical diameter of the cylinder stays the same: D cyl = 0.194 m. The simulations parameters are once again: A = 0.003,

ν = 0.6 × 10 -6 m 2 .s -1 , D = ν, η u = 1 × 10 -2 s and η c = 1.2 × 10 -10 m 2 .s -1 .
The instantaneous concentration field at different times and the slices along the XZ plane are given in Fig. 4.20 for both simulations. The slices of DNSTC are cropped to match the size of the slices of DNSLC.

In order to ensure that no confinement effect appears in simulation DNSLC compared to DNSTC, the mixing zone width as a function of time is plotted for both simulations on Fig. 4.21. Only a very slight difference can be seen between the two simulations. Indeed, DNSLC grows slightly faster than DNSTC. This difference is due to the better resolution, and therefore lesser diffusion, of the simulation in the Large Cylinder configuration. This shows that the Large Cylinder configuration is not vertically confined as it would mean a slower mixing zone growth. Thus, both configurations yield very similar results. As the Large Cylinder configuration is better resolved than the Thin Cylinder one, it is a better candidate for the study of the early times of the instability. However, as it would end up being vertically confined, for the study of the late time dynamics, either the Thin Cylinder configuration, or the Large Cylinder configuration with a computational domain twice as high as shown in Fig. 4.22, would be better choices.

Direct Numerical

Grid definition

As the goal of this study is to assess the impact of the grid on the instability's dynamics, a grid is implemented at the interface between the two fluids using the penalisation method. Hence, a new mask function, called χ G , is created to define the position and width of the threads of a squared threaded grid. In one-dimension, it is defined by either Eq. ( 4.23), giving a centered grid with a thread at the position (x, y) = (0, 0), or by Eq (4.24), giving a non-centered grid (no thread at (x, y) = (0, 0)).

• Centered grid in 1D:

   χ G (x) = 1 if |x| + l/2 -(d + l) × F |x|+l/2 d+l < l, χ G (x) = 0 else (4.23)
where F(x) is the floor function.

• Non-centered grid in 1D:

   χ G (x) = 1 if |x| + d/2 -(d + l) × F |x|+d/2 d+l > d, χ G (x) = 0 else. (4.24)
This definition of the mask function gives χ G = 1 inside the threads of width l and χ G = 0 inside the holes of width d. Thus, inside the threads the flow is penalised with (η u , η c ) = (1 × 10 -2 s, 1.2 × 10 -10 m 2 .s -1 ). Figure 4.23 shows the value of the mask function χ G in 1D (along the x-axis) with the holes being the white spaces between each thread in blue. The centered grid (Fig. 4.23a) is n the left and the non-centered grid (Fig. 4.23b) is on the right. In this figure, the thread diameter is l = 8 mm and the mesh size is d = 16 mm, meaning that the holes are twice the size of the threads, giving a porosity of ϕ = 0.44.

A relevant question is: can we implement the same grid as the ones used in the experiments?

Direct Numerical Simulations

Indeed, we know that the experimental grids have very thin threads and holes, thus we would need a very large resolution to make it possible (a minimum of 2383 3 points to have 10 points in each threads for Grid3). However, this is not only a question of resolution but also of whether the instability would develop in such small numerical holes. To answer this question, the linear stability analysis for the simulations is made in the following section.

Linear stability analysis of the brine-fresh water interface: numerical grid

The linear stability analysis for this brine-fresh water interface taking the viscosity and diffusive effects into account has already been made for the experiments. However, as the diffusion coefficient is much larger in the simulations (D = ν) than in the experiments (D = ν/700), a different result for the dispersion relation of Eq. ( 4.11) is expected. Considering a simulation whose parameters are A = 0.003, G 0 = 9.81 m.s -2 , ν = 1.2 × 10 -6 m 2 .s -1 , D = ν and δ = 3 mm, the growth rate Ω(k) is shown in Fig. 4.24. On the same figure is plotted the growth rate Ω EXP considering a sharp interface and a diffusion coefficient D = ν/700. These correspond to the conditions of the experiments, thus the differences between the simulations Ω and the experiments growth rate Ω EXP are evidenced. The most unstable wavenumber k 0 , the critical wavenumber k c and the minimum wavenumber allowed by the experimental grids k g are also plotted on this figure.

Note that the critical wavenumber is around 2050 m -1 , thus the smallest wavelength that can be used to destabilize the system is λ c = 3 mm which is larger than the biggest mesh size available in the experimental grids. Hence, we cannot numerically implement one of the experimental grids in our simulations due to the much bigger diffusion in the numerical simulations 4.24: Growth rate Ω(k) (Eq. 4.11) (black continuous line) for simulation parameters (ν = D = 1.2 × 10 -6 m 2 .s -1 and δ = 3 mm) compared to the growth rate Ω EXP (k) for experiment parameters (ν = 700D = 1 × 10 -6 m 2 .s -1 and δ = 0 mm). The red lines mark the critical wavenumber k c and the blue lines give the fastest growing wavenumber k 0 . The smallest wavenumber allowed by the experimental grids k g is given in green dashed line (with λ g = 2 mm for Grid3). The common parameters are A = 0.003, G 0 = 9.81 m.s -2 .

Chapter 4. The Rayleigh-Taylor instability through a grid: Influence of the grid which would stabilise the instability.

Grid convergence

In this section, the mesh convergence of the simulations with grid is evaluated. A good way to gauge the mesh convergence is to compare the mixing zone width for simulations with different spatial resolution, meaning different number of points in the X and Z directions. If the mesh convergence is achieved, then the mixing zone width will be the same no matter the resolution. Indeed, the averaging process involved in the computation of L(t) should counteract the differences brought by the modification of the initial conditions when the resolution is changed. However, implementing a grid means there are much less points inside the fluid around the interface which could bring a greater sensitivity to the initial conditions despite the averaging process.

In order to study the sensitivity to the initial conditions and more precisely the effect of the random initialisation, four simulations are made in each resolution, N p X = N p Z = 512 and N p X = N p Z = 1024. Due to time constraint and a desire for quick results, these simulations were made with a 2D geometry. As such, the differences between 2D and 3D simulations need not be forgotten (Chapter 3). The seed of the random function s is changed in each simulation of each resolution as shown in Tab. F. The parameters of the simulations are taken as: A = 0.003, ν = 1.2 × 10 -6 m 2 .s -1 , D = ν and λ peak = 2π/k peak ∼ 11 mm. The grid's mesh size and thread diameter are d = 34 mm and l = 17 mm respectively, giving a porosity ϕ = 0.44.

Series

N Table F: Label (series and numbers), number of points in the X and Z directions N p X and N p Z respectively and seed of the random function s. The quasi-2D simulations have an Atwood number A = 0.003, a viscosity ν = 1.2 ×10 -6 m 2 .s -1 , a diffusion coefficient D = ν and a λ peak = 11 mm. The grid in these simulations have a thread diameter of l = 17 mm and a mesh size of d = 34 mm giving a porosity ϕ = 0.44.

The mixing zone widths of all these simulations are plotted against time in Fig. 4.25a. Note that none of the mixing zone width overlap, they are all very different from one another, meaning that L(t) in this configuration strongly depends on the seed s of the random function used to initialise the simulation. Hence, the number of fluid points being insufficient, the result of the simulation is very sensitive to the initial condition. This dependence on the initial conditions can be very clearly seen on the 2D instantaneous concentration fields shown in Fig. 4.26. Indeed, very different structures can be seen from one seed to another or from one spatial resolution to the other. The range of mixing zone widths covered by the DNSConvA and by the DNSConvB are given in Fig. 4.25b along with their ensemble mean defined by: 1 N ∑ N L(t, s), (4.25) where N is the number of simulations with the same spatial resolution (the same N p,XZ ) and L(t, s) is the seed-dependent mixing zone width. These simulations serve another purpose: to see whether the lack of fluid points would be compensated by the number of simulations. Indeed, the idea is that statistical convergence Chapter 4. The Rayleigh-Taylor instability through a grid: Influence of the grid could be achieved by ensemble averages. However, these ensemble means do not overlap either, meaning that to reach statistical convergence, a much higher number of simulations would probably be needed. For practical reasons of cost and computation time, however, using such a greater number of simulations in order to achieve mesh convergence cannot be done. We expect to have fewer issues achieving convergence in a three-dimensional geometry as more points are available in the fluid. Thus, a few three-dimensional simulations are performed (see Tab. G) with the same grid as the ones presented above (l = 17 mm and d = 34 mm). Moreover, simulations with a smaller perturbation wavelength are made to understand its importance on the mesh convergence. The simulations are made inside a cubic box, periodic along X and Y. Some of the simulations only had threads along the X direction (DNSConv3DA and DNSConv3DC) while others had a real grid, meaning threads along the X and Y directions (DNSConv3DB and DNSConv3DD).

Direct Numerical Simulations

Series N • N p X = N p Y = N p Z λ peak [mm]

Direct Numerical Simulations

The slices along the XZ plane of the instantaneous concentration fields are given at t ∼ 4 s for these simulations in Fig. 4.27. Note that for the smallest initial wavelength, the instability is less developed as expected from the linear stability analysis and Fig. 4.28. Indeed, we can see on Fig. 4.28 that even if both wavelengths are greater than the most unstable wavelength, the smallest one is much closer to the critical wavelength, and thus has a much smaller growth rate than the other. The mixing zone widths are plotted for these simulations in Fig. 4.29. The simulations with threads only in the X direction are closer to mesh convergence than the ones with the threads along both X and Y directions. This is in accordance with the assumption that the lack of points in the fluid prevents the convergence. A slight difference can be seen between the simulations of initial wavelength λ = 5 mm and λ = 11 mm. Indeed, the simulations with λ = 5 mm have a slightly better mesh convergence than their counterparts with λ = 11 mm.

As expected, the mixing zone widths between different resolutions are closer for the 3D simulations than for the 2D ones, thus the convergence is slightly better due to the increase in the number of points. However, it is clear that achieving mesh convergence with a grid at the interface is very difficult as even the 3D simulations are sensitive to the initial conditions.

A way to achieve this mesh convergence is to drastically increase the number of points inside the fluid by decreasing the number of threads implemented. Thus, two periodic 2D simulations are made with only one thread of width 17 mm in the middle of the domain with N p X = N p Z = 512 and N p X = N p 2 = 1024.

The 2D instantaneous concentration fields are shown in Fig. 4.30 for t ∼ 4 s and the mixing zone widths are given in Fig. 4.31.

Note that even if the concentration fields are not exactly the same due to the change of spatial resolution giving slightly different initial conditions, the mixing zone widths overlap perfectly. Thus, mesh convergence is reached in these simulations. As a conclusion, adding a grid at the interface decreases the number of points in the fluid domain and renders the simulations very sensitive to the initial conditions, which means that the mesh convergence is very difficult to achieve. Increasing the porosity of the grid by increasing the mesh size and thus having less threads will help improve the mesh convergence.

Numerical study of the influence of the grid

The main goal is to understand whether having a grid at the interface between the two fluids impacts the dynamics of the Rayleigh-Taylor instability. In this section, numerical simulations of the Rayleigh-Taylor instability through a grid are presented and compared to the experimental results. All of the simulations presented in this section are performed inside a three-dimensional periodic cubic domain.

Horizontal mean concentration profiles and model

In this section, different method of computing the mean horizontal concentration are detailed and a model is created with the purpose of understanding the cause of the jump at the interface of the mean concentration profiles observed in the experiments. With the addition of the grid at the interface, three ways of computing the mean concentration profile exist. The first, giving ⟨C⟩ XY,tot , is to consider every points even those inside the threads of the grid (dashed white in Fig 4.32a). The second, giving ⟨C⟩ XY,wg (wg meaning without grid), is to consider all of the points except the ones inside the threads of the grid (dashed white in Fig 4.32b). Lastly, the third method, giving ⟨C⟩ XY,wgc (wgc meaning without the grid's columns), considers only the points inside the columns delimited by the holes (dashed white in Fig 4.32c). This is summarized in Fig. 4.32 showing in dashed white the points taken into account in the computation of the horizontal mean concentration for each of the three methods. The simulation's parameters considered here are: 512 3 points, A = 0.003, ν = 1.2 × 10 -6 m 2 .s -1 , D = ν, η u = 1 × 10 -4 s, η c = 1.2 × 10 -10 m 2 .s -1 and λ peak = 11 mm.

As expected, each method gives a different horizontal mean concentration profile. Indeed, taking into account the points inside the columns delimited by the threads like in the first two methods presented, clearly shifts the value of the mean concentration towards 0 below the interface and towards 1 above the interface. This means that more pure, unmixed, fluid is taken into consideration in the averaging process, hinting that the instability only develops inside the columns delimited by the holes of the grid. These differences can be seen in Fig. 4.33, showing the three horizontal mean concentration profiles at time t ∼ 5 s for a three-dimensional periodic cubic simulation with a grid mesh size of d = 10 mm and a thread diameter of l = 10 mm.

To prove that the instability does, indeed, develop mainly in the columns delimited by the holes of the grid, at least in the early dynamics, a model can be made. This model is based on the fact that the columns above and below the threads of the grid are filled with pure unmixed fluids, which allows to compute ⟨C⟩ XY,tot knowing ⟨C⟩ XY,wgc .

For this model, an elementary configuration that is a square of side length d + l composed of a hole (square of side d) and two threads of width l (Fig. 4.34) is considered. The surface of 4.3. Numerical study of the influence of the grid the hole is thus d 2 , the surface of the threads is (d + l) 2 -d 2 meaning that the surface of the elementary configuration is (d + l) 2 . Inside the hole, the mean horizontal concentration is ⟨C⟩ XY,wgc while inside the threads, the mean horizontal concentration would be the one taking only into account the columns delimited by the threads ⟨C⟩ XY,g . The value of ⟨C⟩ XY,g is computed based on the hypothesis that, above the thread is a column of pure dense fluid, ⟨C⟩ XY,g = 1, and below the thread is a column of pure lighter fluid, ⟨C⟩ XY,g = 0. In the width of the interface δ, the value of ⟨C⟩ XY,g is considered to grow linearly from 0 to 1. This principle is summarised in Eq. (4.26). Taking into account the proportion of each horizontal mean concentration in the elementary configuration to form the total horizontal mean concentration ⟨C⟩ XY,tot , we get:

       ⟨C⟩ XY,g = 0 for z < -δ 2 , ⟨C⟩ XY,g = z δ + 1 2 for -δ 2 ≤ z ≤ δ 2 , ⟨C⟩ XY,g = 1 for z > δ 2 . ( 4
⟨C⟩ XY,tot = ⟨C⟩ XY,wgc d 2 (d + l) 2 + ⟨C⟩ XY,g (d + l) 2 -d 2 (d + l) 2 . ( 4.27)
Thus:

Model: Total horizontal mean concentration

⟨C⟩ XY,tot = ⟨C⟩ XY,wgc ϕ + ⟨C⟩ XY,g (1 -ϕ), (4.28)
Figure 4.35 shows the results of this model at time t ∼ 5 s in red line while the blue and green lines are the computed horizontal mean concentrations ⟨C⟩ XY,wgc and ⟨C⟩ XY,tot . A very good agreement can be observed between the model and simulation. As time grows, however, the assumption of pure fluids inside the columns delimited by the threads of the grid become less true and thus the model and the computed means tend to depart from each other.

The model does show that the shift seen in the horizontal mean concentrations between the two first methods and the third is due to the presence of pure unmixed fluid trapped above and In green, the total horizontal mean concentration profile ⟨C⟩ XY,tot computed by taking into account every points. In red, the total horizontal mean concentration profile ⟨C⟩ XY,tot computed with the model Eq. (4.28).

below the threads of the grid. Well above or below the threads, the fluid will end up mixing and will not stay pure which will cause the horizontal mean concentrations ⟨C⟩ XY,tot and ⟨C⟩ XY,wg to become closer to ⟨C⟩ XY,wgc . However, the fluid immediately surrounding the threads stays trapped and unmixed, meaning that these mean concentrations will never be exactly the same. Moreover, from the model it can be deduced that when the grid is more porous (when ϕ is high), ⟨C⟩ XY,wgc is closer to ⟨C⟩ XY,tot , which can be very easily explained as less pure fluid is trapped.

The jump in the concentration profile when considering ⟨C⟩ XY,wgc can also be explained by the no-slip boundary condition at the grid. Indeed, the no-slip boundary means that some pure fluid is trapped right next to the threads inside each hole.

Adding a grid at the interface: Importance of k 0 and effect of the porosity

In order to better understand the effect of a grid placed at the interface between two fluids in the development of the Rayleigh-Taylor instability, several three-dimensional simulations are made inside a periodic domain with an Atwood number of A = 0.003, a kinematic viscosity of ν = 1.2 × 10 -6 m 2 .s -1 and a diffusion coefficient of D = ν. The computational domain is a rectangular one with a box that is twice as high as it is large as shown in gray in Fig. 4.22 (L box × L box × L box ), thus the dimensions of the tank, and thus of the volume of fluid, are

L box × L box × 2H cyl ∼ 0.21 × 0.21 × 0.41 m.
This vertically elongated domain was chosen here so that the instability would not be confined too rapidly compared to a cubic box. Moreover, the simulations were not made inside a cylindrical tank as it was shown that the geometry should not impact the dynamics of the instability too much. Hence, the geometry was chosen to maximize the number of points inside the fluid and thus being as close to convergence as possible.

Numerical study of the influence of the grid

Several grids with different porosities, mesh sizes and thread diameters are tested along with several initial wavelength for the perturbation. These varying parameters are given for each simulation in Tab. H. The notation 13/5 in Tab. H refers to the fact that two wavelengths are tested in the DNSGD cases: λ peak = 13 mm and λ peak = 5 mm.

Series N • λ peak [mm] d [mm] l [mm] ϕ DNSGA 1 13 ∅ ∅ ∅ 2 5 ∅ ∅ ∅ DNSGB
A = 0.003, a kinematic viscosity of ν = 1.2 × 10 -6 m 2 .s -1 and a diffusion coefficient of D = ν.
The aim is to compare a Rayleigh-Taylor simulation without grid (DNSGA) and a Rayleigh-Taylor simulation with grid. However, several parameters are to be taken into account when implementing a grid, like its porosity, its mesh size and especially the value of the mesh size compared to the most unstable wavelength λ 0 given by the linear stability analysis.

Simulations with k g < k 0

Firstly, the DNSGB simulations are compared to the DNSGA2 one. These simulations correspond to the grid with the biggest mesh size implemented and to the case that is furthest from the experiments as its mesh size is larger than the most unstable wavelength, d > λ 0 . This is shown in Fig. 4.36 where the growth rate is plotted against the horizontal wavenumber.

In this figure are given in red both the most unstable wavenumber k 0 and the critical wavenumber k c after which the growth rates become negative and the instability is stabilised by diffusion. The wavenumber corresponding to the mesh size of the grid k g = 2π/d is given in blue and the initial perturbation wavenumber k peak is in green. The hatched gray area represents all of the wavenumbers that cannot develop due to the presence of the grid.

As can be seen in Fig. 4.36, the most unstable wavelength is allowed to develop despite the presence of the grid and thus is the wavelength that is expected to develop in the DNSGB simulations.

The slices along the XZ plane of the instantaneous concentration fields are given at time t ∼ 9 s for DNSGA2 and all of the DNSGB simulations in Fig. 4.37.

It is to be noted that the grid giving a porosity of ϕ = 0.83 of DNSGB3 in image 4.37b has very few points inside each thread, which is why it is barely visible and its resolution could very well be questioned. It was verified, however, that no fluid passes through the threads and that

k [m -1 ] -1 0 1 2 Ω [s -1 ] Growth rate Ω k 0 = 470.00 m -1 k c = 1400.00 m -1 k g = 306.77 m -1 k peak,1 = 1318.26 m -1 Figure 4
.36: Growth rate Ω as a function of the wavenumber k for the parameters of DNSGB (Tab. H) with, in red, the most unstable wavenumber k 0 and the critical wavenumber k c , in blue, k g = 2π/d the wavenumber corresponding to the mesh size and, in green, the initial perturbation wavenumber k peak .

The gray hatched area corresponds to all of the wavenumber that are not allowed to develop by the grid. the velocity stays null inside them.

Figure 4.37 shows that, as expected, the more porous the grid, the closer we get to the Rayleigh-Taylor case without grid DNSGA2. Another noticeable feature is that, when the porosity decreases, as the holes are kept of the same width, the threads grow in size (and decrease in number) and the instability is more and more trapped inside the columns delimited by the holes as was understood from the mean concentration model. Moreover, as the instability cannot de-4.3. Numerical study of the influence of the grid velop sideways for the lower porosity, it develops more vertically and thus less mixing happens. That is only true at the beginning of the instability, as we see here that the turbulent structures do not stay trapped inside these holes delimited channels for the entirety of the simulation. Thus, the horizontal scale becomes larger than the grid mesh size in the long term.

The presence of these elongated structures are corroborated by both the integral scale L u and the dimensionality parameter sin 2 (γ) [START_REF] Gréa | The rapid acceleration model and the growth rate of a turbulent mixing zone induced by Rayleigh-Taylor instability[END_REF]). A note of warning, these parameters depend on spectral quantities in which the contribution of the grid cannot be easily removed and can bring large fluctuations. This effect is in all probability negligible for small threads but not for large threads. Thus, L u and sin 2 (γ) aim only at bringing a qualitative confirmation of the observations made. These quantities are plotted in Fig. 4.38a and Fig. 4.38b respectively for each of the simulations presented in Fig. 4.37. The integral scale is defined as:

L u = E 3/2 cin ε (4.29)
where E cin is the kinetic energy and ε its dissipation. The integral scale shows the extension of the structures, and as expected it becomes smaller as the porosity is decreased with for example a factor two between the classical Rayleigh-Taylor case and for ϕ = 0.44. However, it does not say whether the structures become thinner but more vertically elongated or wider but flatter.

In order to access this information, we look at the dimensionality parameter sin 2 (γ) (defined in [START_REF] Gréa | The rapid acceleration model and the growth rate of a turbulent mixing zone induced by Rayleigh-Taylor instability[END_REF]). Indeed, for isotropic structures, this parameter is equal to 2/3, if the structures are more vertically elongated, it tends towards 1, whereas decreases below 2/3 for flat horizontal structures. For the classical Rayleigh-Taylor case, we expect sin 2 (γ) ∼ 0.7 [START_REF] Gréa | The rapid acceleration model and the growth rate of a turbulent mixing zone induced by Rayleigh-Taylor instability[END_REF]) and this is what is observed in Fig. 4.38b. When the porosity is decreased however, we observe larger dimensionality parameters meaning that the structures are indeed more vertically elongated. These confirm the observations made on the instantaneous concentration fields. Moreover, from Fig. 4.38a, we see that the presence of the grid at the interface delays the transition to the selfsimilar regime as L u is supposed to be proportional to t 2 in that regime.

To quantify further these effects, the mean horizontal concentration profiles, the mixing zone widths and the growth rates are plotted on Fig. 4.39. The mean horizontal concentration profiles plotted in Figs. 4.39a and 4.39b were computed through the third method ⟨C⟩ XY,wgc described in Sec. 4.3.1 where only the points inside the columns delimited by the holes of the grid are taken into account. The case DNSGA2, in blue line, has no grid and is noted ϕ = 1 as it could be seen as a completely porous case. These profiles are given at times t ∼ 5 s and t ∼ 9 s to understand the evolution of the profiles with time. Note that for each simulation with a grid, a jump exists in the mean concentration profile, whatever the porosity may be. This is as was seen in the mean concentration profiles of the experiments studied earlier. On these profiles an inversion of slope can be seen just before the interface (z = 0 m), both above and below it, that was not seen in the experiments. This is believed to come from a portion of pure fluid that is trapped by the thread in its vicinity, probably in part due to the cubic geometry of the threads in the simulations. These pouches of trapped fluid are larger when the porosity is decreased. It is clearly visible on Fig. 4.39b that this behavior in the mean horizontal concentration is present for the smallest porosity, ϕ = 0.44, which is due to the larger thread diameter and thus larger pouches. Moreover, looking at the width of the jump at 4.3. Numerical study of the influence of the grid the interface, it becomes smaller with time which is especially visible for the largest porosity, ϕ = 0.83. However, its evolution is slower than in the experiments, and it does not disappear entirely.

The mixing zone width in Fig. 4.39c is computed with two different methods. Indeed, the classical method, the integral method given in Eq. (4.30), does not take into consideration small filaments of fluids that passes through the holes of the grid and thus the resulting mixing zone width is very sensitive to the porosity and holes size, and has a tendency to be smaller than what could be seen from the instantaneous concentration fields.

Mixing zone width : Integral method

L int (t) = 6 ⟨C⟩ XY,wgc (1 -⟨C⟩ XY,wgc )dz.

(4.30)

A way to take these filaments into account is to compute the mixing zone width with the profile method given by Eq. 4.31. This method is based on an arbitrary threshold on the value of the mean horizontal concentration giving the z-coordinates at which the mixing zone commences and finishes in the mean concentration profiles.

Mixing zone width : Profile method

L z (t) = z(⟨C⟩ XY,wgc = 1 -ϵ) -z(⟨C⟩ XY,wgc = ϵ), (4.31)
ϵ = 0.03 being the value of the arbitrary threshold.

The drawback of this second method is its sensitivity to small variations in the concentration profiles (spike in the dashed green curve in Fig. 4.39c) and to the value of the threshold.

In Fig. 4.39c, as expected, L int decreases as the porosity increases and gives much smaller values than L z . The value of the arbitrary threshold was chosen so that the mixing zone widths would be similar for the case without grid with both methods. Thus the value chosen is ϵ = 0.03. Note that when considering the mixing zone width L z they are all very close together, with a tendency to be slightly larger than the one without grid. This is in all probability due to the larger vertical extension of the instability in the presence of a grid discussed earlier. With a less porous grid, the structures are thinner and mix less which means that they grow faster than in a classical Rayleigh-Taylor case, for which the mixing happens sooner.

As there are two ways of computing the mixing zone width, two α can also be computed:

α int coming from the integral method and α z coming from the profile method (Fig. 4.39d). The growth rate considered in the experiments is the one computed from the integral method, α int .

In the case where the two mixing zone widths are similar, α int and α z should also be similar. This can be seen for the DNSGA2 (ϕ = 1) simulation which shows a growth rate around 0.02 for both methods. However, the profile method gives way more fluctuations on the value of the growth rate compared to the integral method. These fluctuations can be very intense (green dashed curve around 6 s) and they make it difficult to understand whether a plateau has been reached or not. As expected from the mixing zone width, α int becomes smaller as the porosity decreases. This stays partially true for α z except for the simulation DNSGB1 with ϕ = 0.44 whose growth rate is larger than the one for the simulation DNSGB2 with a porosity of 0.6 and as large as the one for the porosity 0.83 (DNSGB3). Note that except for this occurrence, both methods of computing α give close values of growth rate despite their mixing zone width being quite different. We have seen in Fig. 4.38a that the simulations with the less porous grids have not reached the self-similar regime. This means that the growth rate α in Fig. 4.39d is not a measure of the self-similar growth of the mixing zone with grid but only an evaluation of the growth rate of the mixing zone.

The mean concentration chosen to compute the mixing zone width and its growth rate was the one which does not take into account the columns of fluids above and below each thread of the grid, ⟨C⟩ XY,wgc . In the case where ⟨C⟩ XY,tot is used instead, a higher effect of the grid would be expected with larger mixing zone width discrepancies between each simulation, and also larger growth rate for the simulations with grids because of the larger density contrast. This behavior is shown in Fig. 4.40.

Indeed, the mixing zone widths for the simulations with grid part ways much sooner (around 1 s) than in Fig. 4.39c (around 6 s) which give a much smaller mixing zone at t = 10 s for simulations DNSGB1 and DNSGB2. Moreover, the growth rates are indeed a bit larger than the ones computed from Fig. 4.39d. However, the differences are not too important and do not change the evolution of the growth rate with the porosity from one simulation to the next. Hence, both mean horizontal concentration can be considered to study the impact of the grid on the dynamics. One only takes into account the regions of the computational domain where the instability is able to develop from the start, thus becoming slightly closer from the Rayleigh-Taylor without grid case. The other takes every point into consideration and is thus closer to what can be computed from the experiments.

As a conclusion, the grid imposes thinner and more elongated structures as the porosity is decreased. As a consequence, the method used to compute the mixing zone width becomes critical if we want to detect these structures. The same goes for the growth rate α as they are even more sensitive to how L is computed. In order to understand the importance of letting the most unstable wavelength λ 0 develop or not, the DNSGD simulations are compared to the DNSGA simulations. The grid that is closest to the grid used in the experiments, from the linear stability analysis point of view, is the grid used in DNSGD. Indeed, this grid keeps the same distance between k g and k c and between k g and k 0 than the experimental Grid3. However, this grid has a mesh size that is smaller than the most unstable wavelentgh, d < λ 0 . This means that the most unstable wavelength cannot develop and the wavelength with the largest growth rate is d: thus a perturbation of wavelength d is expected to develop. This is shown in Fig. 4.41 giving the growth rate as a function of the wavenumber. The most unstable wavenumber k 0 and the critical wavenumber k c are plotted in red, while the grid wavenumber k g = 2π/d and the initial perturbation wavenumbers k peak,1 and k peak,2 are given in blue, orange and green respectively. The gray hatched area gives the wavenumbers that cannot be destabilised due to the grid. The initial perturbation wavenumber k peak,1 was especially chosen to be equal to the most unstable wavenumber. Thus, note that both k peak,1 and k 0 are inside the gray hatched area and thus they cannot develop.

The simulations DNSGA and DNSGD were both made with two different initial perturbation wavelengths, λ peak,1 = 13 mm and λ peak,2 = 5 mm. The instantaneous concentration fields of Fig. 4.42 show at t ∼ 9 s the simulations done with λ peak,1 while the ones of Fig. 4.43 show at t ∼ 9 s the simulations done with λ peak,2 .

Observe that, as expected, the perturbation that develops in simulations DNSGD is of the size of the holes whatever the initial perturbation wavelength may be. Indeed, for both initial perturbation wavelengths, filaments of the hole's diameter develop. This development of the 0 250 500 750 1000 1250 1500

k [m -1 ] -1 0 1 2 Ω [s -1 ] Growth rate Ω k 0 = 470.00 m -1 k c = 1400.00 m -1 k g = 766.94 m -1 k peak,1 = 468.72 m -1 k peak,2 = 1318.26 m -1 Figure 4
.41: Growth rate Ω as a function of the wavenumber k for the parameters of DNSGB (Tab. H) with, in red, the most unstable wavenumber k 0 and the critical wavenumber k c , in blue, k g = 2π/d the wavenumber corresponding to the mesh size and, in green, the initial perturbation wavenumbers k peak,1 and k peak,2 . The gray hatched area corresponds to all of the wavenumbers that are not allowed to develop by the grid. instability is very different from the previous case when we had k g < k 0 . These filaments that develop here provoke very little mixing, which becomes increasingly obvious as the grid's porosity decreases. It is in appearance very far from the known Rayleigh-Taylor instability development which can be seen in both Fig. 4.42a and Fig. 4.43a.

Looking at both DNSGA simulations, the initial perturbation wavelength impacts the dynamics of the instability. Indeed, around the same time, Fig. 4.42a shows a mixing zone width that is more developed than in Fig. 4.43a. This is an expected behavior as more energy is given in DNSGA2 to a mode that has a much smaller growth rate Ω than k 0 and which is quite close to the critical wavenumber k c . This phenomenon can be seen also when comparing Fig. 4.44a and Fig. 4.44b giving the mixing zone widths computed from ⟨C⟩ XY,tot for λ peak,1 and λ peak,2 respectively. Indeed, we see on Fig. 4.44a that the mixing zone width for ϕ = 1 reaches much higher values than on Fig. 4.44b.

This statement is less true for the simulations with grid: they have very similar mixing zone widths despite the change of initial perturbation wavelengths. Once again however, the integral method gives much smaller mixing zone widths for the simulations with grids than the profile method, with a net decrease of L(t) with the porosity.

The growth rate of these 8 simulations are given in Fig. 4.45. Note that the initial perturbation wavelength λ peak has a great influence on the growth rate α both for the simulations with and without grid. Indeed, even if the grid mesh size prevents large wavelengths from developing, the development of the instability is still sensitive to the smaller wavelengths present in the system and their initial energy. As such, the transient regime is shorter with higher values of α reached for λ peak,1 than for λ peak,2 . In the self-similar regime, α for the classical Rayleigh-Taylor is slightly larger for the largest initial perturbation wavelength whereas the contrary is seen for the simulations with grid. The growth rate of the simulations width grid is rather larger when using a smaller initial perturbation wavelength.

Thus, even if it is not clearly visible on the instantaneous concentration fields or on the mixing zone widths, the value of the initial perturbation wavelength has an effect on the instability's dynamics even with a grid that seems to select its own wavelength.

The mean horizontal concentration profiles ⟨C⟩ XY,wgc (z) for the λ peak,1 DNSGA and DNSGD simulations are plotted in Fig. 4.46.

The jumps is clearly still present even if its amplitude tend to become smaller with time. At t ∼ 5 s, the slope inversion is only visible for the simulation with ϕ = 0.83. That is due to the lack of mixing in the simulations with lower porosity: there are only columns of pure fluid ascending and descending. However, as time passes, the mean concentration profile reverts to a linear profile thanks to the turbulence (before and after the jump) for the ϕ = 0.83 case while the slope inversion becomes visible for the other two simulations with grid as mixing occurs (t ∼ 9 s). The filamentation of the instability slowing down the mixing, it also tend to slow down the apparition of the slope inversion in the mean concentration profile. As a conclusion, there seems to be an effect of whether d is larger or smaller than the most unstable wavelength λ 0 . In the case where d < λ 0 (the case of the experiments), only filaments of width d develop. The less porous the grid and the more filamentous the instability becomes. These thinner structures mix badly thus delaying the mixing in our system compared to a classical Rayleigh-Taylor instability.

Effect of the mesh size

Two other mesh sizes were tested with DNSGC and DNSGE that were equal and smaller than λ 0 to understand whether the filamentation of the instability does come from the fact that d < λ 0 . We notice that DNSGC1, having a grid mesh size that only just allows the most unstable wavelength to develop, is less filamented than the other two simulations but more than the DNSB1 simulation shown in Fig. 4.37d. Hence, we can assume that for d > λ 0 , several wavelengths can develop in each hole bringing more mixing and less filamentation whereas for d < λ 0 , only one wavelength can develop in each hole, mixing less and filamenting more. The case d = λ 0 is the transition case between these two limits. In the experiments, d < λ 0 , thus the experimental grid only allows for perturbations of the size of the holes to develop. There is a selection of the initial perturbation wavelength by the grid at the beginning of the instability. Moreover, as long as the grid stays the same, the same wavelength will be selected thus it should help the experiments to be more reproducible.

Numerical study of the influence of the grid

Moreover, the grid hole's width has a great impact on the size of the vertical extension. This is expected as the smaller the hole, the smaller the wavelength of the developing instability and thus, the smaller the growth rate Ω. The smaller Ω, the more time the instability will take to develop. This is why the vertical extension of the structures in the simulations with the smallest mesh size d (DNSGE1 in Fig. 4.47c) is much smaller than the one in the simulation with the largest mesh size (DNSC1 in Fig. 4.47a).

As a conclusion, the most important grid parameter in terms of the structures shapes, sizes and selected wavelength is the mesh size and where it is situated compared to the most unstable wavelength. The porosity of the grid also has an effect and that is to decrease the growth rate of the instability in the self-similar regime and thus, to get further away from the wanted Rayleigh-Taylor instability. The grid, no matter its porosity or mesh size, always seems to have an impact on the mean concentration profiles by creating a jump at the interface and on the instability's Chapter 4. The Rayleigh-Taylor instability through a grid: Influence of the grid dynamics.

Conclusion

Conclusion

A new type of miscible Rayleigh-Taylor experiment was introduced in this chapter with a grid positioned at the interface between the two fluids. The goal of the grid was to stabilize the heavier fluid when the container is filled and to let the instability develop once the tanks containing the lighter and heavier fluids are in contact. Thus, it would give an initial condition that does not create any shear at the interface; hence, the grid could make the initial condition more controllable. However, for this new type of experiment to be a good candidate to investigate the Rayleigh-Taylor instability, we need to make sure that the grid does not influence the dynamics of the instability by its continued presence during the experiments.

Unfortunately, most of the experiments made showed an ascending jet around the center of the cylinder. Thus, only a few number of experiments could be studied to answer this question. From these, however, it was found that the horizontal mean concentration profiles presented a jump at the interface that should not appear in a 'classical' Rayleigh-Taylor instability. Hence, the grid has an effect but it does not mean that it necessarily has an important impact on the dynamics. Indeed, the mixing zone width and growth rate show that something resembling a self-similar regime could be reached with values of the growth rates around 0.034 for some of the experiments, which is an usual value for experiments of the Rayleigh-Taylor instability.

In order to get more data to study the effect of the grid, direct numerical simulations (DNS) using the Stratospec code were performed. Using the penalisation method, different configurations were tested and a grid was implemented at the interface between the two fluids. One of the result of this study is that using the penalisation method to simulate a cylinder inside the cubic computational domain does not impact much the development of the instability despite the no-slip boundaries imposed. Thus, it is perfectly possible to reproduce the geometry of the experiments.

More importantly, it was found that the grid does impact the development of the instability. Indeed, one important result is the striking effect the grid mesh size has on the structures shapes and sizes. The key parameter here is whether the ratio d/λ 0 (λ 0 being the most unstable wavelength given by the linear stability analysis) is larger or smaller than one. For a ratio larger than one, the structures that develop are very close to the ones observed in a Rayleigh-Taylor instability. Whereas, for a ratio smaller than one, the structures become thin filaments of the size of the holes that mix very poorly. In the experiments, this ratio is smaller than one: hence the grid imposes an initial condition of the size of the holes. Contrarily to the numerical simulations, the experimental grid's holes and threads are much smaller and thus mixing occurs very close to the grid which makes the structures appear much less filamentous.

Moreover, as soon as a grid is implemented in the simulations, no matter the value of the ratio d/λ 0 , the jump in the concentration profiles appears. As seen in the experiments, this jump tends to shrink with time: so, as the instability develops, the presence of the grid has less and less impact. A self-similar regime seems to be reached in the numerical simulations as well, but this time with much smaller values of the growth rate than for the 'classical' Rayleigh-Taylor instability. This value is greatly impacted by the porosity of the grid and the more porous the grid, the higher the value of the growth rate α. However, with the apparition of the thin filamentous structures, the quantities of interest, like the growth rate and the mixing zone width, become very sensitive to how they are measured. This adds dispersion to the results meaning Chapter 4. The Rayleigh-Taylor instability through a grid: Influence of the grid we should be careful when interpreting them.

In order to better answer the question of whether the grid impacts the long term dynamics of the instability in the experiments, we could perform simulations with a higher resolution to have finer grids and decrease the diffusion coefficient D imposed in the simulations in order to allow more wavelengths to develop inside each holes and thus recover an earlier mixing.

The Rayleigh-Taylor instability through a grid: Central jet Introduction

In the previous Chapter, a new type of Rayleigh-Taylor experiment, involving a grid at the interface between the two fluids, was described. This previous chapter was centered around the influence of the grid on the dynamics of the instability. However, in 74% of the experiments described, an inhomogeneity was observed: namely an ascending central jet. Thus, the goal of this chapter is to find the cause of the appearance of that jet through direct numerical simulations.

To this aim, the mean azimuthal concentration, given by Eq. (5.1), will often be studied.

⟨C⟩ θZ (r, t) = 1 4πH cyl 2π 0 H cyl -H cyl C(r, θ, Z, t) dθdZ,
(5.1)

with H cyl the height of one of the cylindrical tank, Z is the vertical coordinate, (r =

√ X 2 + Y 2 , θ)
are the radius and angle in the cylindrical coordinate system. An ascending central jet of lighter fluid means, in the plane of the observations, two descending lateral jets of heavier fluid on the sides of the lower cylindrical tank, the only one visible in the experiments. Thus, at the position of these jets (r = R cyl ), the mean azimuthal concentration should increase while it should decrease in the center of the cylinder (r = 0). In the case where no jets appear, the mean azimuthal concentration should be uniform and around 0.5 along the radius of the cylinder for the simulations in which both the upper and lower tanks are taken in consideration.

In the experiments, only a 2D plane for the bottom tank is visible. Thus, the mean concentration is defined as:

⟨C⟩ Z (X, t) = 1 H cyl 0 -H cyl C(X, Z, t) dZ. (5.2)
Hence, if no descending jet is present then ⟨C⟩ Z (X, t) = 0. An example of this behaviour is shown in Fig. 5.1 for EXPa7 (Tab. B) at time t ∼ 7 s. The ascending jet in this experiment happened on the right side of the cylinder and thus there is only one descending lateral jet on the left side.

The figure presents both the cropped instantaneous concentration field of EXPa7 at time t ∼ 7 s and the mean concentration profile as a function of the radius r. As expected the highest value of ⟨C⟩ Z is reached at the position of the descending jet. A dark part of the grid is visible at the top of the image which increases artificially the value of the mean concentration. This is why it never quite reaches 0 even in the right part of the image where almost no heavier fluid can be seen. This chapter articulates as follows: first of all, various intuitive hypothesis are tested to reproduce the ascending central jet, namely modifying the initial conditions. Unfortunately, this is not enough to mimic the phenomenology of the experiments, meaning that more complex mechanisms are at play. Hence, a simplified 'aspiration model' is derived, consisting in a timeevolving forcing velocity field which results from the pressure jump between the two fluids when contact is made.

The cause of the apparition of the central jet

The first intuition one can have is that particular initial conditions may be at the origin of the ascending central jet. Following this idea, direct numerical simulations, using the code Stratospec, are performed, and three different initial conditions are tested, where we vary:

1. the wavelength for the initial perturbation λ peak 2. the curvature of the interface C

3. the amount of time it takes to establish the contact between the two fluids t o

All of the simulations presented in this section are made using 512 3 points in a cylindrical geometry using 2R cyl = 0.194 m and 2H cyl = 0.194 m (Large Cylinder configuration presented in Chapter 4) with an Atwood number A = 0.003, a kinematic viscosity ν = 0.6 × 10 -6 m 2 .s -1 and a diffusion coefficient D = ν. For these simulations, no grid is implemented at the interface.

Wavelength of the initial perturbation?

The ascending central jet visible in the experiment denotes the breaking of the lateral homogeneity and a breaking of the symmetry with respect to the z = 0 plane. Thus, the goal is to know whether initialising our simulations with a perturbation of large wavelength, and thus possibly confining laterally the instability, could help reproduce this phenomenon. Hence, we varied our initial wavenumber k peak given in Eq. (4.19) (and thus our initial wavelength λ peak = 2π/k peak ), in order to find out its importance on the dynamics of the instability.

Changing the wavelength of the initial condition modifies the growth rate, defined in Eq. (4.11), of the early linear regime of the Rayleigh-Taylor instability as shown in Fig. 5.2. It is to be noticed that among the four simulations made, two of them possess a k peak smaller than k 0 , the fastest growing wavenumber, and two of them have a higher one. All of them have a positive growth rate, thus the diffusion cannot stabilize the growth of the instability. However, the diffusion will affect more the simulations whose k peak are bigger than k 0 . We expect the growth rates for the four k peak to be: Ω DNSλ1 = 2.6 s -1 , Ω DNSλ2 = 3.2 s -1 , Ω DNSλ3 = 2.7 s -1 and Ω DNSλ4 = 1.7 s -1 . So the one which should develop the most rapidly is DNSλ2 with k peak = 585.90 m -1 . In all of the DNSλ simulations, we have k peak ≫ k D cyl , where k D cyl = 2π/D cyl is the wavenumber associated with the diameter of the cylinder. Thus, we do not expect these simulations to be laterally confined.

The instantaneous 3D concentration fields and the vertical slices along the XZ plane are shown on Fig. 5.3 for simulations DNSλ1 and DNSλ3 with k peak = 292.95 m -1 and k peak = 1318.26 m -1 respectively, giving an initial perturbation of wavelength 0.021 and 0.005 m. Note that none of these images show an ascending central jet like the one observed in the experiments (Fig. 4.11 in Chap. 4). Indeed, we only see the development of a homogeneous Rayleigh-Taylor instability with different structure sizes and mixing zone widths. The symmetry is preserved between the lighter and heavier fluids.

In order to confirm that no symmetry with respect to the z = 0 plane is broken in the simu- Thus, for two perfectly mixed fluids, the mixing parameter Θ tends towards 1. As demonstrated in [START_REF] Gréa | The rapid acceleration model and the growth rate of a turbulent mixing zone induced by Rayleigh-Taylor instability[END_REF], both α and Θ are related, and in particular, the larger Θ the smaller α. Indeed, when the two fluids mix, their effective density contrast is reduced and so is the production term of kinetic energy, causing the mixing zone to grow less rapidly. Looking at Figs. 5.4c and 5.4d, the initial phase is indicated by Θ tending towards 0, no mixing occurs. Then, the transition to turbulence is seen by the rapid increase (decrease) of Θ (α), the growth rate slows down and mixing occurs.

We see in the insert of Fig. 5.4b and in the early times of Fig. 5.4d that simulation DNSλ2 has the fastest growth rate and thus the fastest growth of the mixing zone in the beginning of the instability, as was expected from the linear stability analysis. However, this does not last very long as the simulation DNSλ1 quickly catches up during the transient regime and has the biggest growth rate before turbulence settles (around 0.16). When the self similar regime is reached, we see that the mixing zone width is still larger for DNSλ1 due to the lead it took during the transient state. But its growth rate is the smallest of the four simulations, consistently with the largest Θ. This indicates that the bigger structures of the simulation DNSλ1 are less sensitive to diffusion, and thus they develop more rapidly in the beginning of the instability and mix less. However, when the transition to turbulence occurs, their development is slowed down meaning that more mixing occurs compared to smaller structures. Note that the mixing parameter reaches a plateau in the self similar regime around 0.8 for all of the simulations which is an expected value for a well developed mixing layer.

The growth rates α of all four simulations reach also a plateau in Fig. 5.4d, when the self similar state is reached with expected values of around 0.02 [START_REF] Dimonte | A comparative study of the turbulent Rayleigh-Taylor instability using high-resolution three-dimensional numerical simulations: The Alpha-Group collaboration[END_REF]). This means that the mixing zone width grows like t 2 , which shows that our simulations are not confined, as expected, neither vertically nor laterally. Indeed, it was shown that when the instability is confined, the mixing zone width time dependency in the self similar regime is affected. Focusing on the lateral confinement, it was found (Debacq et al. [2003] and [START_REF] Lawrie | Turbulent diffusion in tall tubes. I. models for Rayleigh-Taylor instability[END_REF]), that when the amplitudes of the Rayleigh-Taylor instability are large compared to the width of the tank containing the fluids, the self-similar evolution of the mixing zone width becomes proportional to t 2/5 . Another simulation was made, in the configuration Thin Cylinder (see Fig. 4.20a in Chap. 4) where the height of the cylinder is 60 cm and the diameter is 19.4 cm, with an initial wavelength of 13 cm. This value is very close to the diameter of the cylinder so a lateral confinement is expected (k peak ∼ k D cyl ). The 3D instantaneous concentration fields and 2D vertical slices along the XZ plane are shown at time t ∼ 4 s in Fig. 5.5. Large structures are visible, around half the diameter of the cylinder, but no central ascending jet is present. A slice at a later time is also shown to convince ourselves that even later no central jet appears.

Figure 5.4 gives the mixing zone width L, the growth rate α and the mixing parameter Θ as a function of time for the DNSTC simulation. Note that the initial phase, is much longer in this case than the other simulations presented. Thus, the instability grows faster during the transient regime. The cylinder being twice as high as in the previous simulations, the instability is not vertically confined. As for the lateral confinement, the beginning of a self-similar regime seems to be reached at the very end of the simulation. It is difficult to conclude on this point but, whether it is laterally confined or not, no central jet can be seen. As the wavelengths chosen in this study are already much larger than the ones allowed by the experimental grids (Tab. A and Sec. 4.1.3), it was not considered pertinent to try even larger wavelengths. Moreover, the ascending central jet develops rapidly at the beginning of the instability in the experiments. Such a rapid development could not be observed by decreasing even further k peak as the growth 5.1. The cause of the apparition of the central jet rate would become too small.

To conclude, we can safely deduce that the wavelength of the initial perturbation is not the cause of the appearance of the central jet and that the choice of the k peak influences mainly the growth of the mixing zone width of the instability. For the following simulations k peak will be chosen to be 585.90 m -1 .

Curvature of the interface?

In our experiments, the grids are not perfectly planar, they have a small curvature (of the order of a millimeter) towards the bottom cylinder which allows for the evacuation of the air initially trapped underneath the grid. Indeed, if the grid was perfectly planar, the air would stay trapped and would form a bubble between the grid and the lighter fluid compromising the experiment. This small curvature also creates a large scale perturbation, even larger than the ones tested above by varying k peak .

This curvature of the grid implies that the salt water-fresh water interface is slightly curved. The question is: can this curvature have an effect on the homogeneity of the dynamics of the mixing zone? Several simulations are performed, DNSCurv given in Tab. B, in which a curved interface of 1 or 5 mm towards the lighter fluid is implemented as the initial condition. The curvature of the grid is taken into account in these simulations but the grid is not implemented.

Series N • C int [mm] DNSCurv 1 1 2 5
Table B: Label (series and numbers) and the interface curvature C int in millimeter are given for the two simulations considered.

In order to curve the initial interface, the function S(x, y) (see Sec. 4.2.2 in Chap. 4) is redefined as: S(x, y) = ξ(x, y) + a init Rand,(5.4) where Rand is a random function giving perturbations with amplitudes centered around 0 and of maximum values ±a init , and the ξ(x, y) is the paraboloid function defined in Eq. (5.5) As a conclusion, the curvature of the grid does have an effect on the development of the instability but that effect is opposite to the one observed in the experiments. Thus, this curvature cannot be the reason for the appearance of the central jet. In the following, a curvature of 1 mm will sometimes be imposed to assess whether a combination of various effects may trigger or inhibit the ascending central jet.

ξ(x, y) = x b 2 + y b 2 + z 0 , ( 5

The time of contact between the two fluids?

As explained previously, the grid is slightly curved to allow for the air to be evacuated. This means that the contact between the two fluids is not instantaneous everywhere. Indeed, this contact is first made near the center of the grid and then propagates towards the extremities, getting rid of the trapped air. This process, represented in Fig. 5.8b, takes about 0.5 s. If it is not rapid enough it could lead to an asymmetry in time in the development of the instability between the center of the cylinder and the sides.

To determine the importance of the contact time upon the dynamics, various simulations are performed, whose parameters are given in Tab. C. To mimic this contact delay, a small cylindrical solid wall of height L pen,PW = 0.01 m is introduced by penalisation method, with the penalisation parameters η u = 1 × 10 -2 s and η c = 1.2 × 10 -10 m 2 .s -1 , at the center of the domain between the two fluids (Fig. 5.8a). This wall opens from the center at speed v PW so that the contact radius is:

R(t) = v PW t.
(5.7)

For simplicity, a linear dependence is chosen here, but other profiles could be considered (an exponential profile is briefly addressed in Appendix A). At time t o , later referred to as the aperture time, one has R(t o ) = R cyl . Once the wall is open, no grid is present at the interface between the two fluids. Note that two different aperture velocities (v PW = 0.194 m.s -1 and v PW = 0.0485 m.s -1 ) and two different types of interface (planar or curved) were tested.

The instantaneous concentration fields and vertical slices along the XZ plane are given at different times in Fig. 5.9 for DNSContact3 and DNSContact4 (Tab. C). Observe that, in Fig. 5.9d at t = 2.0 s, the Rayleigh-Taylor instability has had the time to start developing at the center Chapter 5. The Rayleigh-Taylor instability through a grid: Central jet of the cylinder before the wall is completely opened (before t o ). An asymmetry is thus created between the center of the cylinder and the sides. However, this effect is eventually overcome by the effect of the curvature of the interface, where the lighter fluid will go up on the sides and the heavier fluid will go down in the center. For an aperture time of 0.5 s, the instability does not have the time to be triggered before the wall is completely opened (Fig. 5.9b at time t = 1.0 s) and thus no asymmetry occurs and the effect of the curvature is the only one visible. Hence, a competition exists between the asymmetry created by the time of contact between the two fluids and the one created by the curvature of the interface. Be that as it may, this effect is only seen for longer aperture time than the ones observed in the experiments. It also appears that it is not a good candidate to explain the apparition of a central jet.

Series N • k peak [m -1 ] C int [mm] t o [s] v PW [m.s -1 ]
In order to confirm the observations made with the instantaneous concentration fields, the azimuthal mean concentration profiles (Fig. DNSContact3 and 4. This large peak, consistent with the ⟨C⟩ θZ profile, is due to the delay of growth between the sides and the center of the cylinder. The mixing parameters on Fig. 5.10c show significant dispersion between the simulations with a curved interface and the one with a planar interface. Indeed, the simulations with a curved interface seem much more sensitive to the aperture time, as can be seen on the red curve representing DNSContact4 which decreases almost linearly until the aperture time is reached. They also reach a much smaller value, around 0.6, while the simulations with a planar interface reach a value closer to 0.75. This seems to indicate that the curvature actually slows down the mixing and the dynamics of the instability.

v 2 z θZ [m 2 .s -1 ] t o = 0.5 s, C int = 0 mm t o = 2 s, C int = 0 mm t o = 0.5 s, C int = 1 mm t o = 2 s, C int = 1 mm (
Hence, having a bigger aperture time than the time of development of the instability coupled with a small curvature at the interface really changes the dynamics of the instability. However, the aperture times during our experiments were closer to 0.5 s than 2 s, so this effect should not be seen experimentally.. It will be taken in consideration in the rest of this work. In any way, it is clear that changing the aperture time does not create the ascending central jet we have observed, as it never breaks the horizontal homogeneity or symmetry.

As a conclusion, none of the initial conditions tested above can create the ascending central jet that is observed in the experiments. Indeed, only the curvature was found to be able to break the horizontal homogeneity of the instability but as it is oriented towards the lighter fluid, it would rather tend to slow down the formation of an ascending central jet. Since none of the above conditions could create the desired jet, simulations with an imposed vertical velocity field at the interface in the early dynamics will be performed. Their goal is to find out what order of magnitude the vertical velocity should have to create the ascending jet, in order to orientate our search for what physical mechanism is behind the jet.

Parabolic vertical velocity

As said above, none of the initial conditions tested had the right effect, so direct numerical simulations with a vertical velocity imposed at the interface in the center of the cylinder are now performed. These simulations use the following initial incompressible velocity field:

           u x = u 0 2z σ 2 e -z 2 σ 2 1 2 -1 4 x 2 +y 2 r 2 0 x, u y = u 0 2z σ 2 e -z 2 σ 2 1 2 -1 4 x 2 +y 2 r 2 0 y, u z = u 0 e -z 2 σ 2 1 -x 2 +y 2 r 2 0 , (5.8)
with r 0 = R cyl /2 = 0.0485 m, σ = 0.001 m the initial width of the interface and u 0 being the initial speed. This forms an initial ascending parabolic vertical velocity at the interface as shown in Fig. 5.11 at the interface z = 0. This velocity field is imposed during the initialisation inside the width of the interface. The goal being to find an order of magnitude for the velocity, the parameter that is varied in the simulations is u 0 as shown in Tab. D. The two interface types (curved and planar) are also tested as it is important to assess whether the velocity imposed can counteract the effect of the curvature.

The instantaneous concentration fields for DNSInitU1 and DNSInitU2 are shown in Fig. 5.12. We see that in both cases, but especially in the DNSInitU1 simulation, the heavier fluid tends to fall on the sides of the cylinder and the lighter one tends to go up in the center. Both the Fig. 5.13c. This means that, for large u 0 the fluids are more rapidly advected but they have less time to mix. This is also visible on the vertical kinetic energy profiles of Fig. 5.13d with a much smaller energy for DNSInitU2 than for DNSInitU1 and DNSInitU3.

Aspiration Model

From these simulations, we can deduce that the physical mechanism responsible for the appearance of the central jet must be able to produce large enough initial velocities at the interface, at least of the order of the centimeter per second. Indeed, the jet in the experiments appears as soon as the two fluids make contact. It seems also safe to say that even if the curvature of the interface has an opposite effect to the desired one, it does not affect much the dynamics when the jet is created. Thus, the effect of the curvature can in the end be neglected.

Aspiration Model

Within the experimental setup, it appears that only one mechanism can create a strong enough vertical velocity field when contact between the two fluid is made: the curvature of the meniscii formed inside the holes of the grid. These meniscii are created by the pressure jump between the atmospheric pressure P a and the salt water pressure P (in the upper cylinder). This pres- sure jump is significant enough to induce a vertical velocity field of order a few cm/s. In the following, an analytical model is derived to relate the induced velocity U to the pressure jump.

It is our assumption that, depending on the sign of the pressure jump at the interface between salt water and fresh water, the fluid will either be sucked (ascending central jet, see Fig. 5.14) or injected (descending central jet) when the contact between the two cylinders is made. Focusing on the ascending jet, as it is what is most observed in the experiments, the denser fluid pressure needs to be smaller than the atmospheric one. This will create a uniform ascending velocity, U, at the contact zone which will end when P = P a .

Calculations of the jet velocity created from the pressure jump

The pressure jump at the grid is given by the Laplace law applied at the interface salt water-air:

P a -P = γ 2 r , (5.9)
with r the radius of curvature of a meniscus. Taking r = 1×10 -3 m, the order of the grid mesh size, and γ = 72.7 × 10 -3 N.m -1 , the airwater surface tension at 20°C, we find P a -P 0 ∼ 145 Pa, P 0 being the initial pressure P at the interface of the upper cylinder.

In Fig. 5.15 a sketch of an experiment in the case of an ascending vertical jet, is represented. It also shows the development of the contact area between the two fluids which is assumed to be a disk of radius:

R(t) = R cyl t t o , ( 5.10)
where R cyl = 0.098 m is the radius of the cylindrical tank and t o ∼ 0.5 s is the aperture time.

Hence, if t ≥ t o , we have R = R cyl and the aperture velocity is R cyl /t o = 0.194 m.s -1 . Then, taking into account the conservation of volume of fluids in the upper tank, we have the following relationship: (5.11) where V is the velocity of the rising free surface at the top of the upper tank (see Fig. 5.15). From the hydrostatic equilibrium, we can evaluate the pressure evolution as: (5.12) where P s is the pressure at the free surface and hh 0 is the displacement of the free surface (h 0

R 2 cyl V = R 2 U,
P -P s = ρ h G 0 (h -h 0 ) = ρ h G 0 t 0 Vdt ′ = ρ h G 0 t 0 R 2 U R 2 cyl dt ′ = ρ h G 0 t 0 t 2 U t 2 o dt ′ ,

Aspiration Model

is the initial position of the free surface at t = 0 s). Assuming a stationary flow and applying Bernoulli law to a stream line, we get: U = 2(P a -P)/ρ l .

(5.13) Thus, we can deduce that the initial jet velocity U 0 = 2 (γ/rρ l ) 1/2 = 0.54 m.s -1 with ρ l = 996 kg.m -3 , r = 1 mm and γ = 72.7 × 10 -3 N.m -1 .

Injecting Eq. (5.13) into Eq. (5.12), we find:

U 2 = U 2 0 -2 ρ h ρ l G 0 t 0 R 2 U R 2 cyl dt ′ .
(5.14)

We have:

Jet velocity    U(t) = U 0 -ρ h ρ l G 0 t 3 3t 2 o for t ≤ t o , U(t) = U 0 -ρ h ρ l G 0 t -2 3 t o for t > t o .
(5.15)

This process ends when P = P a at a time called t P . At this time, the velocity U = 0. In all of the cases considered hereafter, we have t P < t o , thus:

U(t ≤ t P ) = U 0 - ρ h ρ l G 0 t 3 3t 2 o
and U(t > t P ) = 0.

(5.16) The red dashed line represents either the fresh water-air interface or the initial fresh water-salt water interface, the blue line shows the air still traped between the grid and the fresh water being evacuated, and the green dashed line represents the free surface in the upper cylinder. Velocity [m.s

Simulations with the aspiration model

Simulations with the aspiration model

The aspiration model allows to relate a pressure jump with induced vertical velocities. In the simulations, we cannot take into account a difference in pressure, but we can implement these velocity fields and displacements to reproduce its effect. The aim of this section is to implement the aspiration model in the simulations and find out if it is the right explanation for the apparition of the ascending central jet in the experiments.

Implementation of the aspiration model

Three cases are considered for the numerical implementation of the model.

• Case I: presence of a wall at the interface that opens at the speed v PW as described in Sec. 5.1.3. A slight opening in this wall is already present initially parametrised by the initial contact radius R 0 . A two-dimensional sketch of this case is shown in Fig. 5.17a.

• Case II: in addition to the wall opening at the interface, the free surface velocity and displacement described in Sec. 5.2.2 of the aspiration model are considered. The position of the free surface is emulated by the displacement of a wall (top and bottom penalisation layers, Fig. 5.17b). This displaces the interface between the lighter fluid and heavier fluid from the top of the domain, at z = H cyl , to the bottom of it due to the periodic boundaries of the simulation. The velocity U pen = V(t)e z , where e z is the z-direction unit vector, is forced inside these layers.

• Case III: the full model is simulated in this case with the opening of the wall, the displacement of the free surface and the jet velocity U(t) forced inside a penalised layer at the interface inside the opening. This is shown in Fig. 5.17c.

It is to be noticed that Case I is extremely similar to the simulations already made in Sec. 5.1.3 with the only difference being the slight initial opening in the wall at the interface.

For further reference, the gray mask in Fig. 5.17 is named χ and defines a solid wall, whereas the red mask is called χ u and serves as a support for the forced jet velocity. The evolution of Chapter 5. The Rayleigh-Taylor instability through a grid: Central jet both these mask functions with time are shown in 2D in Fig. 5.18 for Case III. We notice the displacement of the top and bottom walls with time that would not be present in Case I and the initial slight aperture R 0 at t = 0 s in the wall at the interface, from which χ u afterwards expands. It is to be reminded that χ u only exists in Case III, thus for Case I and Case II, this red area in Fig. 5.18 would stay empty (whitish). The evolution of both velocities U(t) and V(t) with time and for U 0 = 0.54 m/s are given by Fig. 5.16b. This provides a general idea of their behaviour and thus it will not be repeated here.

For the sake of brevity, these cases were studied in 3D cylindrical configurations in Appendix A. It is found that Case II is enough to create the desired ascending jet and is thus the chosen configuration for the rest of this study as it is a good compromise between simplicity of implementation and reproducing the jet.

Moreover, in this Appendix, convergence studies were performed to determine the optimal numerical parameters for the penalisation. The conclusions are as follows: the filter parameter C s = 128 is the best compromise between sharpness of the mask and Gibbs oscillations and the velocity penalisation parameter should be taken as η u = 1 × 10 -5 s up to t o + 0.5 s and then be released to η u = 5 × 10 -4 s.

Jet velocity and aperture time

It was seen from the Case II simulation in Appendix A that taking as initial jet velocity U 0 = 45 cm.s -1 and an aperture time of t o = 0.5 s gives a too violent effect with an ascending central jet that is confined almost immediately and descending jet velocities that are around twice the values found in the experiments (Chap. 4). The aperture time of 0.5 s was the original choice as it was very close to what was observed in several of the experiments with a contact time of around 0.5 s. However, in those experiments, the lateral jets of heavier fluid descend much closer to the walls of the cylindrical tank than in the simulations. Moreover, here a linear aperture profile (Eq. (5.10)) is chosen for practical reasons which may very well not be the case in the experiments. So by testing different aperture times, the goal here is to find a good compromise (c) Case III in 2D with the wall opening at the interface at speed v PW , the top and bottom walls moving upwards at speed V(t) and a layer at the interface (hatched red) in which the upward jet velocity U(t) is forced. The gray hatched areas represents the walls defined via the penalisation method χ and χ u .

v PW v PW Heavier fluid Lighter fluid (a) v PW v PW V (t) V (t) Heavier fluid Lighter fluid (b) v PW v PW V (t) V (t) U (t)

Simulations with the aspiration model t = 0s

t = 0.1s t = 0.3s t = 0.5s Figure 5.18: Evolution of the mask functions χ (in green) and χ U (in red) with time t in Case III.

between the experimental value of 0.5 s and the position of the lateral descending jets in the simulations. Thus, we need to find the right couple of parameters (U 0 ,t o ) to get as close as possible to the experimental observations.

Several simulations are performed with 512 3 points, A = 0.003 and ν = D = 1.2 × 10 -6 m 2 .s -1 with different (U 0 , t o ) couples. These simulations are given in Tab. E with the chosen initial velocity U 0 , corresponding to a certain radius of curvature r and pressure jump P a -P, and the chosen aperture time t o . Note that for the initial velocities tested here, the pressure jump is very small (around 1 Pa) while it was much larger for U 0 = 45 cm.s -1 (around 100 Pa). Thus, we expect much less violent central ascending jet.

Series

N • U 0 [cm.s -1 ] r [cm] P a -P [Pa] t o [s] DNSU0toA 1 9
The slices along the XZ plane of the instantaneous concentration fields are given at time t ∼ 1 s in Tab. F for DNSU0toA, B and C for each initial velocity and both t o = 0.5 and t o = 0.1 s. As expected, the smaller the initial velocity, the smaller the magnitude of the ascending central jet. Moreover, the smaller the aperture time and the faster the jet develops, the larger it is and thus the closer from the walls the lateral descending jets are.

As the initial velocity increases, the descending lateral jet velocity V DLJ increases. Indeed, considering the cases t o = 0.5 s, we find that, between t = 0 and t = 1 s, V DLJ = 1.76 cm.s -1 for U 0 = 9 cm.s -1 , V DLJ = 1.3 cm.s -1 for U 0 = 4.5 cm.s -1 and V DLJ = 0.97 cm.s -1 for U 0 = 3 cm.s -1 . However, after some time, an optimal velocity seems to be reached as all of the simulations give relatively close values of V DLJ (around 2 cm.s -1 between t = 1 and t = 3 s). At t = 1 s, we are still very close to the initial conditions imposed by the model, thus it is more pertinent to measure the descending lateral jet velocities further in time. All of the V DLJ converge towards the same Chapter 5. The Rayleigh-Taylor instability through a grid: Central jet value of 2 cm.s -1 , which is quite close to the values found in the experiments. The highest U 0 used here (9 cm/s) still gives a central jet that is very easily confined, especially for t o = 0.1 s. On the contrary, the jet develops too slowly for the smallest value of U 0 . Thus, the better solution in our opinion to match the experiments is to use U 0 = 4.5 cm/s as its central jet is not so quickly confined, and the descending lateral jets are qualitatively compatible with the experimental observations. Had we a way to visualise the central ascending jet in the experiments, this decision could be verified by looking at the jet magnitude, velocity and width. As we can not visualise it, this decision can very well be challenged and will remain a hypothesis.

In addition to changing the initial velocity, changing the aperture time also has a direct impact on the development of the descending lateral jets. Indeed, as said previously, the smaller the aperture time, the faster the lateral jets develop. This is visible on Tab. F where two aperture times were tested: t o = 0.5 s and t o = 0.1 s. From this, we found that the aperture time t o = 0.1 s gives better positioned lateral descending jets. In Fig. 5.19 we verify that this position does not drastically change when using an even smaller aperture time of t o = 0.05 s. Note that the width of the lateral descending jets is slightly larger for t o = 0.05 s but their positions do not change much. The good compromise between position of the descending jets and measured experimental aperture time seems to be around t o = 0.1 s. As a conclusion, the best couple of parameters (U 0 , t o ) to be used to match the observations made in the experiments is (4.5 cm/s, 0.1 s).

Addition of a grid on the aspiration model

The ascending central jet seen in most of the experiments can be reproduced thanks to the aspiration model presented previously. Adding a grid at the interface between the two fluids enables us to almost fully copy the experimental conditions. The means of adding that grid in the numerical simulations through the penalisation method is described in Chapter 4 and will not be repeated here.

These simulations are performed in order to find which parameters for the grid porosity ϕ, grid mesh size d, initial perturbation wavelength λ peak (called λ p hereafter) and initial contact radius R 0 , give the closest result to the experimental observations. Twenty-eight high resolution simulations in a vertically extended domain were made with the aspiration model and they are summarized in Fig. 5.20. As shown in the diagrams, two porosities are tested: ϕ = 0.44 and ϕ = 0.6. For each porosity, two mesh sizes are used: d 1 = 20 mm and d 2 = 8.2 mm. For each mesh size, four initial perturbation wavelengths λ peak = λ p are studied (λ p,1 = 19 mm, λ p,2 = 13 mm, λ p,3 = 8 mm and λ p,4 = 5 mm) and almost for each of these, two initial contact radii R 0 are used: R 0,1 = 14 mm and R 0,2 = 42 mm. Some of the parameters effects on the instability dynamics have already been explored in Chapter 4 such as the grid porosity and mesh size. Here only two mesh sizes are studied (d = 20 mm and d = 8.2 mm) as they represent the two possible configurations: d > λ 0 and d < λ 0 with λ 0 the most unstable wavelength given by the linear stability analysis (Eq. (4.11)). The grid with the mesh size d = 8.2 mm is in the same configuration as the experimental grids (Tab. A), the d < λ 0 configuration.

All the simulations presented in this Section are made inside a vertically extended cylindrical geometry with the parameters given in Tab. G. The velocity penalisation parameter is η u = 1 × 10 -5 s until t o + 0.5 s and then η u = 5 × 10 -4 s. No curvature was implemented here as it was seen in Sec. 5.1.4 that it did not change much the velocity or the development of the central jet.

N • of points A ν = D [m 2 .s -1 ] η c [m 2 .s -1 ] t o [s] U 0 [cm.s -1 ] C s 1024 2 × 2048 0.003 1.2 × 10 -6 1.2 × 10 -10 0.1 4.5 128

The influence of the initial perturbation wavelength λ p

We start by investigating is the effect of the initial perturbation wavelength λ p . We saw in the previous chapter (Chap. 4) that λ p tends to have an influence on the development of the instability even with a grid that selects its own fastest growing wavelength. In order to find out if that is still true with the addition of the aspiration model, several initial perturbation wavelengths are studied for one fixed porosity, grid mesh size and initial contact radius. The parameters of the simulations considered here are summed up in Tab. H. The initial perturbation wavelengths were chosen so that two of them are in the inertial regime and two in the diffusion regime. They are given in Fig. 5.21 with the grid mesh size d, the growth rate Ω, the most unstable wavelength λ 0 and the critical wavelength λ c given by the linear stability analysis (Eq. (4.11)).

Series N

• ϕ d [mm] l [mm] λ p [mm] R 0 [

Addition of a grid on the aspiration model

The 2D slices along the XZ plane of the instantaneous concentration fields of the DNSλ p simulations at time t ∼ 3 s are given in Fig. 5.22. Note that all of the simulations are very similar. The biggest difference concerns the small scale instabilities developing near the walls that tend to develop faster for the larger λ p .

To study further the effect of the initial perturbation wavelength, the mixing zone widths computed from ⟨C⟩ XY,tot (Eq. (4.13)) with the integral (Eq. (4.30)) and profile (Eq. (4.31)) methods using ϵ = 0.01 are plotted in Fig. 5.23a. Note that, in the presence of thin jets surrounded by pure unmixed fluids, the notion of 'mixing zone' based on a horizontal homogeneity hypothesis can be questioned. The quantity L z , derived from the profile method, gives the vertical extension of the jet provided ϵ is small enough, here 0.01. The jet here becomes rapidly confined, around t = 4 s. Whereas L int yields the extent of the true mixing zone, slightly biased by the thin vertically elongated jets. Both L z and L int show that the results are independent of the initial perturbation wavelength. This is different from the case of the grid alone studied in Chapter 4 without the aspiration model.

The integral method is used to compute the growth rates α of the DNSλ p simulations given in Fig. 5.23b. A note of warning, the growth rate α is normally only measured in the self-similar regime. This regime is in all probability not reached here as the jet is rapidly confined. Still, since L int essentially captures the true mixing zone extent, one may argue that a sort of selfsimilar regime is eventually reached for the mixing region. But, as explained before, as the measure of L int is biased by the jet, no quantitative results should be drawn from the present α. Qualitatively, just as expected, the growth rates of each λ p are very close and give very similar results. Thus, the aspiration model and the creation of the central jet tend to make the instability independent of the initial perturbation.

The influence of the initial contact radius R 0

We have seen that the initial perturbation wavelength does not play a significant role when the aspiration model is taken into account, thus for the study of the other parameters effects, λ p is fixed at 13 mm. The next study is on the effect of the initial contact radius R 0 . Indeed, as it defines the size of the initial contact zone between the two fluids, we expect this parameter to have an important effect on the development of the jet. Two values are considered, R 0 = 14 mm and R 0 = 42 mm, and the simulations considered are given in Tab. I with their parameters. To illustrate the initial opened portion of the wall compared to the diameter of the cylinder, a sketch is proposed in Fig. 5.24.

The initial contact zone in the experiments is closer to the one given in Table I: Label (series and numbers), porosity ϕ, grid mesh size d, grid thread diameter l, initial perturbation wavelentgh λ p and initial contact radius R 0 for the 2 simulations considered here. mm. However, the same grids as in the experiments cannot be implemented in the simulations in part due to the lack of resolution, and also due to the much larger diffusion coefficient, which would stabilise the instability if such grids were used in the simulations (Section 4.2.3.6 of Chapter 4). Hence the grids used in the simulations have much larger threads and mesh size meaning that we do not obtain the same number of threads inside the initial opening as in the experiments. Indeed, considering an initial opening R 0 = 14 mm in the experiments, around 5 threads would be present inside it with Grid3 (d = 2 mm and l = 0.9 mm) against less than 1 in the simulations. Hence, the question is which is the right number to be preserved between the experiments and the simulations? The size of the initial opening or the number of threads in it?

On the 2D slices along the XZ plane of the concentration fields given in Fig. 5.25 at t ∼ 5 s, we notice that as the initial radius is increased, the central jet width is also increased and its development is slowed down, which is expected as the aspiration applies on a larger initial surface. The descending lateral jets are also very dependent on the initial radius of contact as their position depends on the width of the ascending central jet. Indeed, with a larger R 0 , the descending jets are much closer to the walls. The simulations that seem closer to the experimental observations are the ones made with R 0 = 42 mm. For these simulations, almost 4 threads and 3 holes are present in the initial opening contrarily to the simulations with R 0 = 14 mm where the initial opening just contains the central hole and a part of the two neighboring threads. Thus, it would seem that the most important number to preserve is the amount of threads/holes available in the initial contact zone.

The effects of the initial contact radius on the mixing zone width can be seen in Fig. 5.26. As expected, the profile method shows that the central jet develops faster R 0 = 14 mm. This jet has already reached the top of the computational domain at t ∼ 5 s in Fig. 5.25a, changing the slope of L z . Thus the instability is completely confined very rapidly for R 0 = 14 mm which is not the case for R 0 = 42 mm. The integral method indicates that the mixing zone grows more in the case with the larger R 0 . sizes have already been investigated in Chapter 4. It was found that the important dimensionless ratio is d/λ 0 : if this ratio is larger than 1 then the wavelength that develops is λ 0 whereas if it is smaller than 1, only filaments of the size of the holes of the grid can develop. Does the addition of the aspiration model change this feature?

Addition of a grid on the aspiration model

Series

In order to answer this question, 2D slices along the XZ plane of the instantaneous concentration fields are given in Fig. 5.28 at time t ∼ 5 s. The instantaneous concentration field for the grid whose ratio d/λ 0 is smaller than 1 (d = 8 mm and l = 4 mm) does show a much more filamentous sort of instability with less mixing and a visibly smaller mixing zone width. The size of the grid holes still has a very important impact on the development of the instability and the structures that can be seen. The fact that, despite being accelerated at the same speed, one jet is much less vertically extended than the other is due to the diffusion slowing down the development of the instability in the case of the smallest mesh size and the adherence to the threads.

The lateral descending jets position and velocity point towards the simulations using d = 20 mm being closest to the experiments. However, the closest grid from the experimental ones (Tab. A) is the grid with mesh size d = 8.2 mm (as d < λ 0 ). If we choose to reproduce at best the experimental grids, we cannot recover the right structures due to the much higher diffusion in the simulations. Hence, in order to be less sensitive to diffusion a coarser grid mesh with mesh size d = 20 mm is used hereafter.

The influence of the porosity ϕ

The only parameter left to study is the porosity. From the previous study, the grid mesh size is fixed at d = 20 mm and two porosities are tested ϕ = 0.44 and ϕ = 0.6 by changing the thread diameter l. The parameters of the two simulations are given in Tab. K. In order to increase the porosity of the grid, the grid mesh size d is kept the same, in order to preserve the d/λ 0 ratio, and the thread diameter is decreased, meaning that more threads are present for the same cylinder diameter. Indeed, here, six threads are available for the porosity ϕ = 0.44 while eight threads can be seen for porosity ϕ = 0.6 (Fig. 5.29).

The slices along the XZ plane of the instantaneous concentration fields at time t ∼ 5 s given in Fig. 5.29 show that as the porosity is increased, the merging and mixing of the structures developing in adjoining holes is made much closer to the interface. This is due to the lower height of the threads, meaning that the instability is less confined inside the channels delimited by the holes and can cross to the adjacent hole more easily. Thus, as the porosity is increased, the instability develops more horizontally and less vertically, yielding a smaller mixing zone extension. This also means that the lateral jets descend more slowly as they also expend more horizontally giving a smaller V DLJ : V DLJ ∼ 3 cm/s for ϕ = 0.44 and V DLJ ∼ 2 cm/s for ϕ = 0.6.

In the experiments, the threads diameter is very small meaning that mixing occurs very close to the grid, like what is obtained in the simulation with ϕ = 0.6. Moreover, it has a closer Table K: Label (series and numbers), porosity ϕ, grid mesh size d, grid thread diameter l, initial perturbation wavelentgh λ p and initial contact radius R 0 for the 2 simulations considered here. lateral descending jet velocity than the simulation with porosity ϕ = 0.44. However, both in terms of porosity and descending lateral jets velocity, the case ϕ = 0.44 is much closer to the experimental data and parameters. The parameter that gives the simulation resembling the most the experiments is ϕ = 0.6.

Conclusion

Conclusion

Rayleigh-Taylor experiments with a grid at the interface between the two fluids were made.

For a majority of these experiments an ascending central jet appeared, breaking the horizontal homogeneity of the instability. Direct numerical simulations were performed to understand the causes for the appearance of this jet, and hence to identify the mechanism that should be better controlled to erase it. The initial conditions were first investigated, and various configurations were addressed, like large initial perturbations, a curved interface, and the contact delay between the two fluids, modeled by a withdrawing wall at the interface. However, these effects could not explain the origin of the central ascending jet and the breaking of the horizontal homogeneity. Thus, we had to search for the cause elsewhere. We found that the central jet was very likely created by a jump between the atmospheric pressure and the upper fluid pressure, visible with the meniscii in the grid holes. When the two fluids make contact, the variation of pressure creates an upward velocity at the contact point, thus sucking the lighter fluid inside the heavier one at this point. It also slightly displaces the free surface above the heavier fluid. A model was derived to account for this pressure jump at the interface, which results in imposing a free surface velocity and displacement. Integrating this model into numerical simulations, and accounting for the contact time between the two fluids, we were able to recreate a central ascending jet similar to the one observed in the experiments. This model relies on several physical parameters. First, the Chapter 5. The Rayleigh-Taylor instability through a grid: Central jet contact between the two fluids, which is not instantaneous and homogeneous, is represented by (i) a prescribed time evolution R(t), (ii) an initial opening radius R 0 , and (iii) the aperture time t o . Secondly, the intensity of the jet itself, given by (iv) the upward initial jet velocity U 0 , related to the pressure jump and to the free surface velocity V. The first parameter is not determinant but the three others are. They were tuned to reproduce qualitatively the jets observed in the experiments, and also to produce descending lateral jet velocities which have a magnitude similar to the experimental ones. In addition to these physical parameters, numerical parameters necessary to implement the model through penalisation methods were studied as well: a compromise was found between precision and cost in term of numerical resources.

The final step in the journey to reproduce the experimental results was to implement a grid at the interface, again using penalisation methods. Obviously, by doing so, new parameters come into play. This brings in new challenges, namely reproducing the phenomenology and taking into account numerical limitations, since the experimental grids cannot be reproduced. It was found that, mostly due to the effect of the diffusion, matching the parameters of the experiments like the ratio of mesh size to the most unstable wavelength (given by the linear stability analysis) d/λ 0 , and the initial contact zone radius R 0 , does not necessarily bring us closer to the experimental observations. In particular, for the contact radius, the most important feature to reproduce is the number of threads inside the contact area and not the size in itself. In the end, the parameters chosen in the numerical simulations do not correspond to the parameters of the experiments (like the grid dimensions), but allow to reproduce, at least qualitatively, the overall dynamics with the central and lateral jets.

As perspectives, it could be interesting to perform zoomed simulations on the interface in order to decrease the diffusion coefficient, and therefore to investigate in more detail the instability development around threads which would have the size of the experimental grids. Moreover, within the framework of the aspiration model in the simulations, it was found that even very small pressure differences (< 1 Pa) gave an ascending central jet, thus the pressure needs to be finely tuned to be able to control the jet appearance. The value from which the pressure jump becomes too small to create the jet could be found through additional measurements: with this, a new experimental setup could be proposed to fully erase the jets.

Conclusion Conclusion

The Inertial Confinement Fusion is a very complicated process involving several different physical phenomena. Our interest lies in the interfacial instabilities happening between the ablator (heavy) and the fuel (light). These instabilities are subjected to variable accelerations, thus we studied two elementary configurations in this thesis: a stable configuration under periodic acceleration and an unstable configuration under a constant acceleration.

The first configuration studied is the stable one under periodic acceleration also known as the Faraday instability. It consists of the superposition of two fluids, chosen as miscible for our study, of different densities with the lighter one above. Under a periodic vertical acceleration, the interface between the two fluids can be destabilized. It is a parametric instability involving the resonance between the forcing frequency and the natural frequency of the interface. Both the wavebreaking mechanism leading to turbulence of the Faraday waves and the mode selection mechanism of the primary wave were investigated.

Considering the quantization induced by the geometry of the tank and the large forcing accelerations used in our experiments, several subharmonic modes are given as unstable by the Floquet linear theory. However, only one of them is selected by the primary Faraday wave. It was demonstrated that the selection of the mode results from a nonlinear competition favoring the modes with smaller wavelength. Indeed, when the primary wave amplitude grows, the natural frequencies of the system decrease. This means that, when the amplitude reaches a critical value, the supercritical modes (negative detuning) fall out of resonance with the forcing frequency and are thus damped while the subcritical mode (positive detuning) becomes more resonant. The subcritical modes are thus favored by the selection mechanism.

In order to study the transition to turbulence of the Faraday miscible standing waves, the breaking of these waves must be understood. To this purpose, two theories were constructed. The first is referred to as 'global' and is based on the horizontal averaging process. It evidenced that the secondary instability breaking the wave is due to the oscillations of the primary wave. A criterion giving the critical steepness of the wave at which wavebreaking is expected to occur was derived. However, this breaking only happens at the node of the waves, which cannot be explained by this theory. Thus, a second theory, referred to as 'local' was thought of. It focused on the flow in a local frame attached to the node of the wave which showed the importance of the shear in the development of the secondary instability. Indeed, it was shown that the unstable modes could be either of a parametric resonance type or of a Kelvin-Helmholtz type. To determine to which of these types the mode belongs to, we compared the theories to the experimental and numerical data. The Thorpe displacement method, indicating an overturning, was used to measure the primary wave amplitudes and the wavenumber associated with the wavebreaking. From these comparisons, we found that the secondary instability is of parametric resonant nature, and develops earlier than the Kelvin-Helmholtz type during the growth of the primary wave. Indeed, the results show that the wavebreaking is detected for wave steepnesses ka ∼ 0.75 inside the parametric resonant band. Thus, the wavebreaking comes from a secondary parametric resonant subharmonic instability developing when the primary Faraday wave amplitude reaches a critical value: it is of subcritical nature.

Numerical simulations with 2D initial conditions perturbed along the spanwise direction were performed to show that it is indeed the wavebreaking that leads to turbulence. From these, it was evidenced that the secondary instability developing at the node of the wave is the origin of the final transition to turbulence.

The second configuration studied during this thesis is the unstable one under constant acceleration, called the Rayleigh-Taylor instability. It consists in superposing two fluids of different densities with the heavier one above in a downward gravity field. This configuration is unstable to any perturbation in an inviscid framework without surface tension. Thus, it can be easily understood that well controlled Rayleigh-Taylor experiments are very difficult to perform. This is a true challenge in the understanding of the instability as a number of questions are still raised and, in particular, the dependence on initial conditions.

In the hope of better controlling the initial conditions, we proposed a new type of Rayleigh-Taylor experiment. In this new setup, a grid is positioned at the interface between the two fluids, whose goal is to stabilise the heavier fluid during the filling of the tank; then, when the contact is made, to let the instability develop freely once the contact between both fluids is made. However, since the grid is present at the interface throughout the experiment, it is necessary to ensure that it does not impact the dynamics of the instability.

Both experimental and numerical works were thus done to understand the effect the grid may have on the development of the instability. From these, it was found that the grid can have a significant impact on the shapes and sizes of the developing structures. Indeed, the grid does not let perturbations with wavelength larger than the holes size develop. This means that choosing a grid with a mesh size larger or smaller than the most unstable wavelength given by the linear stability analysis greatly impacts the resulting instability. In the experiments, the grids used did not let the most unstable wavelength develop. This means that thin filamentous structures of the size of the holes of the grid developed. As the grid thread diameters were very small, the mixing of the filaments with their neighbour happened very close to the grid, which allowed to get larger structures resembling the typical Rayleigh-Taylor ones. Thus, the grid does ensure a certain control over the initial conditions as we know that only the small wavelengths are allowed to develop.

It was evidenced that the grid's impact on the dynamics decreases with time. Indeed, with the addition of the grid comes the creation of a jump on the concentration profiles. However, this jump tends to disappear as the instability develops and more mixing occurs. As in the 'classical' Rayleigh-Taylor instability, a self-similar regime seems to be reached with values of growth rates Chapter 5. Conclusion

for the experiments very close to what is obtained in the literature (α = 0.034). The values for α in the simulations are smaller with grid than without, but this value is also impacted by the porosity of the grid, the initial perturbation wavelength imposed, and the way it is computed. Indeed, with the apparition of thin filaments, the growth rate and mixing zone width become very sensitive to the method of computation. The results should thus be looked at with caution.

During the vast majority of the experiments, an ascending central jet developed as soon as the fluids made contact. As we could only see the development of the instability in the bottom cylinder, it was the lateral descending jets that were observed. If we want this new type of experiment to be a good Rayleigh-Taylor experiment candidate, we need to understand what causes the apparition of this jet and suppress it.

Thus, numerical simulations were performed to reproduce at best the conditions of the experiments and determine the origin of the jet. It was found that the pressure jump at the heavy fluid-air interface generates meniscii inside the holes of the grid, so that when the contact is made, an upward vertical velocity is created at the contact point and is accompanied by the displacement of the free surface at the top of the heavy fluid. A model was derived and implemented in the simulations. With it, we were able to recreate the jet and find parameters that would get us close to the experimental observations. To suppress this jet, the pressure inside each fluid would need to be equalized.

For future research, it would be interesting for both the understanding of the grid effect and of the jet, to be able to simulate a finer grid. For this the resolution could be improved and the diffusion coefficient could be decreased. Decreasing the diffusion coefficient would indeed allow much smaller wavelengths to develop, which would mean recovering an earlier mixing with and without jet. Moreover, the exact pressure jump value for which the central jet stops developing could be found.

A.1. The three cases of the aspiration model

A.1 The three cases of the aspiration model

This section presents the three cases described when implementing the aspiration model in Chapter 5. Case I has a wall opening at the interface between the two fluids in order to account for the delay in contact of the two fluids as given by Eq. (5.7). Case II implements on top of that the free surface velocity and displacement (Sec. 5.2.2) and Case III adds the jet velocity U at the interface (Eq. (5.15)). Three simulations are performed, one for each case and a qualitative comparison of the density fields is made by showing the vertical slices along the XZ plane of the instantaneous concentration fields.

N • of points A ν [m 2 .s -1 ] D [m 2 .s -1 ] η u [s]
η c [m 2 .s -1 ] 512 3 0.003 1.2 × 10 -6 1.2 × 10 -6 1 × 10 -4 1.2 × 10 -10

Table A: Parameters of the simulations with the number of points, Atwood number A, viscosity ν, diffusion coefficient D and velocity and scalar penalisation parameters η u and η c respectively.

The scalar remains unpenalised inside χ u (the mask imposing the jet velocity at the interface) while the scalar penalisation parameter inside χ (the mask defining the top/bottom and side walls) is η c and the velocity penalisation parameters are taken in both masks as η u . All parameters are given in Tab. A.

A.1.1 Case I: Wall opening

Case I simulates the Rayleigh-Taylor instability in a cylindrical geometry with a wall at the interface between the two fluids that opens with time, at the speed v PW = 0.194 m/s and aperture time t o = 0.5 s. The initial contact radius is R 0 = 3.5 mm.

As expected from the study made in Sec. 5.1.3, no central ascending jet appears on the vertical slices in Fig. A.1. A slight difference can be seen in the development of the instability when comparing the center and the sides of the cylinder. This may be due to the no-slip boundary conditions of the cylinder walls. This effect is also combined with the fact that the maximum growth rate given by Eq. (4.11), Ω RT , is slightly bigger than the rate of aperture of the wall, Ω o (Ω RT = 2.49 s -1 > Ω o = 2 s -1 ). Thus, the instability starts developing before the wall is completely opened. In any case, this difference is minimal and does not impact much the development of the mixing layer.

Figure 1 :

 1 Figure 1: Schematics of indirect (upper left) and direct drive (upper right) approaches in Inertial Confinement Fusion experiments. The four bottom images describe the implosion process. First, a spherical capsule is prepared at t = 0 with a layer of DT fuel on its inside surface. Second, the capsule surface absorbs energy which ablates the outer materials. Third, the capsule imploses. Lastly, a hotspot of DT is formed surrounded by colder and denser DT fuel. This last image on the bottom right evidences the effects of the hydrodynamics instability. Image taken from Betti & Hurricane [2016].

  Figure1.1: The mechanism of parametric instability illustrated by a sketch of a pendulum at successive phases with the corresponding modulation δg of the apparent gravity (the extremal positions of the pendulum are indicated by dotted lines and its velocity by the blue arrow). The apparent gravity is stronger than average when the pendulum moves downward and is weaker than average when the pendulum moves upward. This is possible when the excitation period is half T, the natural period of the pendulum. Figure greatly inspired by the figure 2 of Benielli & Sommeria [1998].

Figure 1 . 2 :

 12 Figure 1.2: Sketch of the Faraday instability with the heavier fluid in dark blue and the lighter one in lighter blue. The interface deformation ξ due to the vertical oscillations G(t) is shown in dotted line. The direction of the gravity G 0 is indicated with an arrow.

  Figure 1.3: (a) Single-frequency forcing at high viscosity. A pattern of essentially parallel lines is observed. The forcing frequency is 2ω/2π = 80 Hz. (b) Hexagons produced by two-frequency forcing. The forcing acceleration is f (t) = a[cos(χ) cos(4ωt) + sin(χ) cos(5ωt + ϕ)] where ω/2π = 14.6 Hz, ϕ = 75 • and χ = 45 • . The pattern is observed at the primary transition from the flat surface. The mixture used in both cases is composed of 88% glycerol and 12% water with a kinematic viscosity of 1.00 cm 2 .s -1 . Pictures taken from Edwards & Fauve [1994]. (c) Top view of the oscillating Faraday pattern. Picture taken from Domino et al. [2016].

Figure 1 . 4 :

 14 Figure 1.4: (a) Stability boundaries for ideal fluids (ν h = ν l = 0 m 2 .s -1 ). The tongues correspond alternately to subharmonic (SH) and harmonic (H) responses. Fluid parameters are ρ h = 519.933 kg.m -3 , ρ l = 415.667 kg.m -3 , γ = 2.181 × 10 -6 N.m -2 and 2π/ω = 100 Hz. (b) Stability boundaries for the Full Hydrodynamic System (FHS) for viscous fluids of viscosity ν h = ν l = 7.516 × 10 -8 m 2 .s -1 and other parameters are as in (a).Here, a c is the critical excitation amplitude, g is the gravitational acceleration and k is the wavenumber. This figure was taken from[START_REF] Kumar | Parametric instability of the interface between two fluids[END_REF].

  Figure 1.6: Pictures of a Faraday experiment with miscible fluids with F = 0.30, ω = 4.52 rad.s -1 and A = 0.03 (Briard et al. [2020]) at different times: (a) initial configuration at t = 0 s, (b) deformation of the interface between the fluids at t = 46 s and (c) growth of the deformation amplitude at t = 58 s, (d) breaking of the wave at t = 61 s (see Chap. 2), (e) turbulent mixing at t = 111 s and (f) final saturated state at t = 378 s (Gréa & Ebo Adou [2018] and Briard et al. [2020]).

  Figure 1.7: Comparison between the obtained L sat and the prediction L sat , for the three Atwood numbers: A = 0.015 in blue, A = 0.03 in red and A = 0.045 in green. Symbols: □: L sat with probe; △: L sat with image and ×: L sat with DNS. (a) The A-dependence; (b) The ω-dependence. This figure was taken from Briard et al. [2020].

Figure 2 . 1 :Chapter 2 .Figure 2 . 2 :

 21222 Figure 2.1: The breaking of a Faraday wave in the FARAMIX experiment. (a): Time series images from the camera zooming on one wavelength and presenting two oscillation periods of the primary wave. This illustrates the different stages of the wavebreaking with first a blurring of the interface at the node followed by a roll-up. This case corresponds to the EXPb5 experiment which parameters are detailed in Tab. A. (b): Visualisation of the interface at wavebreaking in the direct numerical simulation DNSd3 which parameters are given in Tab. B. The reference frame as well as the acceleration direction are also indicated.

Figure 2 . 3 :

 23 Figure 2.3: Picture of the tank with its dimensions (W = 94.6 cm, H = 67 cm and D = 11 cm). The tank is filled with fresh water coloured with blue dye and salt water below it. Snapshot during the instability.

Figure 2 . 4 :

 24 Figure 2.4: Map in the (ω, F) plane of the configurations explored experimentally in Briard et al. [2020]. Symbols are for various Atwood numbers: A = 0.015 blue squares, A = 0.03 red triangles, and A = 0.045 green circles. One symbol can represent multiple experiments. Open symbols mean that the instability was not triggered. The dashed line indicates the limit of the admissible parameters. This image was taken from Briard et al. [2020].

Figure 2 . 5 :

 25 Figure 2.5: Drawing of a rectangular basin of width W, height H and depth D, filled with two fluids of different densities, the heavier below, without free surface. The amplitude a of the wave is indicated.

Chapter 2 .

 2 .27) The subcritical transition to turbulence of Faraday waves in miscible fluids Contributions of the horizontal walls at z = -H/2 and z = H/2: Considering the bottom wall z = -H/2, the fluid velocities in the viscous layer are:

Figure 2 . 6 :

 26 Figure 2.6: Plot of the three damping coefficients as a function of the wavenumber k for a diffuse interface of width δ = 0.5 cm and an Atwood number A = 0.03. The vertical dashed black line gives a typical wavenumber used in the experiments.

Figure 2 . 7 :

 27 Figure 2.7: Stability diagram for(2.9) in a non dimensional frequency ω/ AG 0 k 2,0 and forcing F plane. The colored regions correspond to the first subharmonic instability band associated with the different modes of the tank (the mode number is indicated in the figure). The diagram is obtained using the damping coefficient γ = γ δ + γ w at three different Atwood numbers and considering an interface thickness δ = 1 cm. The neutral curves (thick black plain lines) have a slight dependence on the Atwood number showing that they are not completely superimposed. The symbols correspond to the parameters taken in the experiment in Tab. A. The shapes indicate the Atwood number and the colors reveal which mode is eventually selected.

2. 3 .

 3 Mode selection mechanism of the Faraday wave

Chapter 2 .FFigure 2 . 8 :

 228 Figure 2.8: Diagram of the process of mode suppression, showing the evolution of the resonance of subcritical (∆ ≥ 0) (a) and supercritical (∆ < 0) (b) modes.

Figure 2 . 9 :

 29 Figure2.9: Parameters of the DNS (symbols) in the stability diagram ω/ AG 0 k 2,0 -F. The colored areas correspond to the subharmonic instability tongues for the modes (2,0) and (3,0) accounting for viscosity and a diffused interface δ = 1.8 cm. The symbols' color indicates which mode emerges from the simulation, the mixed colors express that both modes can be observed. The two series DNSb and DNSc (see Tab. B) are presented here starting from point A. In the DNSb group, the initial amplitude ratio r between modes (2,0) and (3,0) is set at r = 0.5 and the forcing frequency is decreased. In the DNSc group, the frequency is fixed and r is varied. The point corresponding to DNSe is also placed.

Chapter 2 .

 2 The subcritical transition to turbulence of Faraday waves in miscible fluids r Cases with ω varying (a) Cases with r varying

Figure 2 . 10 :

 210 Figure 2.10: Mode selection in 6 DNS of Fig. 2.9. (a) Cases corresponding to DNSc (see Tab. B) with r varying and A = 0.030, ω = 3.07, F = 0.8. The amplitude of the interface deformation is ϵ = 3 cm. (b) Cases corresponding to DNSb (see Tab. B) with ω varying and A = 0.030, F = 0.8, r = 0.5, ϵ = 1.5 cm.We put the slices of width W of the concentration field at four instants starting from the initial interface at t = 0 and ending when wavebreaking occurs; pure fluids remain in white while mixed fluid with C = 0.5 ± 0.15 is in black.

Figure 2 . 11 :

 211 Figure 2.11: Slices of width W of the concentration field at wavebreaking for the series DNSb with ω = 2.4, 2.6 and 2.8 shown as inserts of the Fourier modes of the interface.

  Figure 2.12: Frameworks applied to model the wavebreaking of the primary wave (a) and detailed in Sec. 2.4. For the global approach (b), we consider the stability of the horizontally averaged concentration profiles. For the local approach (c), we study the development of small perturbations of wavenumber k wb at the node of the primary wave.

Chapter 2 .

 2 The subcritical transition to turbulence of Faraday waves in miscible fluids

Figure 2 .

 2 Figure 2.13: Stability analysis of Eq. (2.40) with L(t) = L 0 (1 + β cos(ωt)) and represented in the plane(Ω 2 0B /ω 2 , F) with Ω 0B = (2AG 0 /L 0 ) 1/2 .The two colored areas indicate the first subharmonic tongues with β = 0 and β = 0.7 respectively. The dashed blue lines correspond to the approximation (2.48) while the red dashed dotted ones correspond to(2.49).

2 .

 2 .53) i = 1 or 2 whether we are looking at the fluid above or below the interface, and the continuity of Chapter The subcritical transition to turbulence of Faraday waves in miscible fluids pressure at the interface gives:

Figure 2 .

 2 Figure 2.14: Stability curve of Eq. (2.67) for F = 0.7 and ∆ = 0 in the κ -ka plane or the P -Q plane (insert) corresponding to the classical representation of the Mathieu equation. The blue colored areas show the Kelvin-Helmholtz (KH) and parametric resonance (PR) instability regions. The dashed curve corresponds to P = 0. The area corresponding to P < 0 is located above the dashed curve in the κ -ka representation. The critical wave steepness value indicated by the black dotted curve corresponds to the criterion (2.69).

Chapter 2 .Figure 2 .

 22 Figure 2.15: (a) Visualisation of the interface in experiment EXPa7 (see Tab. A) at the amplitude maximum just before the wavebreaking and compared to a sinusoidal profile (red line). (b) Mean concentration profile from the experiment EXPa7 (dashed blue line) and a sinusoidal interface (black continuous line). (c) L int plotted as a function of |ξ p | at a maximum amplitude just before the wavebreaking for all the experiments of Tab. A. The values for |ξ p | are evaluated from the crest-to-crest distance of the wave measured directly on the images (lighter blue squares) or using the mean concentration profiles (darker blue circles). The two arrows correspond to the EXPa7 case shown in (a) and (b).

  2.16). However, Chapter 2. The subcritical transition to turbulence of Faraday waves in miscible fluids

Figure 2 .

 2 Figure 2.16: Procedure for the wavebreaking detection. Left: Visual criterion from the calibrated camera image. Right: The Thorpe displacements δ T evaluated by sorting the concentration field in each vertical transect of the calibrated image. The colorbar indicates the displacement in [mm].

Figure 2 .

 2 Figure 2.17: Time evolution of the mixing zone width and the mean concentration profiles at different instants extracted from the DNSd2 (see parameters in Tab. B). (a) Evolution of the integral lengths L int and L int,s , computed from the mean and sorted concentration profiles respectively, as a function of ωt. The dashed line corresponds to the integral scale L int,d expressing the thickening of the interface by diffusion only. The star symbol at ωt = 17.75 indicates the wavebreaking detected by the Thorpe displacement. The horizontal lines correspond to the theoretical wavebreaking predictions, here converted in terms of integral length. (b) Mean concentration profiles at different times corresponding to the local maxima of L int in (a) and renormalized by the integral mixing zone width L int . The inserted images illustrate the state of the interface at the same instants.

Figure 2 .

 2 Figure 2.18: Parametric instability bands (colored regions) and experimental data (symbols) corresponding to the wavebreaking in a κ -ka representation. The parameters for the experiments are detailed in Tab. A. The left (a,c) diagrams correspond to the local theory (PR1 instability band). The right (b,d) ones are associated with the global theory. We place the theoretical instability zones with ∆ = 0 in the (a,b) figures and ∆ = 0.04 for the (c,d) figures, all of them with F = 0.3 and F = 0.7. The symbol colors indicate successively (a) the forcing parameter F, (b) the Atwood number A, (c) the detuning ∆ and (d) the primary wave mode in the experiments.

Chapter 2 .Figure 2 .

 22 Figure 2.19: Same as Fig. 2.18 but for the DNS data corresponding to the series DNSa (see Tab. B).

  Figure 2.20: Wave steepness at wavebreaking as a function of the forcing parameter F for the DNSd cases (symbols) with ω = 4.29 rad.s -1 and A = 0.045 (see Tab. B). The dashed and dotted curves correspond respectively to the local and global criteria. Inserts: Visualisation of the concentration field in DNSd3 and 11 just after the wavebreaking.

Chapter 2 . 3 Figure 2 .

 232 Figure 2.21: Visualisation of the interface at different times of DNSd3 (left) and DNSf(right) corresponding to the parameters in Tab. B. Wavebreaking is detected at ωt = 14.5 for both simulations.

Figure 3

 3 Figure 3.1: (a) Rayleigh-Taylor bubble and spike. The penetration distance, from the initial position of the density interface, of the light bubble (ρ 2 ) is written h 1 and the corresponding penetration of the heavy spike (ρ 1 ) is written h 2 . This figure was taken from Andrews & Dalziel [2010]. (b) Visualisation of a vertical slice of the density field in the Rayleigh-Taylor instability mixing layer at late time for an Atwood number of 0.5. The gravity is directed downwards. The heavy fluid is in black while the light fluid is in yellow. This image was taken from Cabot & Zhou [2013].

Figure 3 . 2 :

 32 Figure 3.2: Romanesco broccoli as a self-similar object. Image taken from the Wikipedia page.

Chapter 3 .

 3 Figure 3.3: (a) Comparison of measurement techniques for the Rayleigh-Taylor instability growth rate, α, including the Reynolds number (Re) dependence. The mixing layer thickness (h) and its growth rate ( ḣ) are normalized by the Atwood number (A), gravity (g) and time (t). Figure taken from Cabot & Cook [2006]. (b) Comparison of α b from numerical simulations (points) and experiments in histogram form.LEM is from Linear Electric Motor[START_REF] Dimonte | Density ratio dependence of Rayleigh-Taylor mixing for sustained and impulsive acceleration histories[END_REF]), RR from 'rocket rig'[START_REF] Read | Experimental investigation of turbulent mixing by Rayleigh-Taylor instability[END_REF]), K from Kucherenko et al.[START_REF] Kucherenko | Experimental study of the gravitational turbulent mixing self-similar mode[END_REF]), AS is from Andrews and Spalding[START_REF] Andrews | A simple experiment to investigate two-dimensional mixing by Rayleigh-Taylor instability[END_REF]). Figure taken from Dimonte et al. [2004].

  Figure 3.3: (a) Comparison of measurement techniques for the Rayleigh-Taylor instability growth rate, α, including the Reynolds number (Re) dependence. The mixing layer thickness (h) and its growth rate ( ḣ) are normalized by the Atwood number (A), gravity (g) and time (t). Figure taken from Cabot & Cook [2006]. (b) Comparison of α b from numerical simulations (points) and experiments in histogram form.LEM is from Linear Electric Motor[START_REF] Dimonte | Density ratio dependence of Rayleigh-Taylor mixing for sustained and impulsive acceleration histories[END_REF]), RR from 'rocket rig'[START_REF] Read | Experimental investigation of turbulent mixing by Rayleigh-Taylor instability[END_REF]), K from Kucherenko et al.[START_REF] Kucherenko | Experimental study of the gravitational turbulent mixing self-similar mode[END_REF]), AS is from Andrews and Spalding[START_REF] Andrews | A simple experiment to investigate two-dimensional mixing by Rayleigh-Taylor instability[END_REF]). Figure taken from Dimonte et al. [2004].

3. 3 .

 3 Figure 3.4: (a) Experimental apparatus of the Rocket Rig. Figure taken from Read [1984]. (b) Left: Experimental configuration for the Linear Electric Motor. Right: Closeup view of projectile in acceleration (top) and brake regions (bottom).Figure taken from Dimonte & Schneider [2000].

Chapter 3 .Figure 3

 33 Figure 3.5: (a) Experimental setup of the overturning tank. Figure taken from Andrews & Dalziel [2010]. (b) Tilted-rig experiments with an initial tilt angle of 0.9 • and an Atwood number of A = 0.048 taken at different times: (1) t = 2 s, (2) t = 2.2 s, (3) t = 2.4 s and (4) t = 2.6 s. Figure taken from Andrews & Dalziel [2010].

Figure 3

 3 Figure 3.6: (a) Low aspect-ratio Rayleigh-Taylor box, showing its barrier fully withdrawn. Photo taken from Lawrie [2009] (b) Composite sliding barrier tank with the barrier being withdrawn. The light grey shading represents the nylon fabric wrapped around the dark grey metal plates. (c) Illustration of process by which the fabric (dashed) is removed through the gap between the plates. Images taken from Andrews & Dalziel [2010].

  [2003] to investigate the effect of an initial tilt on the interface,[START_REF] Jacobs | Rayleigh-Taylor instability in complex stratifications[END_REF], Lawrie[2009] and Davies Wykes &Dalziel [2014] to study the development of the Rayleigh-Taylor instability in complex stratifications. Chapter 3. State of the art : The Rayleigh-Taylor Instability (a) (b) Figure 3.7: LIF images of a typical experiment at (a) t = 5 s and (b) t = 10 s. Pictures taken from Dalziel et al. [1999].

Figure 3 .

 3 Figure 3.8: Schematic of the water tunnel experimental setup. Figure taken from Andrews & Dalziel [2010].

Figure 3 . 9 :

 39 Figure 3.9: Mixing layer development in the water tunnel experiment. Upper, clear, heavy water mixes by Rayleigh-Taylor instability with lower, dark, light water generating turbulent mixing. Picture taken from Andrews & Dalziel [2010].

Figure 3 . 10 :

 310 Figure 3.10: Bubble penetration against τ 2 = AG 0 t 2 /H (H × H × H being the computational domain size) for ρ 1 /ρ 2 = 20. The different symbols represent different meshes and even 2D simulations. Figure taken from Youngs [1991].

Figure 3 .

 3 Figure 3.11: 2D density slices showing dominant bubbles at Agt 2 /L ∼ 14 for the simulation code ALEGRA and ∼ 22 for all other simulation codes. Figures from Dimonte et al. [2004].database to investigate several effects such as the evolution of the mixing parameter or the bubble densification shown in Fig.3.11. They found smaller values of the growth rate α in the simulations than in the experiments which they attributed, at least partially, to larger scales initial conditions in the experiments.Very recent work studied the suppression of the Rayleigh-Taylor turbulence by means of a time-periodic acceleration through 3D simulations[START_REF] Boffetta | Suppression of Rayleigh-Taylor turbulence by timeperiodic acceleration[END_REF]). They discovered that an alternating acceleration could relaminarise the turbulence. This periodic vertical acceleration,

Chapter 3 .Figure 3 .

 33 Figure 3.13: Vertical sections of the density field for: (a) a three-dimensional Navier-Stokes simulation (bulk flow), (b) a two-dimensional Navier-Stokes simulation, both at time t = 25, and for (c) a threedimensional Darcy simulation (porous medium) and (d) a two-dimensional Darcy simulation, both at time t = 30. White represents light fluid while black represents heavy fluid and gravity is in the vertical direction. Figure taken from Boffetta & Musacchio [2022].

Figure 4 . 1 :

 41 Figure 4.1: Operating principle of the new type of Rayleigh-Taylor experiment.

Figure 4 . 2 :

 42 Figure 4.2: Sketch of the water-air interface.

Figure 4 . 3 :

 43 Figure 4.3: Inviscid growth rate Ω (black line), as a function of the horizontal wavenumber k, (a) with surface tension γ T (see Eq. 4.8 and (b) without surface tension (Ω = √ G 0 kA) in the case of two uniform inviscid fluids (water and air) of constant density separated by a horizontal boundary. In (a), the critical wavenumber k c (in dashed red line), the wavenumber of the most unstable mode k 0 (in dashed blue line) and the maximum growth rate Ω 0 (in green dashed line) are also plotted.

Figure 4 . 4 :

 44 Figure 4.4: Experimental setup. (a) The two cylindrical tanks are in contact, with H cyl and R cyl the height and radius of a cylindrical tank respectively; the grid is visible in between. (b) LIF method with the bottom cylindrical tank and the laser sheet visible.

Figure 4 . 5 :Figure 4 . 6 :

 4546 Figure 4.5: Process of the upper tank filling. The upper cylinder with the grid fixed at the bottom is visible along with the tube fixed at the top of the cylinder sucking the air, thus sucking the salt water inside the cylinder.

Figure 4 .

 4 Figure 4.7: Instantaneous concentration field of a Rayleigh-Taylor instability beneath the grid in the plane of the laser sheet at t = 2.8 s after contact for experiment EXPb5 with A = 0.0037 and a mesh size d = 1.8 mm (see Tab. B). In dark is the heavier fluid falling into the lighter fluid in gray. The grid is visible at the top of the image.

Chapter 4 .

 4 The Rayleigh-Taylor instability through a grid: Influence of the grid Dispersion relation with viscosity and diffusivity

Figure 4 .

 4 8 also shows the smallest wavenumber allowed by the experimental grid having the

Figure 4 . 8 :

 48 Figure 4.8: Viscous growth rate Ω (Eq. (4.12)) as a function of the wavenumber k for a sharp interface. The critical wavenumber k c is indicated by the red dashed line, the fastest growing wavenumber k 0 is shown in dashed blue line and the smallest wavenumber allowed by the experimental grids k g = 2π/λ g (with λ g = 2 mm) can be seen in dashed green. The parameters are A = 0.0037, G 0 = 9.81 m.s -2 , ν = 1 × 10 -6 m 2 .s -1 and D = ν/700.

Figure 4 . 9 :

 49 Figure 4.9: (a) Growth rate Ω (Eq. (4.12)) as a function of the wavenumber k. The critical wavenumber k c is given in red dashed line, the fastest growing wavenumbers k 0 is shown in dashed blue line and the smallest wavenumber allowed by the experimental grids k g (with λ g = 2 mm) can be seen in dashed green line. The parameters are A = 0.0009 and A = 0.0134, G 0 = 9.81 m.s -2 , ν = 1 × 10 -6 m 2 .s -1 and D = ν/700. (b) Most unstable wavelength k 0 as a function of the Atwood number A. The minimum wavenumber k g for Grid1 and Grid3 are given in dashed cyan and green lines respectively.

Figure 4 .

 4 Figure 4.10: Evolution of the viscous growth rate Ω(k) with the increase of the diffusion thickness of the interface δ (dark area) from 0 to 5 mm. The evolution of the most unstable wavenumber k 0 (blue area), the critical wavenumber k c (red area) are also given and the grid wavenumber k g (dashed green line) are also given. The parameters are A = 0.0037, G 0 = 9.81 m.s -2 , ν = 1 × 10 -6 m 2 .s -1 , D = ν/700.

  A and densities ρ h and ρ l of the heavy and light fluids respectively. Whether the ascending central jet (ACJ) was seen in the experiment or not, along with the velocity of the descending lateral jets (V DLJ ) are given in the last two columns.

Figure 4 . 11 :

 411 Figure 4.11: Instantaneous concentration field of a Rayleigh-Taylor instability beneath a grid in the plane of the laser sheet at times t ∼ 2.3, 3.8 and 5.7 s (from top to bottom) after contact for experiment EXPc2 with A = 0.0036 and a mesh size d = 2.5 mm (see Tab. B). The heavier fluid is in dark and the lighter fluid is in gray. The grid is visible at the top of the image. The dark red arrows show the falling movement of the heavier fluid and the dark blue arrow shows the ascending movement of the lighter fluid.

Figure 4 .

 4 Figure 4.12: Mean cylindrical concentration ⟨C⟩ XY as a function of z for EXPb2, EXPb5, EXPb9, EXPb13 and EXPb14 of Tab. B at times (a) t ∼ 2 s and (b) t ∼ 4 s.

Figure 4 .

 4 Figure 4.13: (a) Mean cylindrical concentration profiles for experiment EXPb5 at times t ∼ 2 s and t ∼ 4 s with the jump in concentration shown at each time by the red lines and, in dashed black line, a typical Rayleigh-Taylor mean concentration profile without grid. (b) Instantaneous concentration field of EXPb5 at times t ∼ 2 s and t ∼ 4 s.

Figure 4 .

 4 Figure 4.14: (a) Mixing zone width L as a function of time t for EXPb2, EXPb5, EXPb9, EXPb13 and EXPb14 of Tab. B. (b) Growth rate α as a function of time for the same experiments. The value of α given in the legend is the mean value computed between the two dotted gray line inside the gray area.

Figure 4 .

 4 Figure 4.15: 2D perturbation spectrum E (Eq. 4.19) as a function of k (continuous blue line) with k peak = 20 m.s -1 (red dashed line) and the infrared slope taken as s i = 1 (orange dotted line).

Figure 4 .

 4 Figure 4.16: Initial condition on the concentration profile C(z, t = 0) as a function of z for σ = 0.001 m (see Eq. 4.18 and Eq. 4.20).

Figure 4 .

 4 Figure 4.17: Sketches of the numerical simulation geometry. (a) Thin Cylinder: height 2H cyl = 60 cm and diameter 2R cyl = 19.4 cm. (b) Large Cylinder: height 2H cyl = 19.4 cm and diameter 2R cyl = 19.4 cm. The darker blue is the heavier fluid while the lighter blue is the lighter fluid. The gray areas correspond to the penalized portions of the cube representing solid walls and the red hatched area corresponds to the penalised portions of the cylinder representing the grid.

Chapter 4 .

 4 The Rayleigh-Taylor instability through a grid: Influence of the grid

Figure 4 .

 4 Figure 4.18 shows the instantaneous concentration fields at time t ∼ 5 s of simulations DNSCubVSCyl1 and DNSCubVSCyl2 (see Tab. D). Aside from the geometrical differences between configurations, both images show a developed Rayleigh-Taylor instability. The penalisation layers do not seem to have a qualitative effect on the development of the instability.

Figure 4 .

 4 Figure 4.18: Instantaneous concentration fields of simulations in (a) a cube (DNSCubVSCyl1), (b) a cylinder (DNSCubVSCyl2), at t ∼ 5 s with A = 0.003 (see Tab. D). The heavier fluid is in red, the lighter one in blue and the pure fluids are transparent.

Figure 4 .

 4 Figure 4.19: (a) Horizontal mean concentration profiles ⟨C⟩ XY at time t ∼ 5 s and (b) evolution of the mixing zone size L with time t for the two simulations DNSCubVSCyl (see Tab. D) presented in Fig. 4.18.

Chapter 4 .

 4 The Rayleigh-Taylor instability through a grid: Influence of the grid Label D cyl [m] 2H cyl [m] k peak [m -1

Figure 4 .

 4 Figure 4.20: (a) and (c) Instantaneous 3D concentration field and (b) and (d) vertical slice along the XZ plane for the configuration with a cylindrical tank of diameter 19.4 cm and height 60 cm (DNSTC (a)(b)) and for the configuration with a cylindrical tank of diameter 19.4 cm and height 19.4 cm (DNSLC (c)(d)). The heavier fluid is in red and the lighter one in blue. The slices (b) are cropped to match the size of the slices (d).

Figure 4 .

 4 Figure 4.21: Mixing zone width L(t) for the simulations with the Thin Cylinder and Large Cylinder configurations (Tab. E).

Figure 4 .

 4 Figure 4.22: Sketch of the Large Cylinder configuration with a twice higher computational domain, giving a cylinder of diameter 2R cyl = 0.194 m and height 2H cyl = 0.388 m.

Chapter 4 .

 4 The Rayleigh-Taylor instability through a grid: Influence of the grid

Figure 4 .

 4 Figure 4.23: Mask function χ G , defining the grid, along the x-axis. The blue colored parts give for which x the mask function is equal to 1. (a) Centered grid and (b) non centered grid.

  k 0 = 470.00 m -1 k c = 1400.00 m -1 k 0,EXP = 1540.00 m -1 k c,EXP = 20480.00 m -1 k g = 3141.59 m -1

Figure

  Figure 4.24: Growth rate Ω(k) (Eq. 4.11) (black continuous line) for simulation parameters (ν = D = 1.2 × 10 -6 m 2 .s -1 and δ = 3 mm) compared to the growth rate Ω EXP (k) for experiment parameters (ν = 700D = 1 × 10 -6 m 2 .s -1 and δ = 0 mm). The red lines mark the critical wavenumber k c and the blue lines give the fastest growing wavenumber k 0 . The smallest wavenumber allowed by the experimental grids k g is given in green dashed line (with λ g = 2 mm for Grid3). The common parameters are A = 0.003, G 0 = 9.81 m.s -2 .

Figure 4 .Figure 4 .

 44 Figure 4.25: (a) Mixing zone widths L as a function of time t for the simulations with the different seed s and N p,XZ . (b) Range of the mixing zone widths L(t, s) taking into account all the simulations with N p,XZ = 512 (blue colored area) and N p,XZ = 1024 (orange colored area) and their ensemble mean (see Eq. (4.25)) in thin lines.

  (series and numbers), number of points in every directions N p X = N p Y = N p Z and initial perturbation wavelentgh λ peak . The 3D simulations have a viscosity of 1.2 ×10 -6 m 2 .s -1 . The grid's parameters in these simulations are: thread diameter l = 17 mm and mesh size d = 34 mm giving a porosity of ϕ = 0.44. The DNSConv3DA and DNSConv3DC have threads only in the X direction while the DNSConv3DB and DNSConv3DD have threads in both X and Y directions.

Figure 4 .

 4 Figure 4.27: Slice along the XZ plane of the instantaneous concentration fields at t = 4 s for (a) DNSConv3DA1, (b) DNSConv3DB1, (c) DNSConv2DC1, (d) DNSConv3DD1, (e) DNSConv3DA2, (f) DNSConv3DB2, (g) DNSConv3DC2 and (h) DNSConv3DD2 of Tab. G The heavier fluid is in red, the lighter one in blue. The black parts show the grid. The upper walls are not visible.

Figure 4 .

 4 Figure 4.28: Growth rate Ω(k) (Eq.(4.11)) along with the most unstable wavelength λ 0 (dashed red line) and the critical wavelength (continuous red line). Both the initial wavelength λ peak of Tab. G are given in dashed blue and orange lines. The mesh size is also indicated in dashed green line.

Figure 4 .Figure 4 . 30 :

 4430 Figure 4.29: Mixing zone width L as a function of time t for the simulations DNSConv3D of Tab. G.

Figure 4 .

 4 Figure 4.31: Mixing zone width L as a function of time t for the 2D periodic simulations with one thread of width 17 mm.

Chapter 4 .Figure 4 . 32 :

 4432 Figure 4.32: Schematics showing on a 2D slice of a three-dimensional periodic cubic simulation at t ∼ 4 s with a grid at the interface of mesh size d = 10 mm and thread diameter of l = 10 mm, giving a porosity of ϕ = 0.25. In hatched white the points taken into account when computing the horizontal mean concentration. (a) First method giving ⟨C⟩ XY,tot . (b) Second method giving ⟨C⟩ XY,wg . (c) Third method giving ⟨C⟩ XY,wgc .

Figure 4 .

 4 Figure 4.33: Horizontal mean concentration profiles ⟨C⟩ XY,tot , ⟨C⟩ XY,wg and ⟨C⟩ XY,wgc at time t ∼ 5 s for a 3D periodic cubic simulation with a grid at the interface of mesh size d = 10 mm and thread diameter of l = 10 mm giving a porosity of ϕ = 0.25.

Figure 4 . 34 :

 434 Figure 4.34: Elementary configuration considered in the horizontal mean concentration model.

  Figure 4.35: Horizontal mean concentration profiles. In blue, the horizontal mean concentration profile ⟨C⟩ XY,wgc computed without taking into account the columns above and below each thread of the grid.In green, the total horizontal mean concentration profile ⟨C⟩ XY,tot computed by taking into account every points. In red, the total horizontal mean concentration profile ⟨C⟩ XY,tot computed with the model Eq. (4.28).

Figure 4 .

 4 Figure 4.37: 2D slices along the XZ plane of the instantaneous concentration fields at time t ∼ 9 s of (a) DNSGA2, (b) DNSGB3 (ϕ = 0.83), (c) DNSGB2 (ϕ = 0.6) and (d) DNSGB1 (ϕ = 0.44) in Tab. H.

Figure 4 .

 4 Figure 4.38: (a) Integral scale L u (t) (Eq. (4.29)) and (b) Dimensionality parameter sin 2 (γ)(t) (Gréa [2013]). Plots made for the simulations DNSGA2, DNSGB1, DNSGB2 and DNSGB3 (Tab. H).

Chapter 4 .Figure 4 .

 44 Figure 4.39: (a) and (b) Horizontal mean concentration profiles ⟨C⟩ XY,wgc (z) at times t ∼ 5 s and t ∼ 9 s respectively. The gray area represents the height of the largest threads used in the simulations considered (height = l = 10 mm). The dashed line represent the value of the arbitrary threshold use for the profile method (ϵ = 0.03). (c) Mixing zone width L(t) and (d) growth rate α(t) computed with the integral method and the profiles method using ⟨C⟩ XY,wgc . The gray region shows the times onto which the growth rate was averaged to give the values in the legend. Plots made for the simulations DNSGA2, DNSGB1, DNSGB2 and DNSGB3 (Tab. H).

Figure 4 .

 4 Figure 4.40: (a) Mixing zone width L(t) and (b) growth rate α(t) computed with the integral method and the profiles method using ⟨C⟩ XY,tot . The gray region shows the times onto which the growth rate was averaged to give the values in the legend. Plots made for the simulations DNSGA2, DNSGB1, DNSGB2 and DNSGB3 (Tab. H).

4. 3 . 4 . 3 . 2 . 2

 34322 Numerical study of the influence of the grid Simulations with k g > k 0

Chapter 4 .Figure 4 .Figure 4 .

 444 Figure 4.42: 2D slices along the XZ plane of the instantaneous concentration fields at time t ∼ 9 s of (a) DNSGA1, (b) DNSGD3 (ϕ = 0.83), (c) DNSGD2 (ϕ = 0.6) and (d) DNSGD1 (ϕ = 0.44) in Tab. H with λ peak = 13 mm.

Figure 4 .

 4 Figure 4.44: Mixing zone width L(t) computed with the integral method and the profiles method using ⟨C⟩ XY,tot for (a) λ peak,1 and (b) λ peak,2 . The gray region shows the times onto which the growth rate was averaged. Plots made for the simulations DNSGA1, DNSGD1, DNSGD2 and DNSGD3 (Tab. H).

Chapter 4 .Figure 4 .Figure 4 .

 444 Figure 4.45: Growth rate α(t) computed with the integral method and the profiles method using ⟨C⟩ XY,tot for (a) λ peak,1 and (b) λ peak,2 . The gray region shows the times onto which the growth rate was averaged. Plots made for the simulations DNSGA1, DNSGD1, DNSGD2 and DNSGD3 (Tab. H).

Figure 4 .

 4 Figure 4.47 shows the slices of the instantaneous concentration fields for DNSGC1, DNSGD1 (with λ peak,1 ) and DNSGE1.

Figure 4 .

 4 Figure 4.47: 2D slices along the XZ plane of the instantaneous concentration fields at time t ∼ 5 s of (a) DNSGC1 (d = λ 0 = 13 mm), (b) DNSGD1 (d = 8 mm) and (c) DNSGE1 (d = 7 mm) in Tab. H with λ peak = 13 mm. The heavier fluid is in red and the lighter one in blue. The walls and grid are shown in black.

Figure 5

 5 Figure 5.1: Instantaneous concentration field beneath the grid in the plane of the laser sheet at time t ∼ 7 s for EXPa7 (with A = 0.0037, grid mesh size d = 1 mm and a porosity ϕ = 0.41) with the mean concentration ⟨C⟩ Z (in red) as a function of the horizontal coordinate X.

kFigure 5 . 2 :

 52 Figure 5.2: Growth rate Ω as a function of the wavenumber k (Eq. (4.11)). The k peak of the DNSλ and DNSTC series of Tab. A are given in dashed lines of color blue, orange, green, red and purple for DNSλ1, DNSλ2, DNSλ3, DNSλ4 and DNSTC respectively. The red continuous line indicates the wavenumber corresponding to the diameter of the cylinder k D cyl .

Figure 5 . 4 :

 54 Figure 5.4: (a) Azimuthal mean concentration profiles ⟨C⟩ θZ (r, t) at time t = 4.5 s, (b) mixing zone width L(t), (c) the mixing parameter Θ(t) and (d) the growth rate α(t) for simulations DNSλ and DNSTC (Tab. A). The insert in (b) is a zoom on the earlier times of the mixing zone width L(t). The dotted gray line in (c) and (d) represents the time from which the mixing parameter and the growth rate mean values, given in the legend, are computed.

  Figure 5.5: (a) Instantaneous 3D concentration fields and (b) 2D vertical slices along the XZ plane for the configuration Thin Cylinder (height 2H cyl = 60 cm and diameter 2R cyl = 19.4 cm) with an initial wavelength of λ = 0.13 m at time t ∼ 4 s. (c) Shows the 2D vertical slice at t ∼ 6 s. The heavier fluid is in red and the lighter one in blue.

  Figure 5.7: (a) Azimuthal mean concentration profiles ⟨C⟩ θZ (r, t) at time t = 4.5 s and (b) mixing zone width L(t) for simulations DNSCurv (see Tab. B).

  Figure 5.8: (a) Sketch of the wall opening with time (top view). (b) Zoom in on the interface in experiment EXPa6 (Tab. B) during the contact of both tanks. It shows around half of the surface of the grid of heavier fluid in contact with the lighter fluid. The initial contact was made on the top right of the image.

  5.10a), the mixing zone widths (Fig.5.10b), the mix-

  Figure 5.10: (a) Azimuthal mean concentration profiles ⟨C⟩ θZ (r, t) at time t = 4.5 s, (b) mixing zone width L(t), (c) mixing rate Θ(t) and (d) vertical kinetic energy profiles ⟨v 2 z ⟩ θZ (r, t) at time t = 4.5 s for simulations DNSContact (Tab. C).

Figure 5 . 11 :

 511 Figure 5.11: Imposed vertical velocity u z at the interface z = 0 as a function of x with the parameter r 0 indicated.

  Figure 5.13: (a) Azimuthal mean concentration profiles ⟨C⟩ θZ (r, t) at time t = 4.5 s, (b) mixing zone width L(t), (c) concentration variance ⟨c 2 ⟩ XY (z, t) at time t = 4.5 s and (d) vertical kinetic energy profiles ⟨v 2 z ⟩ θZ (r, t) at time t = 4.5 s for simulations DNSInitU (see Tab. D).

Figure 5

 5 Figure5.14: Effect of the grid meniscii with the suction case with P ≤ P a . The black dots represents the grid and the red dashed lines the different meniscii. U is the suction velocity.

Figure 5

 5 Figure5.15: Sketch of an experiment in the suction configuration (P a ≥ P) where the meniscii are oriented toward the upper cylinder. It also shows the development of the initial interface between the two fluids. The red dashed line represents either the fresh water-air interface or the initial fresh water-salt water interface, the blue line shows the air still traped between the grid and the fresh water being evacuated, and the green dashed line represents the free surface in the upper cylinder.

Figure 5

 5 Figure 5.16: (a) Displacement of the free surface h(t) and (b) velocities U(t) and V(t) as a function of time for a case with A = 0.003, h 0 ∼ 26 cm and U 0 = 0.54 m/s.

  Figure5.17: (a) Case I in 2D with the wall opening at the interface at speed v PW . (b) Case II in 2D with the wall opening at the interface at speed v PW and the top and bottom walls moving upwards at speed V(t). (c) Case III in 2D with the wall opening at the interface at speed v PW , the top and bottom walls moving upwards at speed V(t) and a layer at the interface (hatched red) in which the upward jet velocity U(t) is forced. The gray hatched areas represents the walls defined via the penalisation method χ and χ u .

  Figure 5.19: Vertical slices of the concentration field along the XZ plane at time t ∼ 1 s for two simulations with U 0 = 4.5 cm/s and (a) t o = 0.1 s or (b) t o = 0.05 s.

Figure 5

 5 Figure 5.20: Diagrams of the simulations performed with the aspiration model and the grid at the interface organized by value of the porosity: (a) ϕ = 0.6 and (b) ϕ = 0.44.

Figure 5

 5 Figure5.21: Growth rate Ω(k) given by the linear stability analysis with the most unstable wavelength λ 0 (in red dashed line) and the critical wavelength λ c (in continuous red line). The grid mesh size d is also given (in dashed blue line) along with the initial perturbation wavelengths λ p studied in the DNSλ p simulations (Tab. H). The gray hatched area are the wavelengths that cannot develop due to the grid.

5. 4 .

 4 Fig. 5.24a with R 0 = 14 Addition of a grid on the aspiration modelSeries N • ϕ d [mm] l [mm] λ p [mm] R 0

  Figure5.24: Initial slight opening (in blue) in the wall at the interface (in gray) for (a) an initial contact radius R 0,1 = 14 mm and (b) an initial contact radius R 0,2 = 42 mm.

Figure 5

 5 Figure5.27: Growth rate Ω(k) given by the linear stability analysis with the most unstable wavelength λ 0 (in red dashed line) and the critical wavelength λ c (in continuous red line). The grid mesh sizes d, studied in the DNSd simulations (Tab. J), are also given (in dashed blue line and dashed green lines).

  Figure 5.28: Vertical slices of the concentration field along the XZ plane at t ∼ 5 s for (a) DNSd1 with d = 20 mm and (b) DNSd2 with d = 8 mm (Tab. J). The heavier fluid is in red and the lighter one in blue. The walls are in black.

  Figure 5.29: Vertical slices of the concentration field along the XZ plane at t ∼ 5 s for (a) DNSϕ1 with ϕ = 0.44 and (b) DNSϕ2 with ϕ = 0.6 (Tab. K). The heavier fluid is in red and the lighter one in blue. The walls are in black.

  t = 0s t = 1s t = 2s t = 3s Figure A.1: Vertical slices of the concentration field along the XZ plane at four different times t for Case I. The heavier fluid is in red and the lighter one in blue.

Table A :

 A Label (series and number), Atwood number A, forcing parameter F and frequency ω considered for the experiments in this work. The wavenumbers k n,m and mode type corresponding to the primary Faraday wave are also indicated. The initial interface thickness δ is either measured by a probe when available or directly from the camera (labelled with * ).

	Chapter 2.
		015 0.30	2.555	13.28	(4, 0)	1.3 *
		0.015 0.30	2.953	16.60	(5, 0)	0.8
		0.015 0.30	3.462	23.25	(7, 0)	1.0
		0.015 0.30	3.924	33.21	(9, 0)	1.1
		0.015 0.40	2.457	13.28	(4,0)	1.7
		0.015 0.40	3.462	23.25	(7,0)	0.6
		0.015 0.40	3.942	29.89	(9,0)	0.6
		0.015 0.50	3.142	19.93	(6,0)	1.4
		0.015 0.50	3.924	26.57	(8,0)	0.5
	EXPb	0.03 0.30	4.524	19.93	(6,0)	0.6 *
		0.03 0.40	3.462	9.96	(3,0)	0.8
		0.03 0.50	4.524	19.93	(6,0)	0.8
		0.03 0.69	4.290	16.60	(5,0)	0.4 *
		0.03 0.69	4.290	16.60	(5,0)	0.5 *
		0.03 0.69	4.290	19.93	(6,0)	0.5 *
	EXPc	0.045 0.50	3.066	6.64	(2,0)	1.2
		0.045 0.69	4.290	13.28	(4,0)	0.3 *
		0.045 0.69	4.290	13.28	(4,0)	0.7 *

The subcritical transition to turbulence of Faraday waves in miscible fluids

Table B :

 B Label (series and number) and parameters in physical units (Atwood number A, forcing parameter F and frequency ω) taken for the direct numerical simulations presented in this work. The cases DNSa, d and e correspond to the wavebreaking detection. The series DNSb and c are dedicated to the competition between mode (2, 0) and

even 2.3. Mode selection mechanism of the Faraday wave

	1.0
	0.8
	0.6
	F
	0.4				(2,0)				(3,0)
	0.2
	1.6 0.0	1.7	1.8	1.9	2.0	2.1	2.2	2.3	2.4	2.5
					ω/ AG 0 k 2,0

Table B :

 B Label (series and number), mesh size d in millimeters, grid porosity ϕ, Atwood number

	EXPa	1	1	0.41 0.0035	1003.1	996.1	✓	1.07
		2	1	0.41 0.0037	1003.3	996.0	✓	1.07
		3	1	0.41 0.0037	1003.4	996.1	✓	1.18
		4	1	0.41 0.0037	1003.4	996.0	✓	1.76
		5	1	0.41 0.0037	1003.4	996.0	✓	0.8
		6	1	0.41 0.0037	1003.4	996.0	✓	1.76
		7	1	0.41 0.0037	1007	999.5	✗
		8	1	0.41 0.0038	1003.5	995.9	✓	1.14
		9	1	0.41 0.0038	1003.5	995.9	✓	1.58
		10	1	0.41 0.0075	1008.1	993.1	✓	2.82
	EXPb	1	1.8	0.48 0.0015	1001.9	998.9	✓	1.91
		2	1.8	0.48 0.0024	1002.9	998.0	✗
		3	1.8	0.48 0.0037	1003.4	996.0	✓	2.01
		4	1.8	0.48 0.0037	1003.6	996.2	✓	1.91
		5	1.8	0.48 0.0037	1003.4	996.0	✗
		6	1.8	0.48 0.0038	1003.4	995.9	✓	1.41
		7	1.8	0.48 0.0038	1003.4	995.9	✓	2.01
		8	1.8	0.48 0.0038	1003.5	996.0	✓	1.91
		9	1.8	0.48 0.0038	1003.3	995.8	✗
		10	1.8	0.48 0.0047	1003.8	994.5	✓	2.82
		11	1.8	0.48 0.0047	1004.1	994.7	✗
		12	1.8	0.48 0.0048	1004.1	994.6	✗
		13	1.8	0.48 0.0060	1007.3	995.2	✗
		14	1.8	0.48 0.0073	1007.9	993.3	✗
		15	1.8	0.48 0.0074	1007.8	993.0	✓	2.35
		16	1.8	0.48 0.0134	1015.1	988.2	✓	2.82
	EXPc	1	2	0.48 0.0009	1001.1	999.2	✓	0.83
		2	2	0.48 0.0036	1003.4	996.3	✓	1.76
		3	2	0.48 0.0123	1015.0	990.4	✓	2.52

  2, in Jause-Labert et al. [2012] and[START_REF] Kadoch | A volume penalization method for incompressible flows and scalar advection-diffusion with moving obstacles[END_REF]. It consists of an additional force in the momentum equation (see Eq. (4.21)) and an effective diffusivity in the advection-diffusion equation (see Eq.(4.22)).

Table E :

 E Label, diameter of the cylinder D cyl , height of both cylinder tanks 2H cyl and initial wavenumber k peak are given in this table for the two simulations considered. The simulations have an Atwood number of 0.003 and 512 3 points.

	t = 0s	t = 0s	t = 0s	t = 0s
	t = 2.0s	t = 2.0s	t = 2.0s	t = 2.0s

t = 4.0s t = 4.0s t = 4.0s t = 4.0s t = 6.0s

Table G :

 G Label

	DNSConv3DA 1	512	11
		2	1024	11
	DNSConv3DB	1	512	11
		2	1024	11
	DNSConv3DC	1	512	5
		2	1024	5
	DNSConv3DD 1	512	5
		2	1024	5

Table H :

 H Label (series and numbers), initial wavelength of the perturbation λ peak , mesh size d of the grid, thread diameter l of the grid and porosity ϕ for the simulations considered. The simulations have an Atwood number of

		1	5	20	10	0.44
		2	5	20	6	0.6
		3	5	20	2	0.83
	DNSGC	1	13	13	7	0.44
	DNSGD 1	13/5	8	4	0.44
		2	13/5	8	2	0.6
		3	13/5	8	0.8	0.83
	DNSGE	1	13	7	3	0.44

5.1. The cause of the apparition of the central jet

	Series	N • Configuration k peak [m -1 ] λ peak [m]
	DNSλ	1	LC	292.95	0.021
		2	LC	585.90	0.011
		3	LC	1318.26	0.005
		4	LC	1757.69	0.004
	DNSTC		TC	47.35	0.13

Table A :

 A Label (series and numbers), configuration Large Cylinder (LC) or Thin Cylinder (TC) initial unstable wavenumbers k peak and the corresponding initial unstable wavelength λ peak are given for the five simulations considered.

Table C :

 C Label (Series and number), initial unstable wavenumber k peak , curvature of the interface C int , the aperture time t o and the aperture velocity v PW for the four simulations considered.

	DNSContact 1	585.90	0	0.5	0.194
	2	585.90	0	2	0.0485
	3	-	1	0.5	0.194
	4	-	1	2	0.0485

R(t) R(t + δt)

Table E :

 E Label (series and numbers), initial jet velocity U 0 , radius of curvature of the meniscii r, pressure jump P a -P and contact time t o for the simulations of Case II with different (U 0 ,t o ) couples.

		2		3.6	4	0.5 0.1
	DNSU0toB	1				0.5
		2	4.5	14.4	1	0.1
		3				0.05
	DNSU0toC	1 2	3	32.4	0.4	0.5 0.1

Table G :

 G Parameters of the 2D simulations with the number of points, Atwood number A, viscosity and diffusion coefficient ν = D, scalar penalisation parameters η c , aperture time t o , initial jet velocity U 0 and filter parameter C s .

Table J :

 J N • ϕ d [mm] l [mm] λ p [mm] R 0 [mm]Label (series and numbers), porosity ϕ, grid mesh size d, grid thread diameter l, initial perturbation wavelentgh λ p and initial contact radius R 0 for the 2 simulations considered here.

	DNSd 1 0.44	20	10	13	42
	2 0.44	8	4	13	42

FARAMIX2 campaign made at GTT (Gaztransport & Technigaz company) on October

with B.-J. Gréa (CEA), A. Briard (CEA) and L. Gostiaux (CNRS, LMFA)

2.6. Conclusion

Acknowledgements

lations and thus, no central jet appears even in the slightest, we give the mean azimuthal concentration profile ⟨C⟩ θZ as a function of the radius r of the cylinder in Fig. 5.4a. This is shown for all four simulations presented in Tab. A at time t = 4.5 s. We notice that all of the profiles stay roughly around 0.5 throughout the length of the cylinder since there is the same amount of heavy and light fluids at any given r. This stays true at all time during the simulations. Thus, the profiles and the instantaneous concentration fields demonstrate that neither the symmetry nor the homogeneity of the problem are broken by the change in the initial perturbation wavelengths.

To better understand the dynamics of the instability in each simulation, the mixing zone L and its growth rate α, defined previously in Chapter 3 (Eqs. (3.1) and (3.2) respectively), and recalled here for clarity: L(t) = 2αAGt 2 and α = L2 8AGL , are investigated in Fig. 5.4 along with the mixing parameter Θ, defined as:

Mixing parameter (see [START_REF] Youngs | Numerical simulation of mixing by Rayleigh-Taylor and Richtmyer-Meshkov instabilities[END_REF])

where 1 L ⟨c ′2 ⟩ XY dZ is the concentration variance, and c ′ = C -⟨C⟩ XY .

Chapter 5. The Rayleigh-Taylor instability through a grid: Central jet t = 0.1s t = 0.1s t = 0.1s t = 0.1s t = 1.5s t = 1.5s t = 1.5s t = 1.5s t = 3.0s t = 3.0s t = 3.0s t = 3.0s t = 4.5s the lighter fluid goes upward on the sides of the cylinder and the heavier fluid goes down in the center, following the curvature of the interface. This effect is particularly visible for DNSCurv2 with C int = 5 mm but it is also present for DNSCurv1 with C int = 1 mm, which is the curvature of the grids in the experiments. Hence, the curvature has an effect on the development of the instability and it does break the horizontal homogeneity. However, this effect is contrary to the one we observed in the experiments where the heavier fluid would fall down on the sides of the cylinder and the lighter fluid would go up in the center. Moreover, as seen in the previous section, the smaller k peak the smaller Ω. As implementing a curvature at the interface can be seen as using particular and very small initial perturbation wavenumber, the effect of the curvature develops over a long period of time wheres the jet develops rapidly in the experiments.

This effect of the curvature on the instability development is visible on the azimuthal mean concentration profiles plotted in Fig. 5.7a. Indeed, the mean concentration is smaller on the sides (r = R cyl ) than in the center of the cylinder (r = 0), meaning that, as shown in the instantaneous concentration fields, the heavier fluid falls down in the middle and the lighter fluid goes up on the sides. For the case with C int = 1 mm, this effect is only just slightly visible as the mean concentration stays very close to 0.5 whereas it is very clearly visible for C int = 5 mm with a mean concentration that reaches 0.8 in the center of the cylinder. The curvature has a significant impact on the mixing zone growth in Fig. 5.7b as L increases much more rapidly for C int = 5 mm than for C int = 1 mm. Note that the mixing zone width is initially bigger for DNSCurv2 even if the same parameters were used: this means that the curvature artificially increases the mixing zone width due to the breaking of the horizontal homogeneity. The azimuthal mean concentration profiles show the effect of the curvature of the interface for DNSContact3 and DNSContact4 as they are close to the ones observed for DNSCurv of Sec. 5.1.2. The effect of the aperture is also visible on the profile of DNSContact4 as the mean concentration is higher in the center (r = 0) and smaller on the sides (r = R cyl ) than for DNSCon-tact3. This comes from the fact that the instability had developed in the center before the sides and thus the sides have a delay to catch up with. The profiles of the DNSContact1 and DNSCon-tact2 oscillate around 0.5 with no sign of a central jet appearing.

The cause of the apparition of the central jet

The delay in contact between the fluids also influences the mixing zone width. Indeed, its growth is delayed when comparing DNSContact1 and DNSContact2 or DNSContact3 and DNSContact4 in Fig. 5.10b which was rather expected. Moreover, it can be observed that this delay has a much greater effect on the dynamics of the Rayleigh-Taylor instability when it is coupled with a curved interface. Indeed, in Fig. 5.10d, the vertical kinetic energy profiles are rather flat and similar for DNSContact1 and DNSContact2, whereas strong peaks appear for Chapter 5. The Rayleigh-Taylor instability through a grid: Central jet horizontal homogeneity and the symmetry with respect to the z = 0 plane are visibly broken in this simulation. This is very close to what was observed in the experiments: the ascending central jet is created. The amplitude of this effect is very dependent on the initial speed u 0 . Indeed, the DNSInitU2 case having a much smaller initial speed than the DNSInitU1 case, the fall of the heavier fluid is less visible and it takes much more time to come into effect. So for small initial speed, the effect is not instantaneous, contrary to the experiments, and develops only after the triggering of the Rayleigh-Taylor instability. This effect is corroborated by the azimuthal mean concentration profiles shown for the three DNSInitU simulations in Fig. 5.13a. Moreover, as can be expected from the instantaneous concentration fields, the mean concentration profile shows less deformation and stays closer to 0.5 for DNSInitU2 than for DNSInitU1. This initial speed u 0 also influences the growth of the mixing zone which increases much more rapidly for the larger u 0 in Fig. 5.13b. Another effect of this difference of initial speed, is the fact that the concentration variance profile gets wider as the initial speed is increased, which is visible when comparing DNSInitU1 and DNSInitU2 in Chapter 5. The Rayleigh-Taylor instability through a grid: Central jet Indeed, taking ρ h = 1003.4 kg.m -3 , U 0 = 0.54 m.s -1 and t o = 0.5 s, we get t P (U = 0) ∼ 0.34 s, so the process ends before the contact is fully made between the two fluids. The equation (5.15) is a general expression that stays true even if t P > t o .

Displacement of the free surface and its velocity

Having the general expression for the jet velocity, we can calculate the displacement h(t) -h 0 of the free surface due to this jet and the velocity V(t) at which this free surface is displaced. Indeed, we know that R 2 cyl V(t) = R(t) 2 U(t), hence:

(5.17) Thus, we have:

(5.18)

We also know that h

(5.19)

In Fig. 5.16 is represented the displacement of the free surface and the velocities U(t) and V(t) as a function of time t in the upper cylinder of 30 cm of height taking A = 0.003, h 0 = 0.26 m and U 0 = 0.54 m.s -1 . In this case, the process ends before the contact is fully made, thus only t ≤ t o is shown. The displacement is shown in dashed red and the free surface position with the air above it is shown in green. We see that in this case the free surface rises of 1.5 cm.

Chapter 5. The Rayleigh-Taylor instability through a grid: Central jet

The influence of the mesh size d

Now that the effect of the initial contact radius has been investigated, it can be fixed at R 0 = 42 mm as it is the value for which the simulations resemble the most the experimental observations. The grid mesh size and thread diameter are the next parameters to be studied by comparing two simulations with different d and l but the same other parameters summarized in Tab. J.

In Fig. 5.27 we show the grid mesh sizes considered here compared to the most unstable wavelength and the critical wavelength, given by the linear stability analysis. These two mesh

A.1.2 Case II: Free surface velocity

In addition of Case I, Case II simulates the displacement of the free surface caused by the pressure jump of the aspiration model. The value of U 0 chosen is 45 cm/s, and gives a radius of curvature of 1.5 mm which is a mean value between the largest and smallest mesh size of the experimental grids used. It is a large value not very realistic but it gives an idea of the most extreme effect that could develop with just the displacement of the free surface. As a reminder, the free surface velocity is given by V(t) = U(t)R 2 (t)/R 2 cyl (Eq. (5.18)). Notice how, on the vertical slice of Fig. A.2, the upper interface between the lighter and heavier fluid is moved to the bottom of the domain at t = 1 s. As previously stated this displacement is particularly important due to the choice of an extreme value of U 0 . This displacement creates an ascending central jet at the center of the cylinder which reaches the upper boundary before the first second of the simulation has passed. For similar Atwood number, the lateral descending jets in the experiments have a velocity of approximately 1.6 cm/s, whereas in this case they have a velocity of around 2.7 cm/s. This shows that the effect created with this initial velocity is much more violent than in the experiments and thus, here we will need to take a smaller velocity U 0 .

A.1.3 Case III: Jet velocity

For this case, we add a penalisation layer at the interface, inside the aperture of the opening wall, in which the jet velocity U(t) is imposed. The same U 0 is used as in Case II.

A.2. Convergence study

The vertical slices along the XZ plane of the instantaneous concentration fields are given at different times in Fig. A.3. Compared to Case II, there are no descending jets on the sides at t = 1 s, they appear later on (see image at t = 2 s). The central jet is also thinner in this case, thus the descending jet are actually much closer to the center than in Case II. This is not what is observed in the experiments where the descending jets are most often falling closer to the side walls of the cylindrical tank.

As a conclusion, implementing only the free surface velocity of the aspiration model is enough to see the apparition of a central ascending jet. Imposing the jet velocity U 0 changes the descending jets position which gives an effect further from the experiments. Thus, Case II of the aspiration model will be implemented in the future simulations.

A.2 Convergence study

In the following, Case II is chosen as the reference configuration to create a central ascending jet and descending jets close to the cylinder walls. Before going on with better resolved simulations, it is necessary to determine optimal parameters in the penalisation layers to ensure that the imposed free surface velocity V is indeed reproduced by the simulation. The numerical parameters to consider are the filtering intensity of the mask C s , the penalisation intensity η u , the initial contact radius R 0 , the form of the contact radius R(t), and possibly the viscosity ν. Both 2D and 3D simulations are performed for this convergence study.

A.2.1 2D simulations

Using 2D simulations, the role of the viscosity and the filter parameter C s , used in the Gaussian filter expression exp(-C s (k/N p ) 2 ) as a means to avoid Gibbs oscillations in the mask function, is investigated (Jause-Labert [2012]). The parameters of the Case II simulations are given in Tab. B. It was seen in the previous section on the implementation of the aspiration model that U 0 = 45 cm/s for the initial jet velocity is an unrealistic choice and gives a too strong jet. Thus, a smaller value should be chosen for this convergence study. Hence, the initial velocity taken in these simulations is U 0 = 4.5 cm/s. This gives much more realistic values of the curvature radius and the jet still exits. The effect of the resolution has also been tested but for the sake of brevity, the 512 2 simulations are not shown as it was found that they yield very similar results compared to the 1024 2 ones.

Viscosity: Two simulations are performed with ν = 2.3 × 10 -6 m 2 .s -1 and ν = 0.6 × 10 -6 m 2 .s -1 . The vertical velocities inside the upper wall of both these simulations are compared to the velocity imposed by the aspiration model given in Eq. (A.1) on Fig.

In the simulations considered here, P a -P = 0 before t o is reached. Thus, the free surface velocity is defined by V th (t) = U(t) t 2 t 2 o before t o and V(t) = 0 after. The vertical velocities given by the simulations are perfectly equal through time. Moreover, the time evolution of the imposed velocity is well recovered by the simulations, despite a slight difference of about 3 × 10 -5 m/s around the maximum value. Thus, the viscosity has no influence on the convergence of the velocity.

Filter parameter:

The filter parameter C s is varied as we suspect that the slight gap observed comes from the fact that for large values of C s , the mask 'leaks' (less sharply defined). As the filter parameter is decreased, the mask becomes more 'abrupt'. However, for too small values of C s , Gibbs oscillations appear. Hence, a compromise has to be found. Three simulations are performed with three different filter parameter: C s = 256, C s = 128 and C s = 64. Again we compare the vertical velocities given by each simulation to the theoretical one in Fig. A.4b. Note that, when the filter parameter is decreased, we get closer to the theoretical value of the vertical velocity. This difference being very small and the filter parameter C s = 64 being a bit small to efficiently remove all Gibbs oscillations for this mesh resolution, C s = 128 will be the filter parameter used hereafter.

A.2.2 3D simulations

The gap between the imposed velocity and the velocity given by the simulations persisting, 3D simulations are performed to check if that is not a 2D vs 3D effect. These 3D simulations are used to understand the effects of the penalisation parameters, the form of the contact radius R(t) and the initial contact radius R 0 . The parameter η u can be seen as a relaxation parameter and it is expected to greatly influence how close to the theoretical input value the velocity inside a penalised layer is. Thus, depending on the value of this parameter, the form and magnitude A.2. Convergence study of V(t) will be changed, which will affect the form and magnitude of the jet resulting from the aspiration visible on the concentration field.

128

Table C: Parameters of the 2D simulations with the number of points, Atwood number A, viscosity and diffusion coefficient ν = D, scalar penalisation parameters η c , aperture time t o , initial jet velocity U 0 and filter parameter C s .

As very similar results were obtained for simulations with 512 2 points compared to simulations 1024 2 points, 3D simulations are performed with 512 3 points. Moreover, we verified that the aperture time has no influence on the convergence of the simulation. Hence, as the time step will be greatly decreased to study the effect of smaller penalisation parameter, the aperture time is taken smaller than in the previous simulations (t o = 0.1 s instead of t o = 0.5 s).

Penalisation parameter and form of the contact radius: Seven simulations given in Tab. D are performed in which the penalisation parameter η u is varied with R 0 = 7 mm.

Table D: Label (series and numbers) and velocity penalisation parameter η u for the Case II simulations with different velocity penalisation parameter. The parameters of the simulations are given in Tab. C with R 0 = 7 mm.

The ones marked η u,1 /η u,2 in Tab. D are simulations that were made using η u = η u,1 up to t = t o + 0.5 s and then η u = η u,2 . This allows for faster simulations as the time step is increased after t o + 0.5 s. It was verified by comparing simulations DNSη u 3 and DNSη u 4 and simulations DNSη u 5 and DNSη u 6 that it gave the same results than when η u = η u,1 throughout the entire simulation, as seen on the instantaneous concentration fields in Fig. A.5.

Two types of aperture of the wall at the interface were used to check the importance of the form of R(t):

• the usual linear aperture: R(t) = R cyl t t o ,

• an exponential one: R(t) = R cyl 1-e -t/to 1-e -1 . For the exponential aperture, the jet velocity and the free surface velocity are written: and

Once again, the pressure equalisation happens before t o is reached, thus only the expressions for t ≤ t o are of interest to us.

The results are given in Fig. A.6 with the vertical velocity inside the upper wall given as a function of time and compared with V th (t) given by Eqs. (5.18) and (A.3).

Note the very visible difference between the simulations when varying η u up to η u = 1 × 10 -5 s and how far they are from the imposed velocity in both the linear aperture and the exponential one. These differences can also be seen on the instantaneous concentration fields with very different central jets shown on Fig. A.7 for the linear aperture. Thus, the type of aperture (linear or exponential) does not seem to impact the convergence of the simulations. A slight difference can still be seen on the simulations DNSη u 6 and DNSη u 7. However, this difference is small enough to consider that our simulations are converged at η u = 1 × 10 -5 s. The reason for these large discrepancies between the imposed V and the numerical results is that for the largest η u , the time step imposed by the penalisation ∆t η is much larger than the time step ∆t U imposed by U 0 . Hence, stronger penalisation parameters are required to reach ∆t η ∼ ∆t U . Then, it was verified that for a small enough η u (like 1 × 10 -4 ), there is no difference by further imposing an even smaller time step.

Initial contact radius:

To further test out the effects of the initial opening on the convergence, two values of R 0 were used for DNSη u 6 with the linear aperture: R 0,1 = 7 mm and R 0,2 = 42 mm. Figure A.8 shows both the vertical velocity and the vertical slices of the instantaneous concentration field for these simulations. Note that the form and extension of the jet is completely different: indeed, a larger initial opening means more surface through which the lighter fluid can be sucked and thus it gives a larger and less extended jet. However, looking at the vertical velocity inside the upper wall (Fig. A.8a), there is very little difference between the two simulations. This means that, changing the initial opening in the wall will not affect much the convergence of the simulation.

A.2.3 Conclusion

As a summary, for this convergence study, the type of aperture is irrelevant and thus the linear one is chosen for the following. The initial opening has not much effect for the convergence, but its effect is more deeply investigated in Section 5.4. Regarding the penalisation parameters: (i) the filtering parameter C s = 128 is kept as a good compromise; (ii) the penalisation intensity η u is the most critical parameter. To make for less time-consuming simulations, η u = 1 × 10 -5 s is imposed initially up to t o + 0.5 s and then this strong constraint is released by decreasing the penalisation parameter to η u = 5 × 10 -4 s.