Chao-Kun Huang

Keywords: Cavitation, Écoulement diphasique, HEM, TTV Cavitation, Two-phase flow, Homogeneous model, Transport equation model

Turbulence and cavitation : applications in the NSMB and OpenFOAM solvers Résumé L'objectif de ce travail de thèse concerne l'étude et la mise en oeuvre de deux modèles de cavitation dans le solveur NSMB (Navier-Stokes-Multi-Blocks): les modèles HEM (Homogeneous Equilibrium Model) et une équation pour le taux de vide: le modèle à transport de taux de vide (TTV). Le phénom ène de cavitation est modélisé par différentes équations d'état de mélange liquide-vapeur (EOS). De s simulations numériques sont réalisées sur des écoulements diphasiques compressibles unidimensi onnels et bidimensionnels avec des conditions d'interface et comparées à des solutions de référence . De plus, la méthode TTV basée sur le taux de vide incluant les termes source pour la vaporisation et la condensation dans le logiciel libre open source OpenFOAM est également présentée sur la géom étrie Venturi pour capturer le phénomène du jet réentrant. La modélisation de la turbulence joue un r ôle majeur dans la capture des comportements instationnaires et un limiteur est introduit pour réduir e la viscosité turbulente afin de mieux prédire la structure à deux phases. Une comparaison de diver s modèles de cavitation couplés avec des modèles de turbulence est étudiée. Les résultats computat ionnels sont comparés aux données expérimentales existantes.

RÉSUMÉ

En général, la cavitation se réfère à des poches de gaz apparaissant dans un écoulement fluide. En d'autres termes, il s'agit d'un phénomène diphasique avec changement de phase. La cavitation se produit lorsque la pression d'écoulement est inférieure à la pression de vapeur saturante. Les structures ainsi formées sont entraînées par l'écoulement et lorsqu'elles atteignent une zone de pression plus élevée, elles se condensent et implosent violement. La cavitation conduit à des pertes importantes de performance de l'installation, à des problèmes d'instabilités de fonctionnement des machines et à l'erosion des parois du composant. C'est ainsi une source de problèmes techniques primordiaux dans le domaine des turbomachines hydrauliques et de la construction navale. Il existe différents types de cavitation selon la configuration d'écoulement, les propriétés du fluide et les géométries. Généralement, il y a quatre types de cavitation de base et c'est-à-dire traveling cavitation, sheet cavitation, cloud cavitation et tip-vortex cavitation. Il est classique de distinguer si l'écoulement est cavité ou non par le nombre de cavitation qui est défini par l'écart adimensionnel entre une pression de référence et la pression de vapeur saturante, noté σ ∞ = (P ∞ -P vap )/(0.5ρ ∞ U 2 ∞ ). P ∞ représente la pression absolue en un point de référence de l'écoulement, P vap est la pression de la vapeur saturante à la température d'essai, ρ ∞ est la masse volumique du liquide et U ∞ est la vitesse de référence.

La prédiction numérique de la cavitation reste un défi pour plusieurs raisons. La modélisation du changement de phase (thermodynamique) et les interactions avec la turbulence n'est pas encour totalement établie. Du point de vue de la modélisation, la grande majorité des codes dédiés à la simulation de la cavitation est basée sur une approche moyennée à la fois pour l'écoulement diphasique et la turbulence. Une hiérarchie de modèles existe, du modèle simple à trois modèles d'équations (un fluide ou modèle homogène) jusqu'au modèle à sept équations (deux fluides) qui restent plus adaptés pour des géométries simples ou des fluides nonvisqueux. Les modèles deux fluids à sept équations sont les plus complets. Dans ce modèle, on suppose que les deux phases coexistent à chaque point du champ d'écoulement et sont exprimées en termes de deux ensembles d'équations de conservation qui développent l'équilibre de masse, de moment et d'énergie pour chaque phase. L'équation de transport pour la fraction de vide est introduite pour décrire la topologie de l'écoulement. Les modèles réduits à six équations sont similaires aux modèles de sept équations à l'exception sans tenir compte de l'équation d'évolution de la fraction de vide. Cependant, ils restent difficile à utiliser en écoulements industruels (turbomachines). La méthode à un fluide, ou méthode homogène, considère les écoulements comme un mélange de deux fluides se comportant comme un fluide qui est semblable au courant monophasé. De cette façon, un seul ensemble d'équations de conservation est employé pour exprimer l'interaction fluide pour le mélange. Compte tenu de sa simplicité et de son faible coût de calcul, la méthode homogène est plus intéressante pour les simulations numériques des écoulements cavitants.

La plupart des phénomènes de cavitation impliquent une turbulence et l'interaction turbulencecavitation est un phénomène sous-connu et documenté (dû notamment à la difficulté d'effectuer v des mesures expérimentales dans les écoulements cavitants). Les effets de la compressibilité sur la turbulence et les effets de la phase dispersée sont également inconnus. La précision numérique de la cavitation turbulente dépend de la modélisation de la cavitation et de la turbulence. Ainsi, le choix d'une modélisation de la turbulence est une question importante pour la simulation de la cavitation. La simulation numérique directe (Direct Numerical Simulation (DNS)) a la capacité la plus élevée de résoudre toutes les échelles de turbulence. Toutefois, il nécessite une résolution de grille très fine et, par conséquent, il est encore assez difficile à appliquer en raison de la consommation élevée de performances informatiques. Bien que la simulation des grands échelles (Large Eddy Simulation (LES)) ait déjà été mise en oeuvre pour les écoulements turbulents de cavitation, les codes habituels sont formulés dans un modèle de Navier-Stokes (RANS) de Reynolds à tensor turbulent par une équation de transport kε (hypothèse de Boussinesq) Entre l'effort de calcul et la précision. Cependant, les modèles standards de viscosité par tourbillons basés sur l'hypothèse de Boussinesq tendent à sur-prédire la viscosité par tourbillonnement qui réduit l'effet du jet re-entrant et de la décomposition de structure biphasée. Ces modèles de turbulence sont inadéquats pour prédire correctement la dynamique des bulles de cavitation. Plusieurs solutions ont été proposées et testées pour réduire la viscosité des turbulences et améliorer le comportement des modèles de turbulence. Reboud a proposé une modification arbitraire en introduisant un limiteur de viscosité de turbulence assigné en fonction de la densité au lieu d'utiliser directement la densité du mélange. Une méthode basée sur le filtre (Filter-based Method (FBM)) qui combine le concept de filtre et le modèle RANS a été étudiée en imposant une échelle de filtre indépendante, généralement la taille de la grille, sur le calcul de la viscosité de Foucault. Une fois que l'échelle de longueur de turbulence est supérieure à la taille du filtre, la viscosité de turbulence peut être réduite par une fonction de filtrage linéaire. L'interaction entre la turbulence et la cavitation en ce qui concerne l'instabilité et la structure du flux est complexe et mal comprise. De plus, il ya moins d'études sur l'influence des modèles de turbulence sur le débit de cavitation. Dans cette étude, la correction de Reboud est mise en oeuvre en trois modèles de turbulence différents et simulée avec différents modèles de cavitation. L'objectif final est de fournir un aperçu de l'interaction entre les modèles de turbulence et de cavitation.

Cette étude présente la mise en oeuvre et la validation des modèles de cavitation développés au LEGI (Laboratoire des Écoulements Géophysiques et Industriels) dans les solveurs NSMB (solveur compressible structuré multiblocks parallèle avec maillage chimère) et OpenFOAM (Open source Field Operation And Manipulation). Les modèles de mélange homogène ou un fluide avec une équation d'état de barotrope effectués au LEGI ont réalisé dans le solveur NSMB. Les modèles proposés ont été validés à l'aide de divers cas de test non invasifs, y compris le problème de mouvement de l'interface, le tube de choc eau-air et le tube d'expansion et l'interaction chocbulle. La possibilité d'obtenir des solutions correctes de ces cas de test a été étudiée. Les résultats obtenus à partir des cas de test indiquent que la mise en oeuvre de ces deux modèles de cavitation ne pouvait malheureusement pas être la panacée et être généralisée pour tous les cas de test. Bien que les validations aient montré la capacité des modèles à simuler le développement de la cavitation, les deux modèles souffrent toujours du problème de l'instabilité numérique. La principale différence entre ces deux modèles est que le modèle à trois équations a l'hypothèse d'un équilibre thermodynamique complet entre les phases; par conséquent, cela pourrait expliquer les écarts existant dans les cas de test ci-dessus. Puisque la mise en oeuvre et la validation dans le solveur NSMB avaient déjà pris trop de temps, afin d'atteindre les objectifs de cette étude, qui sont la turbulence et la cavitation, un autre logiciel open source libre, OpenFOAM, a été adopté pour effectuer les cavitations dans un venturi.

Les modèles à quatre équations qui sont composés de trois lois de conservation pour le vi mélange plus une équation de transport pour le taux de vide dans le solveur OpenFOAM appelée interPhaseChangeFoam est étudié. Une comparaison de divers modèles de cavitation couplés à des modèles de turbulence sur la géométrie Venturi 2D et 3D a été proposée. Le solveur interPhaseChangeFoam a été utilisé pour simuler la poche de cavitation par la formulation de modèles de cavitation à équation de transport à rapport de vide, y compris les modèles Kunz, Merkle et SchnerrSauer. Pour la fermeture de la turbulence, trois modèles sont considérés: le modèle Spalart-Allmaras à une équation, le modèle kε à deux équations et le modèle Menter kω SST. Le limiteur de turbulence Reboud est introduit pour réduire la viscosité turbulente afin de capturer la dynamique du jet ré-entrant. Les résultats numériques ont été comparés à des données expérimentales concernant la ration de vide moyennée dans le temps et la vitesse longitudinale, la pression pariétale, les fluctuations de pression de paroi RMS et la viscosité tourbillonnaire turbulente. Les résultats ont montré que l'utilisation d'un limiteur de turbulence par turbulence permet au modèle de simuler correctement les comportements instables de la feuille, cependant de grandes différences apparaissent entre les modèles et l'effet de la réduction n'est pas assez fort. En général, les trois modèles de cavitation étaient capables de reproduire le phénomène de jet ré-entrant, mais la longueur de la cavité était sur-prédite. Parmi les résultats issus de la simulation qui ont été comparés aux données expérimentales, c'est le modèle de cavitation de Kunz couplé au modèle de turbulence kω SST qui pourrait avoir une meilleure prédiction pour la géométrie Venturi. De plus, l'effet 3D n'a pas beaucoup amélioré la prédiction en fonction des résultats numériques obtenus. Ceci peut être dû au problème d'étalonnage du terme de transfert de masse du taux de condensation et du coefficient de vitesse de vaporisation ou au manque de cohérence thermodynamique. Aussi, l'impact sur la valeur de l'exposant n utilisé dans cette correction doit être étudié. En outre, interPhaseChangeFoam est un solveur incompressible qui est moins capable de résoudre le type de géométrie interne.

vii

INTRODUCTION

Background of cavitation

C

avitation is a phenomenon that occurs frequently in conventional hydraulic components such as pumps, valves, turbines and propellers. Over-speeds imposed by the local geometry, shear phenomena, acceleration or vibration may cause local pressure drops in the fluid. When the flow pressure is less than the vapor pressure of the fluid, there is a partial vaporization and vapor structures arise. The so formed structures are entrained by the flow and when they reach a higher pressure zone they condense and implode violently. Cavitation leads to significant loss of system performance, problems of instability of operation of machines and erosion of the component walls. It is thus a primary source of technical problems in the field of hydraulic turbomachinery, naval propulsion and space as well as in high pressure fuel injection.

However, it should be noticed that in certain cases cavitation has a desired effect, for example, supercavitation for underwater vehicles such as torpedoes. The gaseous cavities enveloping the external body make it possible to reduce the friction drag. In addition, cavitation is used for the purpose of cleaning by the control of erosion.

The mechanisms of the process of cavitation and boiling are similar except that in boiling, the vaporization occurs with only small pressure change. In contrast to boiling, the vaporization in cavitation occurs under only a minor temperature change (Figure 1.1).

In the development of a space launcher, cavitation is one of the most limiting factor generated by the hydraulic because it requires from the design phase the introduction of safety margins resulting primarily from an increase in pressure in the reservoirs. This increase in pressure requires an increase in the wall thickness which generates an increase in the structure. The magnitude of this increase in dry weight is 100 k g for 100 mbar of additional pressure, which corresponds about to 2% of the total weight of the largest telecommunications satellite built.

Cavitation appears in the ergol turbo pumps of the launcher propellant and it generates falls of performances, instability of operation as well as mechanical loads on structures. The consequences can be tragic as the failure of the Japanese H-II launch vehicle in 1999.

As for the shipbuilding industry, cavitation is one of the major constraints in the design of marine propellers. Noise, vibration, erosion as issues resulting of cavitation are very tricky.

The appearance and disappearance of bubbles on the propeller blades create local pressure fluctuations that can be compared to shock waves because of their violence. Moreover propeller produces a rotating flow in its wake. Sections of rudders that are placed behind the propeller are then in incidence and can cavitate violently at high speed. Cavitation is also very energetic and very noisy in the audible range. Depending on the type of cavitation frequencies and very specific signatures appear. This type of nuisance is obviously crucial for military vessels, as brought up to 100 km offshore by poorly controlled cavitation. The determination of cavitation instabilities regime is essential.

In the hydraulic energy field, cavitation is a limiting phenomenon in the design phase of hydraulic machinery (pumps, turbines) and its consequences in terms of erosion of the walls are a very important nuisance (operating range and duration component life). Damage to solid walls (Figure 1.2) is caused by very short pressure spikes (10ns to 1µs), high amplitude (∼ 1GPa), attributed to the impact of pressure waves emitted during the collapse of vapor structures.

Knowledge of the dynamics of pockets is therefore very important. Also operating machinery instabilities related to the hydrodynamic coupling between the inter-blade channels are observed.

Types of cavitation

There exists different patterns of cavitation according to the flow configuration, the properties of the fluid and the geometries. Generally, there are four basic types of cavitation and are described briefly below:

• Traveling cavitation

These bubbles are formed in the zone of low pressure, travel with the flow and implode after when they enter the region of higher pressure. This kind of cavitation is observed particularly in the blades of turbine or propeller (Figure 1.3). • Tip-vortex cavitation At the tips of the rotating blade or wing, the pressure may be very low locally which will generate a filament-looking cavitation (Figure 1.6).

Cavitation inception

It is conventional to distinguish whether the flow is cavitating or not by means of cavitation number, σ ∞ , which is defined as

σ ∞ = P ∞ -P vap 0.5ρ ∞ U 2 ∞ (1.1)
This parameter relates the vapor pressure, P vap , to the free-stream pressure, P ∞ , and the free-stream dynamic pressure.

Once the cavitation number, σ ∞ , is reduced in the flow, cavitation will first be observed to appear at some particular value which can be called the incipient cavitation, σ i .

The pressure coefficient, C P , is given by the relation:

C P = P -P ∞ 0.5ρ ∞ U 2 ∞ (1.2)
Therefore, cavitation number can be compared to the pressure coefficient and the following estimate is considered for cavitation inception, σ i

σ i = -C P min = P min -P ∞ 0.5ρ ∞ U 2 ∞ (1.3)
where C P min is the minimum pressure coefficient.

With these definitions above, it is useful to consider that if P min = P vap or σ ∞ = -C P min , the incipient cavitation occurs which means the limiting regime between the non-cavitating and cavitating flow. If further reduction in cavitation number which implies that σ ∞ < -C P min , the developed cavitation happens with an increase in the size and number of bubbles.

CHAPTER 1. INTRODUCTION

Objectives and organization of this thesis

Cavitation for most engineering applications is turbulent, and the interplay between cavitation and turbulence makes the cavitation dynamics even more complicated, and thus the detail dynamics of the phase change is not well understood. Specific issues to numerical techniques in this type of flow also persist. The objectives of this thesis are to implement several cavitation models in the NSMB solver. The emphasis is placed on the study and implement of the Homogeneous Equilibrium Models (HEM) coupled with a barotropic state law and a void ratio Transport-based Equation Model (TEM). The TEM based method for the void ratio including the source terms for vaporization and condensation in the free, open source software OpenFOAM (Open source Field Operation And Manipulation) is also presented on the Venturi geometry to capture the re-entrant jet phenomenon. For the turbulence closure, a density correction approach proposed by Reboud is imposed to several turbulence models.

Besides the introduction, which presents the background of cavitation and the objectives of the study, the thesis is organized as follows.

In Chapter 2, a literature review for the modeling of two-phase flows is investigated which presents the theory in the modeling of cavitating flow, including the different models used for the present work.

In Chapter 3, the flow solvers, the NSMB and OpenFOAM, used in this study are described, including the essential elements of the governing equations, the modeling concepts and the numerical schemes.

In Chapter 4, different test cases carried out by the NSMB solver are presented together with validations against exact solutions of the Euler equations and the models implemented in the solver.

In Chapter 5, the 2D and 3D Venturi geometry are performed by OpenFOAM with the built-in solver interPhaseChangeFoam coupled with different turbulence models. Validation and comparisons are done with experimental measurements including time-averaged void ratio and velocity profiles, RMS wall pressure fluctuations.

Finally, conclusions and future investigations are discussed in Chapter 6.

C H A P T E R 2

REVIEW OF CAVITATION MODELING

N umerical prediction of cavitation remains a challenge for several reasons. First the modeling of phase transition (thermodynamics) and the interactions with the turbulence is not fully established. In addition, it is a complicated task to deal with the large variations of density between the liquid and vapor phases. Specific issues to numerical techniques in this type of flow also persist. On the issue of numerical architecture (compressible or incompressible low Mach preconditioning extended to variable densities), the question remains open. However, several studies have shown better capture re-entrant jet of cavitation bubbles by compressible codes [START_REF] Venkateswaran | Computation of multiphase mixture flows with compressibility effects[END_REF]Goncalvès et al., 2010a;[START_REF] Park | Pressure-based solver for incompressible and isothermal compressible flows with cavitation[END_REF][START_REF] Skoda | Comparison of compressible explicit density-based and implicit pressure-based cfd methods for the simulation of cavitating flows[END_REF].

Modeling of two-phase flows

In this chapter only the modeling of gas-liquid flows are presented. There exists two main approaches for the gas-liquid flows :

• Direct or interface-based methods

• The averaged or diffusion methods of the interface

Direct resolution methods

The so-called direct resolution methods allow to solve all the spatial and temporal scales of the two-phase flows. These kinds of methods reconstruct the interfaces and describe the propagation of the flow, while solving the Navier-Stokes equations.

There are different ways of representing the spatial and temporal evolution of an interface :

CHAPTER 2. REVIEW OF CAVITATION MODELING

• Front tracking method (Lagrangian)

• Level Set method (Eulerian)

• Volume Of Fluid method (Eulerian)

• Diffuse interface method ( [START_REF] Jamet | Methodes a interfaces diffuses pour la modelisation des ecoulements dihasiques[END_REF])

Because of the existence of various velocities at the interface i.e. liquid phase velocity, vapor phase velocity and interface velocity, phase changes are difficult to be taken into account in these kinds of methods. Moreover, the reconstruction of the interface in three-dimensional flows can be difficult and very time consuming.

The average resolution methods

In most of these problems, it is not necessary and would be extremely difficult to know the instantaneous values of the local variables of the flow due to the limitation of the capabilities of computers and the difficulty in predicting the position of the interfaces. The prediction of "averaged" properties are mostly interested in, such as the pressure drop in a bubble flow, the volume flow rate in a conduit etc... For this purpose, "averaged" forms of the equilibrium equations will be used to predict mean values of the flow parameters which are meaningful and experimentally accessible. Moreover, since the equations of equilibrium appear in the form of partial differential equations, it is desirable that the mean properties and their first derivatives, spatial and temporal, should be continuous.

The presence of interfaces leads to serious difficulties for the mathematical formulation of the problem, in the same way as the shock waves in single phase.

The concept of beginning with these methods is the use of instantaneous conservation laws of fluid mechanics for each phase. The interfaces appear as surfaces of discontinuity for the different properties of the fluid, so the fundamental equilibrium equations are expressed in the form of "averaged interface conditions".

There are many ways to "average". Averaging of conservation laws can be carried out:

• in space

• in time

• statistically from a set of measures

• or by a combination of the preceding ones (space/time, statistics/space...).

MODELING OF TWO-PHASE FLOWS

Spatial averaging has been mainly used in the field of nuclear engineering (average over a section of a pipe). It allowed the development of 1D code for the safety analysis of nuclear reactors by averaging the equations on the section of a pipe.

Similar to the use of the RANS approach for turbulent single-phase flows, the temprol averaging is widely used for two-phase flows, especially if they are turbulent. Indeed, since transport phenomena are highly dependent on local fluctuations of variables, it is easier in this case to link the laws of state and behavior needed to close the problem with experimental measurements [START_REF] Ishii | Thermo-Fluid Dynamics of Two-Phase Flow[END_REF] .

Local time-averaged equations

In single-phase turbulent regime, an approach in the sense of Reynolds averaged which treats the instantaneous Navier-Stokes equations statistically is used. For a steady flow, the overall average of equations (average obtained over a large number of realizations) can be replaced by a temporal averaging (ergodic hypothesis).

In the two-phase flow; the location of the interface is unknown in time and space, the instantaneous equations can not be solved. The equations are averaged by decomposing each variable into an average part and a fluctuating part.

The temporal averaging operator of the instantaneous equations reveals the presence rate α, defined by:

α = T k T (2.1)
which represents the time T k of the presence of the phase k, with respect to a duration T.

After spatial discretization of the computational domain, the presence rate is averaged over each cell and is then expressed as the volume fraction:

α = V k V (2.2)
where V k is the volume of the phase k in a volume mesh V .

The different models

Different classes of models are present in the literature according to the number of conservation laws treated and the assumptions made: equilibrium model/relaxed model, homogeneous model/two-velocity model, two-fluid model/one-fluid model:

• Two-fluid models

The full seven-equation two-phase models proposed by Baer et Nunziato [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF] are the most complete. These models take into account explicitly the nonequilibrium effects between phases (unequilibrium of pressure, velocity and temperature) CHAPTER 2. REVIEW OF CAVITATION MODELING but remain difficult to be used in industrial flows (turbomachinery). A seven-equation model has been used for supercavitation and expansion tube problems by Saurel [START_REF] Métayer | Modelling evaporation fronts with reactive riemann solvers[END_REF][START_REF] Saurel | A multiphase model for compressible flows with interfaces, shocks, detonation waves and cavitation[END_REF]. The two-fluid method remains more suited for inviscid and simple geometries [START_REF] Métayer | Modelling evaporation fronts with reactive riemann solvers[END_REF]Saurel et al., 2008a;[START_REF] Petitpas | Diffuse interface model for high speed cavitating underwater systems[END_REF][START_REF] Zein | Modeling phase transition for compressible two-phase flows applied to metastable liquids[END_REF][START_REF] Saurel | A multiphase model for compressible flows with interfaces, shocks, detonation waves and cavitation[END_REF]Yeom andChang, 2006, 2013].

• One-fluid homogeneous mixture models

The models are composed of three conservation laws written for the mixture and are based on a assumption of non-slip between the phases. With the assumption of thermodynamic equilibrium, the Homogeneous Equilibrium Models (HEM) are constituted. The non-equilibrium effects can be introduced empirically [START_REF] Yoon | BIBLIOGRAPHY Choking flow modeling with mechanical and thermal non-equilibrium[END_REF]. Different equations of state for the mixture have been developed in cavitation in a thermosensitive fluid :

barotropic law [START_REF] Cooper | Analysis of single and two-phase flow in turbopump inducers[END_REF][START_REF] Rapposelli | A barotropic cavitation model with thermodynamic effects[END_REF], algorithm for calculating temperature based on the equality of the free enthalpies between the phases [START_REF] Edwards | Low-diffusion flux splitting methods for real fluid flows with phase transition[END_REF]].

• Reduced models with five equations These models are obtained from a simplification of the complete two-fluid model. The archetype five-equation model is the one of Kapila [START_REF] Kapila | Two-phase modeling of deflagration-to-detonation transition in granular materials: Reduced equations[END_REF] which is composed of two conservation equations for masses, one conservation equation for the mixture momentum, one conservation equation for the mixture energy and one non-conservative equation for the void ration to describe the flow topology. They involve two temperature which makes it possible to reproduce thermal non-equilibrium effects, as proposed in the model of Saurel [START_REF] Saurel | Modelling phase transition in metastable liquids: application to cavitating and flashing flows[END_REF] for cavitation simulation in diesel injectors. Some formulations have been proposed to the simulation of interface between two fluids [START_REF] Allaire | A five-equation model for the simulation of interfaces between compressible fluids[END_REF][START_REF] Kreeft | A new formulation of kapila's five-equation model for compressible two-fluid flow, and its numerical treatment[END_REF][START_REF] Murrone | A five equation reduced model for compressible two phase flow problems[END_REF][START_REF] Tian | [END_REF].

• Relaxed models with four equations

A four-equation model was developed for a flashing flows and ebullition applications :

the Homogeneous Relaxation Model (HRM). It consists of three conservation laws for the mixture and one additional transport equation for the void ratio. The latter contains a relaxation source term. The source term involves a relation time that is the time for the system to regain its thermodynamic equilibrium state. This relaxation time is very difficult to determine and is estimated from experimental data [START_REF] Barret | Schemes to compute unsteady flashing flows[END_REF][START_REF] Downar-Zapolski | The non-equilibrium relaxation model for one-dimensional flashing liquid flow[END_REF]. Another formulation of the relaxation term was proposed by Helluy [START_REF] Helluy | Relaxation models of phase transition flows[END_REF], based on a constrained convex optimization problem on the mixture entropy.

Another four-equation model which is very popular to simulate cavitating flows in cold water has been adapted to cryogenic application [START_REF] Hosangadi | Numerical study of cavitation in cryogenic fluids[END_REF][START_REF] Utturkar | Recent progress in modelling of cryogenic cavitation for liquid rocket propulsion[END_REF][START_REF] Zhang | Computational fluid dynamic study on cavitation in liquid nitrogen[END_REF] by adding a transport equation for the void ratio : the

MODELING OF TWO-PHASE FLOWS

Transport-based Equation Model (TEM). This equation includs a cavitation source term for the modeling of condensation and vaporization. The main difficulty is related to the formulation of the source term and the tunable parameters involved for the vaporization and condensation process. The calculation of the void fraction by an additional transport equation including the source terms for vaporization and condensation processes is increasingly used for this model. In this case, the term of mass transfer between phases must be treated explicitly. Several empirical formulations have been proposed to simulate cavitating flows [START_REF] Ahuja | Simulations of cavitating flows using hybrid unstructured meshes[END_REF][START_REF] Wang | Large eddy simulation of a sheet/cloud cavitation on a {NACA0015} hydrofoil[END_REF][START_REF] Merkle | Computational modeling of the dynamics of sheet cavitation[END_REF][START_REF] Singhal | Mathematical basis and validation of the full cavitation model[END_REF][START_REF] Venkateswaran | Computation of multiphase mixture flows with compressibility effects[END_REF][START_REF] Vortmann | Thermodynamic modeling and simulation of cavitating nozzle flow[END_REF][START_REF] Wu | Time-dependent turbulent cavitating flow computations with interfacial transport and filterbased models[END_REF][START_REF] Morgut | Comparison of mass transfer models for the numerical prediction of sheet cavitation around a hydrofoil[END_REF][START_REF] Kunz | A preconditioned navier-stokes method for two-phase flows with application to cavitation prediction[END_REF][START_REF] Senocak | A pressure-based method for turbulent cavitating flow computations[END_REF][START_REF] Hosangadi | Numerical study of cavitation in cryogenic fluids[END_REF] but still suffer from a calibration problem and thermodynamics inconsistency [START_REF] Goncalvès | Constraints on equation of state for cavitating flows with thermodynamic effects[END_REF]. Different sets of parameters are presented in [START_REF] Utturkar | Recent progress in modelling of cryogenic cavitation for liquid rocket propulsion[END_REF][START_REF] Frikha | Influence of the cavitation model on the simulation of cloud cavitation on 2d foil section[END_REF][START_REF] Agnieszka | Review of numerical models of cavitating flows with the use of the homogeneous approach[END_REF].

The different classes of models are summarized in

The two-fluid model

This model is about the Navier-Stokes equations for those phases involved. Here the case of two phases is considered, where k is the phase index, k=1, 2. This gives the following six conservation equations :

∂α k ρ k ∂t + ∇. α k ρ k u k = Γ k (2.3) ∂α k ρ k u k ∂t + ∇. α k ρ k u k ⊗ u k = -∇(α k p k ) + ∇.(α k τ k ) + α k ρ k F k + M k (2.4) ∂α k ρ k E k ∂t + ∇. α k ρ k E k u k = -∇. α k q k -∇. pu k + ∇. τ k .u k + α k ρ k F k .u k + Q k (2.5)
CHAPTER 2. REVIEW OF CAVITATION MODELING E = e + 1 2 u 2 is the specific total energy. Γ k , M k , Q k are the source terms relating to transfers of mass, momentum and energy between phases. They represent the interfacial effects and must be modeled.

M k = M Γ k + P kI ∇α k + F d k (2.6)
The term M Γ k represents the momentum transfer due to the mass transfer. F d k corresponds to the interfical friction force exerted on the phase k. P kI is the pressure of phase k at the interface.

Q k = H Γ k -p kI ∂α k ∂t + F d k .u kI + Q kI (2.7)
H Γ k = L vap Γ k represents the energy transfer due to the mass transfer, where L vap is the latent heat of phase change. Q kI corresponds to the interfacial heat transfer. u kI is the vector of velocity of phase k at the interface.

In addition:

2 k=1 M k = M m = 0 and 2 k=1 Q k = Q m = 0 (2.8)
It should notice that these two terms are not necessary equal to zero although they are generally be taken like that. Indeed due to the variation of the curvature of the interface, the momentum and the energy provided by one phase are not equal to those received by the other.

The one-fluid model

This model, also known as homogeneous mixture approach of two-phase flow consists in writing the averaged Navier-Stokes equations for a "mixing" fluid. It is assumed that the two phases move at the same velocity (i.e. neglecting the drag term between phases). The exchanges and the unequilibrium between phases are then no longer directly modeled, but it is possible to represent them in the closure of the system. Actually, the equation of state of the mixture may introduce a difference at the saturation point (for example, the barotropic law).

A physical property of the mixture is defined by a weighting of the void ratio to its value between each phase. For the weighting of the extensive properties, the density will be used.

ρ m = αρ V + (1 -α) ρ L and ρ m e m = αρ V e V + (1 -α) ρ L e L (2.9)
The conservation equations are as follows :

∂ρ m ∂t + ∇. ρ m u m = 0 (2.10) ∂ρ m u m ∂t + ∇. ρ m u m ⊗ u m = -∇(p m ) + ∇.(τ m ) + ρ m F m (2.11) ∂ρ m E m ∂t + ∇. ρ m E m u m = -∇. q m -∇. pu m + ∇. τ m .u m + ρ m F m .u m (2.12)
It can be observed that the energy required for phase change, the latent heat, does not appear explicitly in the energy conservation equation. In fact, this term is treated implicitely for the mixture.

Four-equation models

These models are intermediate models between one-fluid and two-fluids ones. It consists of solving the conservation equations for the mixture plus a continuity equation for one phase. This makes it possible to treat the mass transfer term explicitly.

∂ρ m ∂t + ∇. ρ m u m = 0 (2.13) ∂ρ m u m ∂t + ∇. ρ m u m ⊗ u m = -∇(p m ) + ∇.(τ m ) + ρ m F m (2.14) ∂ρ m E m ∂t + ∇. ρ m E m u m = -∇. q m -∇. pu m + ∇. τ m .u m + ρ m F m .u m (2.15) ∂α 1 ρ 1 ∂t + ∇. α 1 ρ 1 u 1 = Γ 1 (2.16)
There exists different models according to the modeling of the mass exchange term between the phases. This assumption has the effect of decoupling the mass conservation equation and the momentum conservation equation with the energy conservation equation. In fact, the temperature no longer appears in the first two equations,therefore it has no more influence on the other physical properties.

The equations of state

Tait law

For the case of a slightly compressible flow it is possible to take into account the compressibility of a fluid by the relation :

∆P = c 2 ∆ρ
Tait law : ρ ρ re f = [ n] P + P 0 P re f + P 0 where ρ re f and P re f are reference density and pressure. For water, P 0 = 3 × 10 8 and n = 7.

It is the formulation used by [START_REF] Venkateswaran | Computation of multiphase mixture flows with compressibility effects[END_REF][START_REF] Pouffary | Simulation numérique d'écoulements 2D/3D cavitants, stationnaires et instationnaires[END_REF] to take into account the compressibility in the pure phases for the modeling of cavitation. The speed of sound c is a given value for each phases.

Perfect gas law

This state law allows to model a large number of gases with a good approximation: PV = nRT avec R=8.314 J/(K.kg).

It is also written in the form: P = ρrT where r = R/M = C p -C v (=287 SI unit for air).

According to the internal energy : P ρ, e = (γ -1)ρe where γ = C p C v is the ratio of specific heats. With Joule's law : ∆e = C v ∆T and ∆h = C p ∆T where C v and C p are constants.

There is also the semi-perfect gas law, which defines C p (T) and C v (T) no longer to be constant, but by using polynomial laws as a function of temperature.

Van der Waals law

This law was first introduced by van der Waals in 1873. It contains two constants a and b which are calibrated on the behavior of the fluid at the critical point. It represents one of the first state laws for real gases.

P + a v 2 (v -b) = rT where v is the specific volume (2.17)
This law produces a negative sound speed (dP/dρ < 0) in the phase transition zone (unstable thermodynamic equilibrium).

Stiffened gas law

This low is detailed in [START_REF] Rolland | Modélisation et résolution de la propagation de fronts perméables[END_REF]. It is valid for a large number of fluids, and is sometimes used for solids : P ρ, e = (γ -1)ρ(eq) -γp ∞

The term (γ-1)ρ(eq) represents the intermolecular repulsive effect. The term -γp ∞ represents the molecular attraction which is responsible for the cohesion of liquids or solids. This term is null for the perfect gas state law.

It is set for each fluid by the constants γ and p ∞ (q=0). In the phase change the parameter q, which refers to the energy of the fluid at a given reference state, is non-zero.

The heat capacities are constants in the approximation of stiffened gas law. In the same way as for the perfect gas law, a semi-stiffened gas law makes it possible to define C v and C p by polynomial laws as a function of temperature.

Several sets of parameters for cold water have been proposed as shown in Table 2.2 : [START_REF] Saurel | A multiphase godunov method for compressible multifluid and multiphase flows[END_REF] 4.4 6 × 10 8 0 -1625 Barberon et Helluy [START_REF] Barberon | Finite volume simulation of cavitating flows[END_REF] 3 8.533 × 10 8 -0.1148 × 10 7 4200 1569 Paillere et al. [START_REF] Paillere | On the extension of the ausm+ scheme to compressible two-fluid models[END_REF] 2.8 8.5 × 10 8 0 4186 1486 Le Metayer et al. [START_REF] Metayer | Elaborating equations of state of a liquid and its vapor for two-phase flow models[END_REF] 2.35 10 9 -0.1167 × 10 7 4268 1300 Chang et Liou [START_REF] Chang | A robust and accurate approach to computing compressible multiphase flow: stratified flow model and AUSM+-up scheme[END_REF] 1.932 1.1645 × 10 9 0 8095 1487 Table 2.2: Parameters of the stiffened gas law for cold water by different authors

Authors γ P ∞ (Pa) q (J/kg) C p (J/K.kg) c (m/s) Saurel et Abgrall

Tamman law

This law is equivalent to the stiffened gas law :

P + P c = ρ L K(T + T c )
The use of parameters P c , K, T c , is another formulation but is equivalent to those of stiffened gas law q, P ∞ and γ.

Mie-Grüneisen type law

This law is written as : P(ρ, e) = P ∞ (ρ) + Γ(ρ)ρ ee re f (ρ)

where

Γ = 1 ρ ∂p ∂e ρ
is the coefficient of Grüneisen and P ∞ (ρ) is given as a function of the fluid.

The stiffened gas law is obtained with the assumption of low density variations from the Mie-Grüneisen law. For isentropic evolutions, it becomes the Tait law. Another particular case : if P ∞ is null, then the perfect gas law is obtained.

Benedict-Webb-Rubin law

To get as close as possible to the representation of real gases, there are even more complex form of state laws such as the Redlich-Kwong-Soave equation or the Benedict-Webb-Rubin equation [START_REF] Benedict | An empirical equation for thermodynamic properties of light hydrocarbons and their mixtures: methane, ethane, propane and n-butane[END_REF].

The Benedict-Webb-Rubin law is written as :

P = RT d + d 2 RT (B + bd) -A + ad -aαd 4 - 1 T 2 C -cd 1 + γd 2 exp -γd 2
With P the pressure, R the perfect gas constant, T the temperature, d the molar density, and a, b, c, A, B, C, α, γ the empirical parameters. This law is for example used to represent refrigerants. It is used to characterize hydrogen in the formulation "condensable fluid" in the code Fine T M /Turbo.

Presentation of different models of cavitation

In this section, a review of various models available in the literature that describe the phenomena of cavitation with or without the consideration of thermodynamic effect is presented.

In cold water, or more generally for a non-thermosensitive fluid, the dynamic and thermal phenomena are decoupled. The energy equation is therefore not necessary.

In contrary, in thermosensitive fluid, it is necessary to include the equation of energy.

Models with the mixture state law

These are models with three equations (or two equations without the energy) for which the phase change is controlled by a state law. There are several types of closure relations to link the two phases in the literature :

• Sinusoidal barotropic law [START_REF] Delannoy | Two phase flow approach in unsteady cavitation modelling[END_REF] • Logarithmic barotropic law [START_REF] Schmidt | A fully compressible, two-dimensional model of small, high-speed, cavitating nozzles[END_REF][START_REF] Moreau | A numerical study of cavitation influence on diesel jet atomisation[END_REF][START_REF] Xie | Isentropic one-fluid modelling of unsteady cavitating flow[END_REF] • Saurel's equilibrium law [Saurel et al., 1999] • Tabulated state law [START_REF] Ventikos | A numerical study of the steady and unsteady cavitation phenomenon around hydrofoils[END_REF][START_REF] Clerc | Numerical simulation of the homogeneous equilibrium model for two-phase flows[END_REF] • Equilibrium law based on free enthalpy [START_REF] Edwards | Low-diffusion flux splitting methods for real fluid flows with phase transition[END_REF] • Polynomial law (of degree 5) [START_REF] Song | Current status of cfd for cavitating flows[END_REF] • Barotropic law "Italian" [START_REF] Rapposelli | A barotropic cavitation model with thermodynamic effects[END_REF][START_REF] Sinibaldi | A numerical method for 3D barotropic flows in turbomachinery[END_REF].

• State law based on entropy [START_REF] Barberon | Finite volume simulation of cavitating flows[END_REF] • Mixture of stiffened gas law [START_REF] Goncalves | Numerical simulation of cavitating flows with homogeneous models[END_REF] a/ Sinusoidal barotropic law

The barotropic model existing in Fine T M /Turbo was developed by the successive theses of Coutier [START_REF] Coutier-Delgosha | Modelisation des écoulements cavitants: etude des comportements instationnaires et application tridimensionnelle aux turbomachines[END_REF]] and Pouffary [START_REF] Pouffary | Simulation numérique d'écoulements 2D/3D cavitants, stationnaires et instationnaires[END_REF]. It was originally proposed by Delannoy et

Kueny [START_REF] Delannoy | Two phase flow approach in unsteady cavitation modelling[END_REF]. This law relates the pressure to the density by a sinusoidal relation :

ρ = ρ L + ρ V 2 + ρ L -ρ V 2 sin p -p vap c 2 min 2 ρ L -ρ V (2.18)
c min represents the minimum speed of sound in the mixture. This law introduces a small nonequilibrium effect on the pressure. The unequilibrium is controlled by the value of c min .

b/ Schmidt's barotropic law

From the integration of the Wallis mixture speed of sound which is the propagation velocity of acoustic waves without mass transfer, Schmidt [START_REF] Schmidt | Cavitation in Diesel fuel injector nozzles[END_REF] proposes a barotropic law in the form of :

P = p sat + ρ V c 2 V ρ L c 2 L ρ V -ρ L ρ 2 V c 2 V -ρ 2 L c 2 L ln ρ V c 2 V ρ L + α ρ V -ρ L ρ L ρ V c 2 V -α ρ V c 2 V -ρ L c 2 L (2.19)
This expression is used in [START_REF] Moreau | A numerical study of cavitation influence on diesel jet atomisation[END_REF][START_REF] Dumont | Modélisation de l'écoulement diphasique dans les injecteurs Diesels[END_REF] to simulate the cavitation of diesel in the injectors of piston engine. A modified version was proposed by [START_REF] Xie | Isentropic one-fluid modelling of unsteady cavitating flow[END_REF] in order to avoid the appearance of negative pressure.

c/ Saurel's equilibrium law

For compressible flows, Saurel [Saurel et al., 1999] uses the Tait law for the liquid and the perfect gas law for the vapor to calculate the pressure in each phase. The mixture is assumed to be in kinematic and thermodynamic equilibrium. In this way, there is a logarithmic relation to connect P and T in the form of :

ln(P/P 0 ) = k a k (T/T 0 ) k (2.20)
The densities of each phase are given by polynomial functions of the temperature. The void ratio is defined as :

α = ρ -ρ Lsat (T) ρ V sat (T) -ρ Lsat (T) (2.21)
d/ Edwards equilibrium law [START_REF] Edwards | Low-diffusion flux splitting methods for real fluid flows with phase transition[END_REF] propose an equilibrium model to simulate twophase octane flows. The pure phases are governed by Sanchez-Lacombe's law. Thermodynamic equilibrium is defined by the equality of free enthalpies (g = h -T s) between phases : g L = g V .

The iterative resolution of this equation makes it possible to determine the vapor pressure P vap (T). The void ratio is calculated by : α

= ρ-ρ Lsat (T) ρ V sat (T)-ρ Lsat (T)

e/ Rapposelli's barotropic law

Using thermal analysis on a bubble, a relation between the speed of sound in the two-phase mixture and the temperature can be obtained [START_REF] Rapposelli | A barotropic cavitation model with thermodynamic effects[END_REF]. It is possible to find a law between the density and the temperature by integrating the speed. This law has been used for the calculation of hydrofoil in non-viscous flow.

The speed of sound is expressed as the relation :

1 ρ c 2 = 1 ρ ∂ρ ∂p ∼ = 1 -α p 1 -ε L p ρ L c 2 L + ε L g * p c p η + α γ V p (2.22) In this expression, γ V = C p V
Cv V and ε L represent the liquid fraction participating in the heat exchanges with the vapor and :

ε L = α 1 -α 1 + δ T R 3 -1 (2.23)
where δ T R is a controlled parameter obtained from calibration of the model from experimental results. The other parameters are as follows :

For cold water : g * = 1.67; η = 0.73; P c = 221.29 10 5 Pa For nitrogen : g * = 1.3; η = 0.69; P c = 3.4 10 6 Pa f/ State law based on entropy [START_REF] Barberon | Finite volume simulation of cavitating flows[END_REF] proposed to calculate the entropy of the mixture to evaluate the pressure and the temperature. The pure phases are both governed by the stiffened gas law. The specific entropy of the mixture is maximal at thermodynamic equilibrium. During the process of maximization the entropy can be determined when equilibrium is reached and then also for the pressure P = T ∂s ∂v , where v is the specific volume.

g/ Mixture of stiffened gas law

With the assumption of thermal and mechanical equilibrium, an expression for the pressure and the temperature can be deduced as follows [START_REF] Goncalves | Numerical simulation of cavitating flows with homogeneous models[END_REF] :

P ρ, e, α = (γ(α) -1)ρ(e -q(α)) -γ(α)P ∞ (α) 1 γ(α) -1 = α γ V -1 + 1 -α γ L -1 and ρ q(α) = αρ V q V + (1 -α)ρ L q L P ∞ (α) = γ(α) -1 γ(α) α γ V P V ∞ γ V -1 + (1 -α) γ L P L ∞ γ L -1 T ρ, h, α = h -q(α) C p(α) with ρC p(α) = αρ V C p V + (1 -α)ρ L C p L
The void ratio is computed with saturation values of densities :

α = ρ-ρ Lsat ρ V sat -ρ Lsat .
An extension version considering thermodynamic effects for thermosensible fluids is proposed in [START_REF] Goncalves | Numerical study of cavitating flows with thermodynamic effect[END_REF] by introducing a linear variation relation of P vap , ρ L and ρ V with the temperature.

However, this law failed to obtain reasonable results for Venturi case of 4 degree.

Models with four equation, transport-based equation models (TEM)

In these models, a conservation equation for one of the phases is added by means of the source term S which models the mass exchange between the phases. There are different formulations for the source term (more or less empirical constants) :

• Merkle's model [START_REF] Merkle | Computational modeling of the dynamics of sheet cavitation[END_REF]]

• Kunz's model [START_REF] Kunz | A preconditioned navier-stokes method for two-phase flows with application to cavitation prediction[END_REF] • Senocak and Shyy model [START_REF] Senocak | A pressure-based method for turbulent cavitating flow computations[END_REF] • Saito's model [START_REF] Saito | Numerical analysis of unsteady vaporous cavitating flow around a hydrofoil[END_REF] • Vortmann's model [START_REF] Vortmann | Thermodynamic modeling and simulation of cavitating nozzle flow[END_REF] • Utturkar's model [START_REF] Utturkar | Recent progress in modelling of cryogenic cavitation for liquid rocket propulsion[END_REF] • Hosangadi and Ahuja model [START_REF] Hosangadi | Numerical study of cavitation in cryogenic fluids[END_REF] • Goncalvès model [START_REF] Goncalvès | Numerical study of expansion tube problems: Toward the simulation of cavitation[END_REF] • Source term based on the simplified Rayleigh-Plesset equation a/ Merkle's model (1998)

The model proposed by [START_REF] Merkle | Computational modeling of the dynamics of sheet cavitation[END_REF]] is one of the first models that uses the mass conservation equation for the vapor phase to simulate the cavitation.

The equation solved for the vapor phase is as follows:

∂x V ∂t + u.∇x V = - x V τ V = x L τ L (2.24)
where x V and x L are the mass fractions of the vapor and liquid phases respectively (αρ

V = x V ρ).
The source term is defined as:

1 τ V = 0 when P < P vap 1 kτ re f P-P vap q when P > P vap τ L is defined in the same way for condensation.

τ re f = L re f
U re f is the reference time scale of the fluid, and k is a constant with the value around 10 -3 . The parameter q is not specified in the article [START_REF] Merkle | Computational modeling of the dynamics of sheet cavitation[END_REF]] but seems to be the reference dynamic pressure q = 0.5ρU 2 re f .

b/ Kunz's model (2000)

Kunz's model [START_REF] Kunz | A preconditioned navier-stokes method for two-phase flows with application to cavitation prediction[END_REF] is based on an empirical source term split into two contributions for the evaporation and condensation process :

∂α L ∂t + ∇.(α L u) = ṁ+ + ṁ- (2.25)
This model is implemented in the IZ code [START_REF] Coutier-Delgosha | Simulation of unsteady cavitation with a two-equation turbulence model including compressibility effects[END_REF][START_REF] Rolland | Modélisation et résolution de la propagation de fronts perméables[END_REF][START_REF] Patella | Numerical model to predict unsteady cavitating flow behavior in inducer blade cascades[END_REF]. The evaporation and condensation source terms are given as following expressions :

ṁ-= C dest ρ V α L M in(0, P -P vap ) ρ L (ρ L U 2 re f /2)t ∞ and ṁ+ = C prod ρ V α 2 L (1 -α L ) ρ L t ∞ (2.26)
where t ∞ is the relaxation time, C dest and C prod are the constants to be calibrated.

The condensation rate is modeled as being proportional to the liquid volume fraction and the amount by which the pressure is below the saturated vapor pressure. For the evaporation rate, a simplified Ginzburg-Landau relationship is used.

c/ Senocak and Shyy model (2001)

Senocak et Shyy [START_REF] Senocak | A pressure-based method for turbulent cavitating flow computations[END_REF] try to eliminate the empirical constants by adopting from Kunz's model. It is carried out by the idea of introducing the normal interfacial velocity.

However there will be a problem of locating the interface arises. This difficulty is overcome by the calculation of the density gradient. In this way, a fictitious interface is obtained because of modeling effort inside it (see Figure 2.1). The mass transfer source terms are as follows : Saito [START_REF] Saito | Numerical analysis of unsteady vaporous cavitating flow around a hydrofoil[END_REF]] uses a mass transfer equation for the vapor phase. The system is closed by the modeling of the source term and a mixture state law. The mixture state law is determined by the weighting of each phase form the Tamman law for the liquid phase and the perfect gas law for the vapor phase respectively. :

ṁ-= ρ V α L M in(0, P -P vap ) ρ V (U V ,n -U I,n ) 2 (ρ L -ρ V )t ∞ and ṁ+ = (1 -α L )Max(0, P -P vap ) (U V ,n -U I,n ) 2 (ρ L -ρ V )t ∞ (2.
1 ρ = 1 ρ L (1 -x) + 1 ρ V x or ρ = P P + P c K (1 -x) P(T + T c ) + rx P + P c T (2.28)
The vapor pressure is given by an empirical formula as a function of the temperature. The mass transfer source term is proportional to the pressure difference, P vap -P, as well as the inverse of the square root of the saturation temperature.

ṁ =          ṁ+ = C e Aα (1 -α) ρ L ρ V P * V ap -P 2πRT S if P < P * vap ṁ-= C c Aα (1 -α) P * vap -P 2πRT S if P ≥ P * vap
where T S is the saturation temperature and A = C a α (1 -α)

A denotes the interfacial area concentration in the vapor-liquid mixture.

The saturation vapor pressure of cold water is given by the empirical formula as :

P * vap = 22.13 × 10 6 exp 1 - 647.31 T 7.21379 + 1.152 × 10 -5 -4.787 × 10 -9 T (T -483.16) 2
The parameters C a , C c and C e are empirical constants.

e/ Vortmann's model (2003)

A rate equation for vapor quality x is formulated by Vortmann [START_REF] Vortmann | Thermodynamic modeling and simulation of cavitating nozzle flow[END_REF] as :

∂x ∂t + u.∇x = (1 -x)K l→v -xK v→l (2.29)
The terms K l→v and K v→l mean the probabilities of phase change from liquid to vapor and from vapor to liquid respectively. These terms integrate the Gibbs free energy and involve the relaxation time set to 10 -4 s kg/m 3 . The vapor pressure is supposed to be constant.

f/ Utturkar's model (2005)

The previous IDM model is adapted by [START_REF] Utturkar | Recent progress in modelling of cryogenic cavitation for liquid rocket propulsion[END_REF] to take the thermodynamic effects into account. The new model is then called Mushy Interfacial Dynamics

Model. The original model without thermodynamic effects uses the averaged interface coditions of a liquid-vapor interface to construct the mass transfer source term. This approach is justified by the authors that the sheet cavitation of the cold water contain a significant void ratio. Starting from the analysis of Hord [START_REF] Hord | Cavitation in liquid cryogens, vol 4, combined correlations for venturi, hydrofoil, ogives and pumps[END_REF] for the composition of cryogenic sheet cavitation which describes the vapor zone as the mixture zone with lower void ratio, a model using the averaged CHAPTER 2. REVIEW OF CAVITATION MODELING interface conditions between the liquid and the mixture is formulated.

The mass transfer source terms are given below :

ṁ-= ρ L α L M in(0, P -P vap ) ρ i (U m,n -U I,n ) 2 (ρ L -ρ V )t ∞ and ṁ+ = ρ L (1 -α L )Max(0, P -P vap ) ρ j (U m,n -U I,n ) 2 (ρ L -ρ V )t ∞ (2.30) if α L ≥ 0.99 ρ i = ρ m and ρ j = ρ m otherwise ρ i = ρ V and ρ j = ρ L
This model is valid for the cavitating flows of cold water. As for the sharp IDM model:

U m,n = u.n with n = ∇α L
|∇α L | For the steady calculation, the normal component of interfacial velocity, U I,n is equal to zero.

Numerical simulations are presented for the analysis of the thermodynamic effects for 2D turbulent liquid nitrogen around a warhead. The pressure profiles at the wall in the sheet show a good qualitative behavior of the model, and the void ratio inside the sheet is significantly decreased in comparison with the calculations for cold water.

g/ Hosangadi and Ahuja model (2005)

Hosangadi et Ahuja [START_REF] Hosangadi | Numerical study of cavitation in cryogenic fluids[END_REF] use the source term based on the one of Merkle [START_REF] Merkle | Computational modeling of the dynamics of sheet cavitation[END_REF]]. The formulation has been implemented within a 3D unstructured code CRUNCH. To our knowledge, this is the first one to calculate the performance drop of a liquid cryogenic inducer (LH2). The transport equation of the vapor phase is shown as below :

∂ρ V α ∂t + ∇.(ρ V αu) = m t (2.31)
m t is the mass transfer source term :

m t = ṁ-αρ V + ṁ+ (1 -α)ρ L with : ṁ-=   0 P < P vap 1 τ V U re f L re f P-P vap 1 2 ρ L U 2 re f P > P vap   and ṁ+ =   0 P > P vap 1 τ L U re f L re f P-P vap 1 2 ρ L U 2 re f P < P vap   τ V
and τ L are the time constants for liquid reconversion and vapor formation respectively which were set at 0.001 s. Fluid thermodynamic properties of each phase are calculated from the NIST Chemistyr WebBook (http://webbook.nist.gov).

h/ Goncalvès model (2013)

Goncalvès [START_REF] Goncalvès | Numerical study of expansion tube problems: Toward the simulation of cavitation[END_REF] presents the first version of transport equation which has a form that includes two quantities not used before: the speed of sound, c, and the Wallis mixture speed of sound [START_REF] Wallis | One-dimensional two-phase flow[END_REF], c wall is . The Wallis speed of sound is expressed as a weighted harmonic mean of speeds of sound of each phase:

1 ρ c 2 wall is = α ρ v c 2 v + 1 -α ρ l c 2 l (2.32) 2.

MODELING OF TWO-PHASE FLOWS

The void ratio equation can be expressed as :

∂α ∂t + u ∂α ∂t =   ρ l c 2 l -ρ v c 2 v ρ l c 2 l 1-α + ρ v c 2 v α   =K ∂u ∂x +   c 2 v α + c 2 l 1-α ρ l c 2 l 1-α + ρ v c 2 v α   =1/ρ I the interfacial density ṁ (2.33)
where ṁ is the mass transfer between phases and ρ I is the interfacial density. The term K involves the speed of sound of pure phases and it reflects the effects of changes in volume of each phase. By assuming that the mass transfer is proportional to the divergence of the velocity, the mass transfer ṁ is expressed as :

ṁ = ρ l ρ v ρ l -ρ v 1 - c 2 c 2 wall is div V (2.34)
Computations are performed by Goncalvès and Charrière [START_REF] Goncalvès | Modelling for isothermal cavitation with a four-equation model[END_REF] for several cases including an underwater explosion with cavitation, bubble collapse by a pressure wave and the 8 degree Venturi geometry.

I/ Source term based on the simplified Rayleigh-Plesset equation model

These models, unlike the bubble tracking approach, do not solve the complete Rayleigh-Plesset equation [START_REF] Brennen | Cavitation and bubble dynamics[END_REF] which is written as :

R d 2 R dt 2 + 3 2 dR dt 2 = 1 ρ P vap + P g -P ∞ (t) - 2σ R - 4µ R dR dt (2.35)
This equation describes the evolution of a spherical bubble in an infinite domain of liquid. R is the radius of the bubble, P vap the saturation pressure, P g the pressure of dissolved gas, the last two terms represent the surface tensions and the viscous effects respectively.

• The CAVKA code (CAVitation KArlsruhe)

This code was developed at the University of Karlsruhe by Sauer [START_REF] Sauer | Unsteady cavitating flow -a new cavitation model based on modified front capturing method and bubble dynamics[END_REF][START_REF] Schnerr | Physical and numerical modeling of unsteady cavitation dynamics[END_REF]. It models the cavitation by a void ration equation where the source term S is obtained from a number of nuclei, a characteristic radius of these nuclei and the radius growth rate according to the simplified Rayleigh-Plesset equation.

S = n 0 4πR 2 1 + n 0 4 3 πR 3 dR dt (2.36)
n 0 represent the bubble number density (10 8 by default), and the radius growth rate is expressed as :

dR dt = ± 2 3 P -P vap ρ 1 (2.37) CHAPTER 2.

REVIEW OF CAVITATION MODELING

The expression for the mass transfer source term is derived from the relation between the void ratio and the bubbles :

α = V V V cell = N B 4 3 πR 3 V V + V L = n 0 4 3 πR 3 1 + n 0 4 3 πR 3 (2.38)
where N B represents the number of bubbles in the computational cell and V cell is the volume of the computational cell. A comparison of the results between the CAVKA code and another CATUM code for the unsteady non-viscous flow on 2D hydrofoils is given in [START_REF] Schnerr | Shock and wave dynamics of compressible liquid flows with special emphasis on unsteady load on hydrofoils and a cavitation in injection nozzles[END_REF].

A model with thermodynamic effects was proposed by the same authors [START_REF] Sauer | Unsteady cavitating flow -a new cavitation model based on modified front capturing method and bubble dynamics[END_REF]. They propose to associate the thermodynamic properties of each phase, especially for the vapor pressure, with the temperature obtained from the resolution of the mixture energy equation.. Simulation of unsteady cavitating flow of hot water has been carried out in a 2D nozzle.

• The codes of Fluent and ACE+

The model of Singhal [START_REF] Singhal | Mathematical basis and validation of the full cavitation model[END_REF][START_REF] Dular | Experimental evaluation of numerical simulation of cavitating flow around hydrofoil[END_REF][START_REF] Zhang | Computational fluid dynamic study on cavitation in liquid nitrogen[END_REF][START_REF] Tani | Development and validation of new cryogenic cavitation model for rocket turbopump inducer[END_REF] where the phase change is managed by the conservation equation of vapor phase :

∂x V ρ m ∂t + ∇.(x V ρ m u) = R e -R c
(2.39)

x represents the vapor, R e and R c are the evaporation and condensation source terms respectively :

R e = Ce k σ ρ L ρ V 2 3 P vap -P ρ L 1 -x V -x g if P < P vap (2.40) R c = Cc k σ ρ L ρ L 2 3 P vap -P ρ L x g otherwise (2.41)
Ce and Cc are the empirical constants, k is the local kinetic energy, σ is the surface tension, and x g is the mass fraction of dissolved gases. According to Singhal [START_REF] Singhal | Mathematical basis and validation of the full cavitation model[END_REF],

Ce=0.02 and Cc=0.01.

An extened version taking into account thermodynamic effects is proposed by [START_REF] Tani | Development and validation of new cryogenic cavitation model for rocket turbopump inducer[END_REF]. The formulation shows the vapor pressure P vap (T) varies with the temperature which is derived from an analytical calculation based on the B-factor theory.

• The CFX code

The model contains the TEM model [START_REF] Bouziad | Physical modelling and simulation of leading edge cavitation, application to an industrial inducer[END_REF][START_REF] Mejri | Comparison of computational results obtained from a homogeneous cavitation model with experimental investigations of three inducers[END_REF]. The form of the source term S is based on the simplified Rayleigh-Plesset equation :

S =        S V = F V N V ρ V 4πR 2 B 2 3 |Pvap-P| ρ L if P < P vap S L = F C N C ρ V 4πR 2 B 2 3 P vap -P ρ L if P > P vap
(2.42)

MODELING OF TWO-PHASE FLOWS

N V and N C represent the number of bubbles for different modeling of physical phenomenon (vaporization or condensation).

N V = (1 -α) 3α d 4πR 3 B and N C = 3α 4πR 3 B (2.43)
F V and F C are the empirical constants and represent different time scales for the vaporization and condensation processes : F V =50 and F C =0.01.

In addition, R B is the initial radius for the bubbles and x d is non-condensable gases which provide nucleation sites for the cavitation process : R B = 10 -6 m and x d = 10 -5 by default.

• The Star-CD code

This code contains the Rayleigh-Plesset model. However the expressions of the sources terms are unknow, they are not explained in the articles [START_REF] Kimura | Numerical simulation for unsteady cavitating flow in a turbopump inducer[END_REF][START_REF] Ugajin | Numerical simulation of cavitating flow in inducers[END_REF].

Remark:

The free, open source software OpenFOAM is a one-fluid RANS solver which is developed under the pressure-based schemes (SIMPLE and PISO). Several TEM models are implemented in OpenFOAM [START_REF] Erney | Verification and validation of single phase and cavitating flows using an open source CFD tool[END_REF], including the Kunz's model, Merkle's model and Sauer and Schnerr model. The validation was carried out in a flow behind a hemisphere and two different hydrofoils (NACA0012 and NACA66).

Saurel's model with five equations

This model is composed of four conservation equations (two mass balance for each pure phase, one mixture momentum , one mixture energy) plus a non-conservative equation for the void ratio.

The inviscid formulations are written as :

∂αρ V ∂t + div αρ V u = ṁ (2.44) ∂(1 -α)ρ L ∂t + div (1 -α)ρ L u = -ṁ (2.45) ∂ρu ∂t + div ρu ⊗ u + grad P = 0 (2.46) ∂ρE ∂t + div u(ρE + P) = 0 (2.47) ∂α ∂t + u.∇α = α(1 -α)(ρ L c 2 L -ρ V c 2 V ) αρ L c 2 L + (1 -α)ρ V c 2 V div (u) (2.48) + α(1 -α) αρ L c 2 L + (1 -α)ρ V c 2 V Γ V α + Γ L 1 -α Q + ρ V c 2 V α + ρ L c 2 L 1-α c 2 V α + c 2 L 1-α ṁ
The mass transfer term ṁ and the heat transfer term Q between phases are expressed as :

CHAPTER 2. REVIEW OF CAVITATION MODELING ṁ = νρ(g L -g V ) and Q = H(T L -T V )
These terms involve free Gibbs enthalpy and the temperature of pure phases as well as two relaxation coefficients ν and H.

The property

Γ = 1 ρ ∂p ∂e ρ
is the Grüneisen coefficient of pure phases governed by the stiffened gas law. With the assumption of mechanical equilibrium, it is possible to calculate the mixture pressure :

p ρ, e, α, ρ V , ρ L = (ρe -αρ V q V -(1 -α)ρ L q L ) -α γ V P V ∞ γ V -1 + (1 -α) γ L P L ∞ γ L -1 α γ V -1 + 1-α γ L -1 (2.49)
The relaxation coefficients n and H are unknown and very difficult to determine in practice. To remove this uncertainty, the resolution is carried out in two steps :

1. It is assumed that the thermodynamic equilibrium is reached, that is to say the relaxation coefficients are taken as infinite. This allows to determine the equilibrium void ration α eq . 2. Then the whole system is solved. The mass exchange term between phases is evaluated by :

ṁ = νρ(g L -g V ) = ∂αρ V ∂t = (αρ V ) eq -αρ V ∆t
. The same goes for the heat transfer term Q.

The model has been tested on the problems of 1D expansion tube, 2D supercavitating flow around an obstacle and 2D Venturi flow corresponding to a fuel injector (dodecane).

Two-fluid models a/ Grogger and Alajbegovic model (1998)

This two-fluid model has been initially applied to the cavitation in cold water in 2D and 3D Venturi [START_REF] Grogger | Calculation of the cavitating flow in venture geometries using two fluid model[END_REF]. It has been applied more recently to a three-phase (liquid, vapor and air) turbulent cavitating flows simulation for the high-pressure swirl injector geometry (3D) [START_REF] Alajbegovic | Three-phase cavitating flows in high-pressure swirl injectors[END_REF]. This model is based on the resolution of mass and momentum conservation equations for each phase (four equations in two phases and six equations in three phases). It has the advantage of modeling the interactions of momentum between phases.

The system solved in the two phases is shown as below (k = L, V ) :

∂α k ρ k ∂t + ∇. α k ρ k u k = Γ k (2.50) ∂α k ρ k u k ∂t + ∇. α k ρ k u k ⊗ u k = -α k ∇p + ∇.(α k τ k ) + F d k + u k Γ k (2.51)
The mass transfer source term is based on the bubble growth rate from the simplified Rayleigh-Plesset equation :

Γ L = ρ L N4πR 2 ∂R ∂t = -Γ V (2.52)
The bubble radius R is defined in function of the void ratio α according to the following expres-

sion : R = 1 2 6α Nπ 1/3
N represents the bubble number density : -α) for α > 0.5 where N 0 is the initial bubble number density and is set to 10 12 .

N = N 0 for α ≤ 0.5 2 N 0 -1 (1
The bubble growth rate is modeled by :

∂R ∂t = ± 2 3 p -p V ρ L
The interfacial momentum transfer term is based on the effects of drag force of a sphere and turbulent dispersion forces :

F d L = c TD ρ L k L ∇α + c D ρ L |u rel | u rel 2 A = -F d V with u rel = u L -u V (2.53)
The term A represents the surface of a spherical bubble. C TD is the turbulent dispersion coefficient and k is the turbulent kinetic energy. The drag coefficient c D is given by the relations : This model is based on the seven-equation model of [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF],

c D = 24/R e(1 + 0.15R e 0 .687) if R e P <
which uses the transport equation of the volume fraction α 1 to close the two-fluid model with six equations [START_REF] Saurel | A multiphase godunov method for compressible multifluid and multiphase flows[END_REF]. It was originally proposed without considering the mass transfer. In the thesis of Le Métayer [START_REF] Rolland | Modélisation et résolution de la propagation de fronts perméables[END_REF], a phase change term is introduced to the model. It allows to solve the complex problems such as flows with three components of velocity.

The seven equations of Baer-Nunziato without phase change are as follows (k=1,2) :

∂α k ρ k ∂t + ∇. α k ρ k u k = 0 ∂α k ρ k u k ∂t + ∇. α k ρ k u k ⊗ u k = -∇(α k P k ) + ∇.(α k τ k ) + α k ρ k F k + F d k + P I ∇α k ∂α k ρ k E k ∂t + ∇. α k ρ k E k u k = -∇. α k q k -∇. p k u k + ∇. τ k .u k + P I u I .∇α k + µP I P k -P k + α k ρ k F k .u k + F d k .u I + Q kI ∂α 1 ∂t + u I .∇α 1 = -µ (P 1 -P 2 ) CHAPTER 2.

REVIEW OF CAVITATION MODELING

The index I represents the interface. The term µ (P 1 -P 2 ) which represents the production of the volume fraction α 1 , is equal to the pressure difference between phases multiplied by a coefficient µ that controls the velocity at which pressure equilibrium is reached. This term induces µP I (P k -k ) due to the interfacial pressure work to energy conservation equations.

To close the system, the interfacial pressure P I is represented by the most compressible phase pressure and the interfacial velocity u I is represented by the less compressible phase velocity in the original model. Subsequently, Saurel et Abgrall [START_REF] Saurel | A multiphase godunov method for compressible multifluid and multiphase flows[END_REF] estimate the interfacial pressure by the mixture pressure :

P I = 2 k=1 α k P k
After solving an inert Riemann problem i.e. without phase change, the study carried out by Le

Métayer consists of considering the mass transfer term by the Rankine-Hugoniot relations across the shock front. Different models are proposed. The most complex one deals with three velocities, one for each phase plus one for the front. More details about this model are available in the references [START_REF] Saurel | A multiphase godunov method for compressible multifluid and multiphase flows[END_REF][START_REF] Saurel | A multiphase model for compressible flows with interfaces, shocks, detonation waves and cavitation[END_REF][START_REF] Rolland | Modélisation et résolution de la propagation de fronts perméables[END_REF][START_REF] Bibliography Abgrall | Efficient numerical approximation of compressible multi-material flow for unstructured meshes[END_REF]].

These models have been used and validated on a 1D inviscid flow expansion tube problem.

c/Saurel and Le Métayer model with ten equations (2001)

Saurel et Le Métayer [START_REF] Saurel | A multiphase model for compressible flows with interfaces, shocks, detonation waves and cavitation[END_REF] propose a multiphase model composed of five equations for each phase. It is able to deal with a wide range of applications with the very general formulation (interfaces between compressible materials, homogeneous two-phase flows, the problems of shocks and cavitation). This approach is based on the one proposed by Baer et

Nunziato [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF]. The system of equation for each phase k is written as :

∂α k ρ k ∂t + ∇. α k ρ k u k = m k ∂α k ρ k u k ∂t + ∇. α k ρ k u k ⊗ u k = -∇(α k P k ) + m k u I + F d k + P I ∇α k ∂α k ρ k E k ∂t + ∇. α k ρ k E k u k = -∇. [P k u k ] + P I u I .∇α k + m k E kI -µP I (P k -P k ) + F d k .u I + Q kI ∂α k ∂t + u I .∇α k = µ (P k -P k ) + m k ρ X ∂N k ∂t + ∇. (N k u k ) = Ṅk
with the average interface conditions :

k m k = 0 k m k u I + F d k + P I ∇α k = 0 k P I u I .∇α k + m k E kI -µP I (P k -P k ) + F d k .u I + Q kI = 0
The conservation law of mass, momentum and energy, appear the terms related volume fraction transport equation α k of phase k, as well as the number density of the individual entity N k composing phase k (ex : number of bubbles for a liquid-gas flow mainly liquid). The source term Ṅk of the equation represents the phenomena of breakup or coalescence.

The interfacial velocity u I as well as the interfacial pressure p I of the system are defined as :

u I = k α k ρ k u k k α k ρ k and p I = k α k P k + ρ k (u I -u k ) 2 (2.55)
The momentum transfer term due to the drag force between phases is modeled by the velocity relaxation term :

F d k = λ k u k -u k
The mass transfer term is obtained from the interface-averaged equations :

m 1 = -Q 1I + Q 2I E 1I -E 2I (2.56)
The energy transfer at the interface is provided by empirical correlations :

Q kI = h k T kI -T k A ex (2.57)
where A ex is the exchange interfacial area.

For spherical entities (bubbles, drops) :

A ech = N 1 4πR 2
The heat exchange coefficient is defined as :

h k = λ k N u

2R

The Nusselt number is defined from the Reynolds number and the Prandtl number :

N u = 2 + 0.6R e 0.5 P r 0.33 (2.58)

The µ (P k -P k ) and µP I (P k -P k ) terms are related to the pressure relaxation process and are controlled by the value of µ.

When mass transfer occurs an extra source term, m k ρ X , is present in the transport equation of the volume fraction, where ρ X represents the density of the less compressible phase. The main purpose of this term is to separate the mass transfer and the acoustic propagation during numerical resolution.

d/ Model of Saturne code of EDF

Mimouni et al. [START_REF] Mimouni | Modelisation et simulation des ecoulements cavitants par une approche diphasique[END_REF] present the simulation of cavitation carried out with the NEPTUNE code which is originally developed by EDF and CEA for the study of two-phase flow with ebullition. This code uses a two-fluid approach and allows to simulate a large number of flow configurations : gaseous or liquid phase with solid particles, liquid phase and vapor , as well as different flow topologies : disperesd phases, flows with a continuous phase... The mass transfer term is obtained from the interface-averaged equations :

Γ L = -Γ V = q LI + q V I h V I -h LI A I (2.59)
where q kI represents the interfacial heat flux in the phase k, h kI the enthalpy of the phase k at the interface and A I = 6α d the exchange interfacial area (α is the void ration and d the averaged bubble diameter taken equal to 0.1mm).

The enthalpies of each phase at the interface are assumed to be saturated. The interfacial heat flux is modeled as follows :

q kI = c kI (T sat (P) -T k ) (2.60)
With the specific heats as follows :

c LI = N u L λ L d and c V I = αρ V C p V ∆t
The Nusselt number of liquid is modeled as : The code has also been tested on a case of cavitation behind an orifice (EPOCA test case).

N u L = 2 + 0.6R e 1/2 P r 1/3 L with R

Summary

The vast majority of computer codes dedicated to the simulation of cavitation is based on an averaged approach for both the two-phase flow and the turbulence. With proper averaging, the mean values of fluid motions and properties can be obtained. Within the averaged model family, there are different approaches according to the physic assumptions made on thermodynamics equilibrium and slip condition between phases. This has resulted in the development of various systems ranging from seven (two-fluid) to three (one-fluid) equations only. The two-fluid approach is the most complete and is also known to be a real challenge for numerical simulation due to the complicated characteristics of the equation system and the troublesome non-conservative terms.

SUMMARY

On the other hand, the one-fluid method, or homogeneous method considers the flows as a mixture of two fluids behaving as one that is similar to the single-phase flow. In this way, only one set of conservation equation is employed to express the fluid interaction for the mixture. Because of the difficulty of modeling nonequilibrium thermodynamics pattern during a phase transition, the existing models have systematic use of mechanical equilibrium assumptions (single pressure model) and thermal (single temperature model). Besides, vaporization and condensation processes are assumed to be instantaneous. Then, this method cannot reproduce strong thermodynamic or kinetic non-equilibrium effects. Considering its simplicity and low computation cost, the homogeneous method is more attractive for numerical simulations of cavitating flows. On the assumptions of velocity equilibrium and pressure equilibrium for each phase from the full models, a reduced model five-equation model is obtained. This model is capable of modeling the mass and energy transfer terms between phases and taking the thermo non-equilibrium into account. By assuming the the velocity, pressure and thermal equilibrium between phases, a four-equation model can be expressed. With an additional transport equation, usually the void ratio, the mass transfer between phases can be modeled. The main problem of this model is the formulation of the source term and the tunable parameters involved for the cavitating processes. With the assumption of complete thermodynamic equilibrium between phases, that is local temperature, pressure and free Gibbs enthalpy are supposed to be in equilibrium, the three-equation models or Homogeneous Equilibrium Models (HEM) are derived. The most difficult part for this approach is to define a proper equation of state (EOS) for the thermodynamic behavior of the mixture to close the system.

The one-fluid method (homogeneous method) has received more attention up to now, because of its lower computational cost and easier coupling with turbulence models. Among the various existing models, the main differences count on the relation between the pressure and density fields. This coupling is generally treated through a barotropic equation of state, or developed by the framwork of transport-based equation method. In this study, the HEM based formulations coupled with a barotropic state law [START_REF] Goncalves | Numerical simulation of cavitating flows with homogeneous models[END_REF]] and a four-equation based model completed with a void ratio transport equation [START_REF] Goncalvès | Numerical study of expansion tube problems: Toward the simulation of cavitation[END_REF] Immersed Boundary Method (IBM) have been successfully implemented in NSMB. The details about the solver can be referred to the NSMB Handbook [Vos et al., 2013]. The cavitation models and numerical method that are implemented in this study will be presented in the following sections.

Governing Equations

The governing equations of the HEM based formulations coupled with a barotropic state law and a four-equation based model completed with a void ratio transport equation are presented below.

A three-equation model

The movement of fluids is governed by three basic physical conservation equations: mass, momentum and energy conservation equations. The homogeneous mixture approach is used to model two-phase flows. The phases are assumed to be sufficiently well mixed and the disperse particle size are sufficiently small thereby eliminating any significant relative motion. The phases are strongly coupled and moving at the same velocity. In addition, the phases are assumed to be in kinematic and thermodynamic equilibrium: they share the same pressure, temperature and velocity . The evolution of the two-phase flow can be described by the conservation laws that employ the representative flow properties as unknowns just as in a single-phase problem. The void fraction α is introduced to characterize the volume of vapor in each cell: α = 1 means that the cell is completely filled with vapor; inversely, a complete liquid cell is represented by α = 0.

The density ρ, the center of mass velocity u and the internal energy e for the mixture are defined 3.1. NSMB by [START_REF] Ishii | Thermo-Fluid Dynamics of Two-Phase Flow[END_REF]:

ρ = αρ v + (1 -α)ρ l (3.1) ρu = αρ v u v + (1 -α)ρ l u l (3.2) ρ e = αρ v e v + (1 -α)ρ l e l (3.3)
where the subscripts v and l are the vapor and liquid phase respectively.

The inviscid compressible Navier-Stokes equations in 3-D Cartesian corrdinates (x, y, z) can be expressed in conservative form as:

∂ ∂t (W) + ∂ ∂x ( f ) + ∂ ∂y (g) + ∂ ∂z (h) = 0 (3.4)
where t denotes the time.

The state vector W is given by:

W =          ρ ρu ρv ρw ρE          (3.5)
and the convective fluxes are defined as:

f =          ρu ρu 2 + P ρuv ρuw u(ρE + P)          , g =          ρv ρvu ρv 2 + P ρvw v(ρE + P)          , h =          ρw ρwu ρwv ρw 2 + P w(ρE + P)          (3.6)
Here u, v and w are the Cartesian velocity components, P is the pressure and E is the total energy.

The specific total energy E can be expressed in terms of the specific internal energy e and kinetic energies as:

E = e + 1 2 u 2 + v 2 + w 2 (3.7)
It is sometimes useful to recast the energy equation in terms of enthalpy. The specific total enthalpy is given by:

H = γ γ -1 P ρ + 1 2 u 2 + v 2 + w 2 (3.8)
Also, the specific total energy can be written as:

E = 1 γ -1 P ρ + 1 2 u 2 + v 2 + w 2 (3.9)
From the Equations of 3.8 and 3.9, E = H -P ρ is obtained. To close the system of equations the pressure P must be related to the state vector W. This relation depends on the model used to describe the thermodynamic properties of the gas. The difficulty with this homogeneous approach is to specify an equation of state (EOS) that covers all possible fluid states: pure phases (incompressible region) and two-phase mixture (compressible region).

For the pure phases, the convex stiffened gas EOS [START_REF] Metayer | Elaborating equations of state of a liquid and its vapor for two-phase flow models[END_REF]] is used:

P(ρ, e) = (γ -1)ρ(e -q) -γP ∞ (3.10) P(ρ, T) = ρ(γ -1)C v T -P ∞ (3.11) T(ρ, h) = h -q C p (3.12)
where γ = C p /C v is the polytropic coefficient, C p and C v are thermal capacities,h the enthalpy, q the energy of the fluid at a given reference state and P ∞ is a constant reference pressure. The speed of sound c is given by:

c 2 = γ P + P ∞ ρ = (γ -1)C p T (3.13)
In terms of computational methods, the application of a compressible formulation to simulate low speed cavitating flows results in poor convergence and erroneous calculations. To achieve this goal, a preconditioned method is necessary. The preconditioning matrix proposed by Turkel [Guillard andViozat, 1999] [Turkel, 1987] is used in this research (see Appendix A).

For the two-phase mixture, a sinusoidal barotropic law [START_REF] Delannoy | Two phase flow approach in unsteady cavitation modelling[END_REF]] is applied:

P(ρ, α) = P vap + ρ sat l -ρ sat v 2 c 2 min Arcsin(A(1 -2α)) (3.14) T(ρ, h) = h l -q l C p l = h v -q v C p v = h -q(α) C p (α) (3.15)
This law is characterized by its maximum slope 1/c 2 min . The quantity c min is an adjustable parameter of the model, which can be interpreted as the minimum speed of sound in the mixture. With this barotropic law, there is no coupling with the temperature and the cavitation phenomenon is assumed to be isothermal. In the original approach, pure phases are considered as incompressible and the speed of sound is infinite in each phase. In order to join compressible pure phases, a constant A, close to 1, is introduced to avoid infinite value of speed of sound. The speed of sound can be computed by:

c 2 = ∂P ∂ρ s = ∂P ∂ρ T = Ac 2 min 1 -A 2 (1 -2α) 2
(3.16)

NSMB

A four-equation model

Consider the five-equation Kapila model [START_REF] Kapila | Two-phase modeling of deflagration-to-detonation transition in granular materials: Reduced equations[END_REF] by assuming the liquid is at its saturation state, a four-equation model is obtained. The model consists of three conservation laws for mixture quantities as in Equation (3.4) and an additional equation for the void fraction α.

The void fraction equation can be expressed as:

∂α ∂t + u ∂α ∂x + v ∂α ∂ y + w ∂α ∂z = K ∂u ∂x + K ∂v ∂y + K ∂w ∂z (3.17) K =   ρ l c 2 l -ρ v c 2 v ρ l c 2 l 1-α + ρ v c 2 v α   (3.18)
The term K involves the speed of sound of pure phases and it reflects the effects of change in volume of each phase.

To compute the pressure and the temperature, the convex stiffened gas EOS is used for the pure phases as presented above for the HEM model. For the two-phase mixture, an expression for the pressure and the temperature can be deduced from the thermal and mechanical equilibrium assumption [START_REF] Saurel | Modelling phase transition in metastable liquids: application to cavitating and flashing flows[END_REF] on the basis of the stiffened gas EOS. These expressions are available in all possible fluid states along with the function of the void fraction and the mass fraction of gas Y = αρ v /ρ:

P ρ, e, α, Y = γ (α) -1 ρ (e -q (Y )) -γ (α) P ∞ (α) (3.19) 1 γ (α) -1 = α γ v -1 + 1 -α γ l -1 (3.20) q (Y ) = Y q v + (1 -Y ) q l (3.21) P ∞ (α) = γ (α) -1 (α) α γ v γ v -1 P v ∞ + (1 -α) γ l γ l -1 P l ∞ (3.22)
By assuming the thermal equilibrium between phases, the mixture temperature is expressed as:

T ρ, h, Y = h -q (Y ) C p (Y ) with C p (Y ) = Y C p v + (1 -Y ) C p l (3.23)
Without mass transfer, the propagation of acoustic waves follows the Wood or Wallis speed of sound. This speed c wall is is expressed as a weighted harmonic mean of speeds of sound of each phase:

1 ρ c 2 wall is = α ρ v c 2 v + 1 -α ρ l c 2 l (3.24)

Numerics

The numerical simulations are carried out using an explicit time integration and based on a finite-volume discretization. The convective flux through the cell interface is computed with a HLLC scheme [START_REF] Batten | On the choice of wavespeeds for the hllc riemann solver[END_REF][START_REF] Toro | Restoration of the contact surface in the hll-riemann solver[END_REF]. The method considers two averaged intermediate states U * L and U * R separated by the contact wave of speed S M . The numerical flux Φ i,l at cell interface l can be expressed as:

Φ i,l (U L ,U R ) =                  G (U L ) n i,l if S L > 0 G U * L n i,l if S L ≤ 0 ≤ S M G U * R n i,l if S M ≤ 0 ≤ S R G (U R ) n i,l if S R > 0 (3.25)
where S L and S R are referred to the speeds of the smallest and largest waves at the cell interface.

The normal velocity component V n = V.n. The left (K = L) and right (K = R) states of the variables U * K , and corresponding fluxes G U * K , are defined by:

U * K =             ρ * K ρu * K ρv * K ρw * K ρE * K α * K             = 1 S K -S M             ρ K S K -V n K ρu K S K -V n K + (P * -P K ) n i,l ρv K S K -V n K + (P * -P K ) n i,l ρw K S K -V n K + (P * -P K ) n i,l ρE K S K -V n K + P * S M -P K V n K α K S K -V n K             (3.26) G U * K n i,l =             ρ * K S M ρu * K S M + P * n i,l ρv * K S M + P * n i,l ρw * K S M + P * n i,l ρE * K S M + P * S M α * K S M             (3.27)
where the pressure P * is given by:

P * = P L + ρ L V n L -S L V n L -S M = P R + ρ R V n R -S R V n R -S M (3.28)
And the contact-wave speed S M is defined by:

S M = P R -P L + ρ L V n L S L -V n L -ρ R V n R S R -V n R ρ L S L -V n L -ρ R S R -V n R (3.29) 3.2. OPENFOAM
The HLLC solver requires the estimates of wave speeds S L and S R in the Riemann problem.

A direct and simple wave speed estimation is used:

S L = Min V n L -c L , V n R -c R ; S R = Max V n L + c L , V n R + c R (3.30)
The non-conservative term is discretized following the idea of Daude et al. [START_REF] Daude | Numerical experiments using a hllc-type scheme with ale formulation for compressible twophase flows five-equation models with phase transition[END_REF]. The integral term is approximated with the following relation:

C i B (U) divV dS = Bi l∈∂C i ∂C i V.n i,l dl (3.31)
where Bi is some average of B on cell C i . Here, Bi = B (U i ) is used.

the cell interface value u i,l is expressed as:

u i,l (U L ,U R ) =                  V L .n i,l if S L > 0 S L -V n L S L -S M S M if S L ≤ 0 ≤ S M S R -V n R S R -S M S M if S M ≤ 0 ≤ S R V R .n i,l if S R > 0 (3.32)
Details on the numerical method implemented in the NSMB solver can be referred to the textbook from Blazek [START_REF] Blazek | Computational fluid dynamics: Principles and applications[END_REF] and the course from Goncalvès [START_REF] Da Silva | Résolution numérique des équations d'euler 1d[END_REF]

OpenFOAM

Open Unlike many other commercial CFD packages, OpenFOAM does not have a graphical user interface to help the user in preparing a case. In order to run an OpenFoam case, it involves typically the manual preparation of all the required input and mesh files. The minimum required set of files are system, constant and time directories as shown in Table 3.1.

• system folder

The system directory contains run-time control and solver numerics. The decomposePar file describes how the computational domains to be subdivided for multiple processors. The controlDict file contains general simulation settings such as the time step, duration and data saving interval. The fvSchemes file includes the definition of the numerical schemes implemented to discretize the equations and interpolate the solutions. The fvSolution files defines the parameters used to solve the flow equations and the residual tolerance.

• constant folder

The constant directory contains physical properties, turbulence modeling properties, advanced physics and so on. The polyhedral mesh information is placed in the polyMesh directory.

• 0 folder

The 0 directory contains several individual files for every relevant flow quantity including both the initial and boundary conditions. When a new timestep is computed and stored, it is written as a time directory in the case name directory. The time directories contains the solution and derived fields and are created by the solver according to the setting of saving frequency.

Turbulence Closures

Most cavitation phenomena involve turbulence and the turbulence-cavitation interaction is an under-known and documented phenomenon (due in particular to the difficulty of performing experimental measurements in cavitating flows). Compressibility effects on turbulence and the effects of the dispersed phase are also unknown. The numerical accuracy for turbulent cavitation depends on both cavitation and turbulence modeling. Thus, the choise of a turbulence modeling is an important issue for the simulation of cavitation. Direct numerical simulation (DNS) has the highest capability of resolving all turbulent scales. However, it requires a very fine grid resolution and therefore it is still pretty hard to be applied because of the high consuming of computer performance. Although the Large Eddy Simulation (LES) has already been implemented for the turbulent cavitating flows [START_REF] Wang | Large eddy simulation of a sheet/cloud cavitation on a {NACA0015} hydrofoil[END_REF][START_REF] Huang | Large eddy simulation of turbulent vortex-cavitation interactions in transient sheet/cloud cavitating flows[END_REF][START_REF] Gnanaskandan | Large eddy simulation of the transition from sheet to cloud cavitation over a wedge[END_REF], the usual codes are formulated in a Reynolds-averaged Navier-Stokes (RANS) to tensor turbulent closure model by a transport equation kε (Boussinesq hypothesis) considering the balabce between the computational effort and accuracy. This hypothesis suggests that the turbulent shear or Reynolds stresses could be replaced by the product of the mean velocity gradient and a "turbulent or eddy-viscosity", µ t . In this way, the stress tensor and heat flux vector of the set of transport equations contain additional terms due to the Reynolds stresses.

The Reynolds stresses -ρu i u j (i, j = x, y, z) need to be modeled to close the system of equations. The Boussinesq's gradient transport hypothesis for turbulence closure by using the eddy viscosity concept as the following equation:

-ρu i u j = µ t ∂u i ∂x j + ∂u j ∂x i - 2 3 δ i j µ t ∂u k ∂x k + ρk (3.33)
The eddy viscosity µ t is not a fluid property, but a property that depends on the local turbulence structure. The variable k is the turbulent kinetic energy, defined as k = 1 2 u i u i . The second

The Spalart-Allmaras model

For a one equation turbulence model, the state vector contains six unknowns. The last variable is either ν or F = k 2 with corresponding source terms S ν or S F ,

∂ ∂t (W) + ∂ ∂x ( f -f v ) + ∂ ∂y (g -g v ) + ∂ ∂z (h -h v ) = S (3.34) W = ρ, ρu, ρv, ρw, ρE, ν or F (3.35)          f = ρu, ρu 2 + p, ρuv, ρuw, u(ρE + p), u ν or uF , g = ρv, ρvu, ρv 2 + p, ρvw, v(ρE + p), v ν or vF , h = ρw, ρwu, ρwv, ρw 2 + p, w(ρE + p), w ν or wF . (3.36)            f v = 0, τ xx -ρu x u x , τ x y -ρu x u y , τ xz -ρu x u z , (τU) x -q x , D νx or D F x ,
g v = 0, τ yx -ρu y u x , τ yy -ρu y u y , τ yz -ρu y u z , (τU) y -q y , D νy or D F y ,

h v = 0, τ zx -ρu z u x , τ z y -ρu z u y , τ zz -ρu z u z , (τU) z -q z , D νz or D F z .
(3.37)

   S ν = (0, 0, 0, 0, 0, P ν -Φ ν) ,
S F = (0, 0, 0, 0, 0, P F -Φ F ) .

(3.38)

The Reynolds stresses and heat diffusion are calculated using the Boussinesq approximation.

D denotes the different diffusion terms and Φ the destruction terms.

A transport equation for the turbulent viscosity ν is assembled, using empiricism and arguments of dimensional analysis, Galilean invariance, and selective dependence on the molecular viscosity [START_REF] Spalart | A one-equation turbulence model for aerodynamic flows[END_REF]:

D ν Dt convection = c b1 S ν production + 1 σ [∇ • ((ν + ν)∇ ν) + c b2 (∇ ν) 2 ] di f f usion -c w1 f w (r)( ν d ) 2 dissi pation (3.39)
The eddy viscosity is defined as:

µ t = ρ ν f v1 ≡ ρν t (3.40) 3.3. TURBULENCE CLOSURES
To ensure that ν equals K yu τ in the log layer, in the buffer layer and viscous sublayer, the damping function f v1 is defined as:

f v1 = χ 3 χ 3 + c v1 3 (3.41)
as function of the totally local variable λ:

χ ≡ ν ν (3.42)
The function S must be modified to maintain its log-layer behavior ( Ŝ = U τ /(K d)) all the way to the wall:

S = S 1/2 + ν (kd) 2 f v2 (3.43)
which is accomplished with the help of the function f v2 :

f v2 = 1 - χ 1 + χ f v1 (3.44)
The destruction term should vanish in the outer region of the boundary layer. Spalart-Allmaras proposed the function:

f w (r) = g[ 1 + c 6 w3 g 6 + c 6 w3 ] 1/6 (3.45)
with the argument r:

r = ν (kd) 2 S (3.46)
Both r and f w equal 1 in the log layer, and decrease in the outer region.

g = r + c w2 (r 6 -r) (3.47)
the function g is merely a limiter that prevents large values of f w .

The constants of the Spalart-Allmaras model are:

c b1 = 0.1355, c b2 = 0.622, c w2 = 0.3, c v1 = 7.1, σ = 2 3 , c w1 = c b1 k 2 + (1 + c b2 ) σ , c w3 = 2 (3.48)
The kε model Two-equation turbulence models are widely used, as they offer a good compromise between numerical effort and computational accuracy. The kε model requires the solution of transport equations for the turbulent kinetic energy k and the turbulent dissipation rate ε, the same for kω models. This basis model is typically a "high Reynolds number" model. To take into account for the interaction between turbulence and fluid viscosity, many different low-Reynolds number versions have been implemented. These versions differ in the form of the source terms, in the surface boundary conditions imposed, in the values of closure coefficients and also in the form of the damping functions. They both solve an equation for the isotropic component of the turbulent dissipation. The proposal to use this modified dissipation variable is due to Jones and Launder [START_REF] Jones | The prediction of laminarization with a two-equation model of turbulence[END_REF] who cited decisive computational advantages.

The model uses the following transport equations:

Turbulent kinetic energy:

∂ρk ∂t + ∂ρU j k ∂x j = ∂ ∂x j (µ + µ t σ k ) ∂k ∂x j + τ i j ∂U j ∂x j -ρε (3.49)
Dissipation rate:

∂ρε ∂t + ∂ρU j ε ∂x j = ∂ ∂x j (µ + µ t σ ε ) ∂ε ∂x j + ε k C ε1 τ i j ∂U j ∂x j -C ε2 ρε (3.50)
The eddy viscosity is specified as:

µ t = ρC µ k 2 ε (3.51)
relating the variables k and ε via a dimensionless constant C µ .

The constants of the model are:

C µ = 0.09, C ε1 = 1.44, C ε2 = 1.92, σ k = 1.0, σ ε = 1.3 (3.52)

The Menter kω SST model

The k-ω model is also a two-equation turbulence model. Instead of using the turbulent disspiation rate ε as the second transported variable, it uses the specific turbulent dissipation rate ω = ε k . The standard Wilcox [START_REF] Wilcox | Reassessment of the scale-determining equation for advanced turbulence models[END_REF] kω model is extremely accurate and robust in the near-wall region in comparison to the more commonly used k-ε models, but it suffers from high dependency on the chosen inlet freestream turbulence properties. Menter [START_REF] Menter | Influence of freestream values on kω turbulence model predictions[END_REF][START_REF] Menter | Zonal two equation kω turbulence models for aerodynamic flows[END_REF]

D ρ k Dt = τ i j ∂U j ∂x j -β * ρ k ω + ∂ ∂x j (µ + σ k µ t ) ∂ k ∂x j (3.53)
Specific dissipation rate:

D ρ ω Dt +U j = γ ν t τ i j ∂U j ∂x j -β ρ ω2 + ∂ ∂x j (µ + σ ω µ t ) ∂ ω ∂x j + 2(1 -F 1 ) ρσ ω 2 1 ω ∂ k ∂x j ∂ ω ∂x j (3.54)
with ν t = k/ ω. The new constants are generated using the relations:

φ = F 1 φ 1 + (1 -F 1 )φ 2 (3.55)
where φ 1 and φ 2 represent constants in the original Wilcox model and transformed kε respectively. The φ 1 constants (Wilcox kω) are:

σ k1 = 0.5, σ ω 1 = 0.5, β 1 = 0.0750 β * = 0.09, κ = 0.41, γ 1 = β 1 /β * -σ ω 1 κ 2 / β *
The φ 2 constants (standard Launder-Sharma kε) are:

σ k2 = 1, σ ω 2 = 0.856, β 2 = 0.0828 β * = 0.09, κ = 0.41, γ 2 = β 2 /β * -σ ω 2 κ 2 / β *
The other definitions are given by where y is the distance to the closest wall and CD kω is the positive part of the cross-diffusion term in the blended specific turbulence dissipation equation:

F 1 =
CD kω = max 2 ρσ ω 2 1 ω ∂ k ∂x j ∂ ω ∂x j ; 10 -20 (3.58)
The SST model is identical to the BSL model except the set of constants φ 1 and the definition of the eddy viscosity. The new constants (SST φ 1 ) are:

σ k1 = 0.85, σ ω 1 = 0.5, β 1 = 0.0750, a 1 = 0.31 β * = 0.09, κ = 0.41, γ 1 = β 1 /β * -σ ω 1 κ 2 / β *
In this case the eddy viscosity is defined as:

ν t = a 1 k max(a 1 ω; |Ω|F 2 (3.59)
|Ω| is the norm of the mean vorticity vector but the norm of the rate of strain tensor is now preferred [START_REF] Menter | Ten years of industrial experience with the sst turbulence model[END_REF]. F 2 is given by

F 2 = tanh(ar g 2 2 )
ar g 2 = max 2 k 0.009 ωy ; 500ν y 2 ω (3.60)

Menter also recommends to limit the production in the turbulent kinetic energy equation to ten times the dissipation to avoid excessive turbulence production, e.g. near stagnation points.

However, the standard eddy-viscosity models based on the Boussinesq hypothesis tend to overpredict eddy-viscosity that reduce the effect of re-entrant jet and two-phase structure shedding [START_REF] Sorgüven | Modified kω model for simulation of cavitating flows[END_REF][START_REF] Li | A modified sst kω turbulence model to predict the steady and unsteady sheet cavitation on 2d and 3d hydrofoils[END_REF]. These turbulence models are inadequate to correctly predict the dynamics of cavitation bubbles. Several solutions have been proposed and tested to reduce the eddy viscosity and improve the behavior of turbulence models. Reboud [START_REF] Reboud | Two phase flow structure of cavitation: experiment and modeling of unsteady effects[END_REF]] proposed an arbitrary modification by introducing an eddy viscosity limiter assigned as a function of density, f (ρ), instead of using the mixture density directly in the computation of the turbulent viscosity for the kε turbulence model.

µ t = f ρ C µ k 2 ε with f ρ = ρ v + ρ m -ρ v n ρ l -ρ v n-1 (3.61)
where n is a parameter set to 10.

The density function f (ρ) will be equal to ρ l and ρ v in the regions with pure liquid and vapor, but decreases rapidly in the region with a mixture of liquid and vapor.

A filter-based method (FBM) which combines the filter concept and the RANS model was investigated [Wu et al., 2005] [Tseng andWang, 2014] by imposing an independent filter scale, usually the grid size, in the computation of the eddy viscosity. Once the turbulence length scale is larger than the filter size, the eddy viscosity can be reduced by a linear filter function. These corrections have shown some success, but do not take into account the dynamics of small scales [START_REF] Coutier-Delgosha | Numerical simulation of the unsteady behaviour of cavitating flows[END_REF][START_REF] Goncalvès | Unsteady simulation of cavitating flows in venturi[END_REF][START_REF] Goncalvès | Numerical study of unsteady turbulent cavitating flows[END_REF].

The interplay between turbulence and cavitation regarding the unsteadiness and the structure of the flow is complex and not well understood. Moreover, there are less studies about the influence of the turbulence models on cavitating flow. In this study, the Reboud correction is implemented into three different turbulence models and simulated with different cavitation models. The end goal is to provide an insight into the interaction between the turbulence and cavitation models. to solve these test cases.

Interface movement in a uniform pressure and velocity flow

A discontinuity of volume fraction movement between two fluids in a uniform pressure and velocity flow at 100 m/s is considered. The discontinuity separates two nearly pure fluids from each other and is initially located at x = 0.5 m in a one meter length tube. There is liquid water in the left chamber and air in the right chamber. The uniform pressure is set equal to P = 10 5 Pa.

The fluid properties and initial condition for this test are given in Table 4.1.

WATER-AIR MIXTURE SHOCK TUBE

Water-air mixture shock tube

This test case, as proposed in [START_REF] Ansari | Numerical simulation of compressible two-phase flow using a diffuse interface method[END_REF], considers a one meter long shock tube containing two chambers which involves a discontinuity of the volume fraction at the location of x = 0.7 m. Each chamber contains a nearly pure fluid where the left chamber is filled with high-pressure fluids and the right one is filled with low-pressure fluids. The initial velocity is equal to 0 m/s. The fluid properties and initial condition for this test are given in Table 4.2.

Computations have been performed with a mesh of 1000 cells and with a time step ∆t = 10 -7 s. in [START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF]. In addition, the three-equation model was not able to predict well the phenomena.

Water-air mixture expansion tube |u| = 2 m/s

An expansion tube problem is considered with an initial velocity discontinuity located at the middle of the tube. This test consists in a one meter long tube filled with liquid water at atmospheric pressure and with density ρ l = 1150 kg/m 3 . A weak volume fraction of vapor α = 0.01 is initially added to the liquid. The initial discontinuity of velocity is set at 0.5 m, the left velocity is -2 m/s and the right velocity is 2 m/s. The solution involves two expansion waves. As gas is present, the pressure cannot become negative. To maintain positive pressure, the gas volume fraction increases due to the gas mechanical expansion and create a pocket. Liquid water is expanded until the saturation pressure is reached and then evaporation appears and quite small amount of vapor is created.

The parameters of the stiffened gas EOS and saturation values for densities are given in Table 4.3. The quantities have been evaluated with a saturation table at the reference temperature. The vapor pressure P vap = 51000 Pa. The solution obtained with the three-equation model is presented at time t = 3.2 ms in Figure 4.5. Results are compared with the two-fluid solution computed in [START_REF] Zein | Modeling phase transition for compressible two-phase flows applied to metastable liquids[END_REF]. The mesh contains 1000 cells. The time step is set to ∆t = 10 -8 s. The approximate HLLC Riemann solver and the four-equation model were not able to provide a solution and the Jameson-Schmidt-Turkel scheme [START_REF] Jameson | Numerical solution of the Euler equations by finite volume methods using Runge-Kutta time stepping schemes[END_REF] is used instead. Computation results of the void fraction and pressure profiles show large discrepancy with the two-fluid solution.

|u| = 100 m/s

The same conditions are used except regarding velocities which are set to u = -100 m/s on the left, and u = 100 m/s on the right. This case is stiffer than the previous one because of the high value of the initial velocity and evaporation is much more intense resulting in a large cavitation pocket where the gas volume fraction is close to 1. Computations are performed on a 1000-cell mesh with the time step set to ∆t = 10 -8 s.

The result obtained with the three-equation model is presented at time t = 1.5 ms in Figure 4.6. The void fraction profile is in good agreement with the two-fluid solution whereas the pressure simulated by the three-equation model is not able to capture the pressure drop inside the cavitation pocket. Again, the approximate HLLC Riemann solver and the four-equation model failed to solve this case.

Water-air shock bubble interaction

A cylindrical air bubble with an initial diameter D 0 = 6 mm is immersed in a water pool. Due to the symmetry of the problem the calculations are performed in a half-domain above the axis. The center of the bubble is located at (9, 0) mm in the computational domain of size 24 × 12 mm. The bubble is collapsed by a normal wave moving at Mach 1.72, initially located at abscissa x sh = 4 mm. The schematic diagram of the test case is given in Figure 4.7. The initial and post-shock conditions are shown in Table 4.4. The boundary conditions are the following: the top and bottom boundaries are assumed to be a wall and a symmetry axis, respectively. The left and right sides are assumed to be non-reflecting.

Simulations are performed using an uniform mesh composed by 1200 × 600 nodes and a time step ∆t = 10 -10 s. The number of nodes for a bubble diameter is 300.

The time evolution of the density gradient modulus (Schlieren-type representation) is plotted in Figure 4.8 from time t = 0.6 µs to t = 2.9 µs. After the water shock wave has collided with the bubble, a strong rarefaction wave is reflected backwards from the interface, and a weak shock wave is transmitted inside of the bubble (time t = 1.1 µs). At time t = 1.7 µs the incident water shock has traversed almost the full cavity width. The interaction between this shock and the expansion waves originating at the bubble surface has resulted in significant weakening and curvature of the incident shock. The shock inside the bubble propagates more slowly. Due to the pressure difference between both sides, the bubble is asymmetrically contracted and spreads laterally in the process. This change in shape is driven by vorticity generated at the edge of the bubble due to the passage of the wave which induces a jet of water along the axis of flow

Summary

To sum up, the validation of numerical method implemented in the NSMB solver for capturing the phenomenon of cavitation is illustrated in this chapter. There are four different test cases investigated and discussed with both the three-equation and four-equation models coupled with the HLLC scheme. The reuslts obtained from the previous test cases indicate that the implementation of these two cavitation models unfortunately could not be the cure-all and be generalized for all the test cases. In other words, it infers that there is the numerical instability for the implementation in the NSMB solver.

Since the NSMB is a huge solver, when the two cavitation models are implemented into the solver, there are about 20 % (one thousand) of subroutines which require to be modified or be coded. Therefore, it remains difficuilities and challenges for the modelling based on the viewpoint of the CFD. However, in order to achieve the academic goal of this study, turbulence and cavitation, another free open source software, OpenFOAM, is adopted. With the built-in solver of the OpenFOAM, interphaseChangeFoam, the phenomenon of cavitation will be studied with the 4°Venturi geometry and will be presented in the next chapter.

Venturi 2D

Experimental conditions

The Venturi type test section of the CREMHyG cavitation small tunnel was sized and designed to

Mesh and computational set-up

The grid is a H-type topology. It contains 251 nodes in the flow direction and 81 nodes in the orthogonal direction. A special contraction on the mesh is applied in the main flow direction just after the throat to better simulate the two-phase flow area, as illustrated in Figure 5.3. Various computations were performed by varying the cavitation model and the turbulence model.

An overview that includes all test cases of the Venturi 2D case is shown in The empirical values of the three cavitation models in OpenFOAM solver are specified in

C c = 10 C v = 8000 Merkle U ∞ = 10.8(m/s) t ∞ = 0.023(s) C c = 80 C v = 0.001 SchnerrSauer n = 1.6 × 10 13 (m -3 ) dnuc = 2 × 10 -6 (m) C c = 1 C v = 1
Table 5.3: Empirical values of the cavitation models

Results for different turbulence models

The calculations are done by using the three cavitation models that are implemented in the solver interPhaseChangeFoam of OpenFOAM with three different turbulence models. The Reboud correction is applied to these turbulence models with the exponent value of n = 10. All numerical values are obtained by a time-averaged statistical treatment over a simulation time of 5 s. Figure 5.9 illustrates the evolution of the void ratio and longitudinal velocity profiles for the numerical and experimental results from station 1 to 5 with Merkle cavitation model. At station 1, close to the throat, the cavity thickness is similar by all simulations and is underestimated. Computations with the KEReboud and KWReboud turbulence models over-estimate the maximum value of the void ratio: the discrepancy with the experimental value is around 7% and 10%. On the other hand, computation with the SAReboud turbulence model under-estimates the maximum value of the void ratio (around 7%). At station 2, all computations capture well the cavity thickness with an increase of it from station 1. However, computations with the KEReboud and SAReboud turbulence models under-estimate the maximum value of the void ratio. The discrepancy with the experimental value is about 10%. For the SAReboud turbulence model, the void ratio value at the wall is extremely low in relation to the experimental value. On the contrary, the computation with the KWReboud turbulence model is over-predict the the maximum value of the void ratio (around 4%). For the velocity profiles, at station 1 and 2, all computations over-estimate the maximum value of the longitudinal velocity. The re-entrant jet phenomena are observed with all turbulence models in station 1 and the KEReboud and SAReboud turbulence models at station 2.

From station 3, the re-entrant jet becomes noticeable from the velocity measurement. The KEReboud and SAReboud turbulence models reproduce well the effect of recirculating behavior with a re-entrant jet extending through half the sheet thickness at station 3, 4 and 5. However, simulation with the KWReboud turbulence model does not present the thickness of the recirculating area and there is no re-entrant jet. The computation with the KEReboud turbulence model captures better the intensity of the recirculating zone near the wall at station 3 and 4. At station 5, both the KEReboud and SAReboud turbulence models have the similar intensity prediction. As for the void ration profiles, from station 3 to 5, computations with the KEReboud and SAReboud turbulence models provide a better prediction of the void ratio values and the thickness of the cavitation compared to the KWReboud one. In addition, the wall value of the void ratio given by the SAReboud turbulence model at station 3 and 4 has good agreement with the experimental value. The computation with the KWReboud turbulence model largely overestimates the void ratio value.

The dimensionless wall pressure profiles are plotted in Figure 5.10 versus the distance xx inl et . For all computations, the pressure remains at an almost constant value P vap in the cavity. Computations with the KEReboud and SAReboud turbulence models show that the recompression starts from station 4 and re-compress slowly afterward in comparison with the experimental data. On the other hand, the KWReboud turbulence model keeps the pressure equal to the vaporization pressure far downstream. The re-compression is lower for all simulation downstream station 5 in relation to the experimental data.

The Root Mean Square (RMS) wall pressure fluctuations are plotted in Figure 5.11 versus the distance xx inl et . All models predict different behavior for the pressure fluctuations. The .13 illustrates the evolution of the void ratio and longitudinal velocity profiles for the numerical and experimental results from station 1 to 5 with SchnerrSauer cavitation model. At station 1, inside the attached cavity sheet, all computations capture the similar cavity thickness but is under-predicted. Computations with the KEREboud and KWReboud turbulence models over-estimate the maximum value of the void ration: the discrepancy with the experimental value is 6% and 10%. At station 2, all simulations present well the cavity thickness with a correct estimation of it. Computations with the KEReboud and SAReboud turbulence models under-estimate the maximum value of the void ratio (around 10%) whereas the KWReboud turbulence model over-estimate it (around 4%). For the velocity profiles, at station 1 and 2, all computations over-predict the maximum value of the longitudinal velocity. Moreover, the re-entrant jet phenomena are obtained with all turbulence models which are not observed in the experimental value and the effect is even stronger with the KWReboud turbulence model.

From station 3 to 5, the re-entrant jet becomes noticeable from the velocity measurement. The KEReboud and SAReboud turbulence models capture well the effect of recirculating behavior with a re-entrant jet except that the thickness of the recirculating area is over-estimated. In addition, the intensity of the recirculating zone captured by the KEReboud turbulence model is in good agreement with the experiment whereas it is under-estimated by the SAReboud turbulence model.

The KWReboud turbulence model is also capable of predicting the thickness of the recirculating area but largely over-estimate the intensity of the recirculating zone. As for the void ratio profiles, at station 3, all computations provide a better prediction of the thickness of the cavity sheet as compared to the experimental data. The void ratio values are over-predicted by all computations especially for the KWReboud turbulence model. At station 4 and 5, the void ratio values and thickness of the cavitation are over-estimated by all computations. For the SAReboud turbulence model, the void ratio value at the wall is in better agreement with the experimental value.

The dimensionless wall pressure profiles are plotted in Figure 5.14 versus the distance xx inl et . For all computations, the pressure remains at an almost constant value P vap in the cavity. Computations with the KEReboud and SAReboud turbulence models show that the re-compression starts from station 4 and re-compress slowly after in comparison with the experimental data. On the other hand, the KWReboud turbulence model keeps the pressure equal to the vaporization pressure far downstream. The re-compression is lower for all simulation downstream station 5 in relation to the experimental data. Figure 5.17 illustrates the evolution of the void ratio and longitudinal velocity profiles for the numerical and experimental results from station 1 to 5 with the kε turbulence model with the Reboud correction. At station 1, the similar cavity thickness is estimated by all models and is under-predicted. Computations with the Merkle and SchnerrSauer cavitation models over-estimate the maximum value of the void ratio (around 6%) wheres the Kunz cavitation model under-estimate it (around 6%). At station 2, all models estimate well the cavity thickness but under-estimate the void ratio composition. At station 3, all models estimate still well the cavity thickness but over-estimate the void ratio composition. At station 4 and 5, all models over-estimate both the cavity thickness and the void ratio composition. For the velocity profiles, at station 1 and 2, all models present the recirculating area which is not observed in the experiment.

Further downstream the re-entrant jet phenomenon is well observed in the experiment. From At station 4 and 5, all computations over-predict the cavity thickness and composition. For the velocity profiles, at station 1 and 2, computations with all cavitation models present the reentrant phenomenon which is not observed in the experiment. At station 3, 4 and 5, all cavitation models capture the re-entrant phenomenon but the recirculating area and the intensity of the recirculating zone near the wall are under-estimated.

The dimensionless wall pressure profiles are plotted in Figure 5.26 versus the distance x -

Venturi 3D

Computations of the 3D Venturi geometry are performed with Kunz cavitation model with the KEReboud and KWReboud turbulence models. The 3D mesh was built by extruding the 2D mesh with 81 nodes in the cross direction z. As the test section is a square, the same grid evolution was applied in the y and z directions. A view of the mesh is presented in 5. Comparisons between 2D and 3D simulations concern time-averaged quantities extracted on the mid-span plan. Figure 5.30 presents the evolution of the void ratio and longitudinal velocity profiles for the numerical and experimental results from station 1 to 5 with Kunz cavitation model. At station 1, close to the throat, the cavity thickness obtained by all simulations is underestimated. The results of 2D and 3D predicted by the KWREboud turbulence model are almost the same. At station 2, all computations capture well the cavity thickness with an increase of it from station 1. However, all computations under-estimate the maximum value of the void ratio.

With the 3D effect, the void ratio value at the wall obtained by the KWReboud turbulence model is extremely low in relation to the experimental value. For the velocity profiles, at station 1 and 2, all computations over-estimate the maximum value of the longitudinal velocity. In addition, there shows the re-entrant jet phenomena with the KEReboud turbulence model which is not observed in the experiment.

For station 3, 4 and 5, the re-entrant jet becomes noticeable from the velocity measurement.

Simulation with the KWReboud turbulence model in 2D geometry does not present the thickness of the recirculating area at station 3 and the thickness of the recirculating area is under-estimated.

On the contrary, the computation in 3D geometry gives better prediction in it. Compare to the KEReboud turbulence model in 2D geometry, the computation with 3D geometry over-predicts the thickness of the recirculating area. The intensity of the recirculating zone near the wall is under-estimated by all calculations no matter in 2D or 3D. Regarding the void ratio profiles, an over-estimation for both the void ratio values and the thickness of the cavitation is observed by all computations. The result obtained with the KWReboud turbulence model in 3D geometry is in better agreement with the experimental data.

The dimensionless wall pressure profiles are plotted in Figure 5.31 versus the distance xx inl et . It showslarge discrepancy between the 2D and the 3D computations. The KWReboud turbulence model on 3D geometry presents an almost constant value of pressure in the cavity, and the re-compression is well captured downstream station 5 compared to the experimental data. The KEReboud turbulence model on 3D geometry re-compresses slowly and has an underestimation of the pressure downstream the cavity.

The Root Mean Square (RMS) wall pressure fluctuations are plotted in Figure 5.32 versus the distance xx inl et . For the 2D computations, the peak of fluctuations is obvious although the peak location and amplitude are higher and lower for calculations with the KEReboud and the KWReboud turbulence models in comparison with the experimental data. In addition, the decrease of the fluctuations level in the wake of the cavitation sheet is better captured by the KWReboud turbulence model both on 2D and 3D geometry. This is not the case for the KEReboud turbulence model since the level of pressure fluctuations is largely overestimated in the wake. Figure 5.33 compares time-averaged profiles of the viscosity ratio µ t /µ, at the five station, obtained with Kunz cavitation model. At station 1 and 2 inside the cavity, the effect of 3D geometry is less important. Among the results obtained from the simulation which were compared to the experimental data, it is the Kunz's cavitation model coupled with the kω SST turbulence model that could have a better prediction for the 4°Venturi geometry. Nevertheless, the cavity length of all the models unfortunately was over-predicted by all the simulation. In addition, the 3D effect did not much improve the prediction either according to the obtained numerical results. The main reason for the cause might be that the incompressible solver was not suitable for this kind of internal geometry. The three-equation model was closed with a sinusoidal barotropic EOS. The mixture of stiffened gas EOS was applied to the four-equation model. The proposed models were validated through various inviscid test cases including the interface movement problem, water-air shock tube and expansion tube and shock-bubble interaction. The capability to obtain correct solutions of these test cases has been investigated. The reuslts obtained from the test cases indicate that the implementation of these two cavitation models unfortunately could not be the cure-all and be generalized for all the test cases. Although the validations showed the ability of models to simulate the cavitation development, the two models are still suffered from the issue of numerical instability. The main difference between these two models is that the three-equation model has the assumption of complete thermodynamic equilibrium between phases; therefore, it could explain the discrepancies existed for the test cases above. Since the implementation and validation in the NSMB solver had already taken too much time, in order to achieve the purposes of this study, which are the turbulence and cavitation, another free open source software, OpenFOAM, was adopted to perform the cavitating flows on the Venturi geometry. to reduce the eddy viscosity in order to capture the re-entrant jet dynamics. Results showed that the use of an eddy-viscosity limiter lets the model correctly simulate unsteady behaviors of the sheet, however large discrepancies occur between models and the effect of reduction is not strong enough. Generally the three cavitation models were able to reproduce the re-entrant jet phenomenon but the cavity length was over-predicted. This might be due to the calibration problem of the mass transfer term of the condensation rate and vaporization rate coefficient or the lack of thermodynamic coherence. Also, the impact on the value of the exponent n used in this correction is necessary to be investigated. Besides, the interPhaseChangeFoam is an incompressible solver which is less capable of solving the type of internal geometry. Hence, the speed of sound can be several orders of magnitude higher in the liquid phase than in the two-phase mixture. Thus for low-speed applications, the numerical method must be able to properly and efficiently simulate both incompressible and compressible flow areas. In terms of computational methods, the application of a compressible formulation to simulate low speed cavitating flows results in poor convergence and erroneous calculations. To achieve this goal, a preconditioned method is necessary. It is based on the modification of the derivative term by a premultiplication with a suitable preconditioning matrix. The physical acoustic waves are replaced by pseudo-acoustic modes that are much closer to the advective velocity, reducing the stiffness and enhancing the convergence. Therefore, the preconditioned method can provide both efficiency and accuracy over a wide range of Mach numbers.

The preconditioned 1-D Euler equations with the primitive variables W = (P, u, s), where s is the entropy, can be expressed as:

P -1 e ∂W ∂t + A e ∂W ∂x = 0
where,

P e =     β 2 0 0 0 1 0 0 0 1     , P -1 e =     1 β 2 0 0 0 1 0 0 0 1     , A e =     u ρ c 2 0 1 ρ u 0 0 0 u     1 β 2 P t + uP x + ρ c 2 u x = 0 u t + uu x + 1 ρ P x = 0 s t + us x = 0
The equations can be rewritten in vector form as:

    1 β 2 0 0 0 1 0 0 0 1         P u s     t +     u ρ c 2 0 1 ρ u 0 0 0 u         P u s     x = 0
with the eigenvalues of u and λ ± = 1 2 u(1 + β 2 ) ± (β 2 -1) 2 u 2 + 4β 2 c 2 . For the conservative variables w c = (ρ, ρu, ρE), the corresponding form is:

P -1 c ∂w c ∂t + A c ∂w c ∂x = 0
where the preconditioning matrix

P -1 c = (∂w c /∂W)P -1 e (∂W/∂w c ) = R -1 P -1 e R and A c = ∂F c ∂x = R -1 A e R

Differential of the primitive variables with respect to the conservative variables

Velocity

du = 1 ρ d(ρu) - u ρ dρ therefore,          ∂u ∂ρ = -u ρ ∂u ∂ρu = 1 ρ ∂u ∂ρE = 0 Energy e = E - u 2 2 ⇒ de = dE -udu
also, with d(ρE) = ρdE + Edρ and ρudu = ud(ρu) -u 2 dρ, we obtain that

de = 1 ρ d(ρE) - E ρ dρ - u ρ d(ρu) + u 2 ρ dρ = 1 ρ d(ρE) - u ρ d(ρu) + u 2 -E ρ dρ therefore,          ∂e ∂ρ = u 2 -E ρ ∂e ∂ρu = -u ρ ∂e ∂ρE = 1 ρ Pressure d (
+ hdρ = Adρ + BdP ⇒ dP = ρT B ds + h -A B dρ = ∂P ∂ρ s dρ + ∂P ∂s ρ ds ⇒    ∂P ∂ρ s = c 2 = h-A B ∂P ∂s ρ = ρT B dP = ρT B ds + h -A B dρ = ρT B ds + c 2 dρ ⇒ ds = B ρT dP - Bc 2 ρT dρ with dP = 1 B d(ρE) - u B d(ρu) + u 2 2 -A 1 B dρ therefore, ds = 1 ρT d(ρE) - u ρT d(ρu) + u 2 2 -A 1 ρT - Bc 2 ρT dρ ⇒            ∂s ∂ρ = u 2 2 -A -Bc 2 1 ρT = 1 ρT u 2 2 -h where ∂P ∂ρ s = c 2 = h-A B
   =      (γ -1) u 2 2 -h + c 2 γ-1 -u(γ -1) γ -1 -u ρ 1 ρ 0 1 ρT u 2 2 -h -u ρT 1 ρT      R -1 = ∂w c ∂W =    
    =      1 c 2 0 - (γ-1)ρT c 2 u c 2 ρ - (γ-1)uρT c 2 h+ u 2 2 c 2 ρu - γ-1 c 2 ρT u 2 2 + h -c 2 γ-1     
The preconditioning matrix

P -1 c = I d + γ -1 c 2 1 β 2 -1      u 2 2 -h + c 2 γ-1 -u 1 u 2 2 -h + c 2 γ-1 u -u 2 u u 2 2 -h + c 2 γ-1 H -uH H     
with the definition of total enthalpy H = h+ u 2 2 and c

2 γ-1 = h-A, we obtain that u 2 2 -h+ c 2 γ-1 = u 2 2 -A P -1 c = I d + 1 β 2 -1 1 h -A     u 2 2 -A -u 1 u 2 2 -A u -u 2 u u 2 2 -A H -uH H    

B.3 Computational results

All cases for the BSCW were run using three meshes, namely the coarse, medium and fine resolution multi-block structured meshes supplied by the workshop organizing committee. Different turbulence models, the Spalart-Allmaras model with Quadratic Constitutive Relation, 2013 version (SA QCR 2013) [START_REF] Mani | Predictions of a supersonic turbulent flow in a square duct[END_REF], the Menter Shear Stress Transport model (kw SST) [START_REF] Menter | Influence of freestream values on kω turbulence model predictions[END_REF] and the Chien k-epsilon model (kec) [START_REF] Chien | Predictions of channel and boundary-layer flows with a low-reynolds-number turbulence model[END_REF], were applied in the simulation.

Figure 1

 1 Figure 1.1: Phase diagram of water.

  Figure 1.2: Damage of vane by cavitation.

Figure

  Figure 1.3: Traveling cavitation.

Figure

  Figure 1.4: Sheet cavitation.

Figure 1 . 5 :

 15 Figure 1.5: Cloud cavitation.

  Figure 1.6: Vortex cavitation.

From

  the thermodynamic point of view, two state variables are sufficient to represent the thermodynamic state of a fluid. The main relationships existing in the literature are : This assumption leads to a very simplified state law : ρ = ρ 0 and C p = C v = C which are the specific heats at constant pressure and constant volume respectively. This equality leads to the following relation between the internal energy and the CHAPTER 2. REVIEW OF CAVITATION MODELING temperature : de = CdT

  27) where U V ,n = u.n with n = ∇α L |∇α L | The normal interfacial velocity, U I,n , is zero in steady calculation. This model is called Sharp Interfacial Dynamics Model (IDM).

Figure 2

 2 Figure2.1: Representation of a vapourous cavity[Senocak et Shyy, 2004]

  1000 and c D = 0.44 otherwise (2.54) R e and R e P are the Reynolds number relative to the flow and the Reynolds number relative to the bubbles respectively. This is one of the first cavitation model to simulate the slip between phases. However, the evolution of this model by solving the energy conservation equations seems not yet to be achieved. b/ Saurel and Le Métayer model with seven equations (2003)

  e the Reynolds number based on the radius of the bubble : R e = U V -U L d ν L and P r L the Prandtl number of liquid : P r L = ν L a L In these equations, λ L is the thermal conductivity of liquid , C p V the specific heat of vapor at constant pressure, ∆t the iterative time step, ν L the kinematic viscosity of liquid and a L the thermal diffusivity of liquid. The momentum transfer is based on a term due to mass transfer and a drag term. The results are presented to a 3D cold water flow through a diaphragm. The calculations are compared with the test data from the case named SUPER MOBY DICK. This study makes it possible to bring a new modeling of the interfacial transfer for the cavitation based on those of ebullition.

  are implemented and tested with the interface movement and shock-bubble intercation for the NSMB solver. In addition, the transport equation based method for the void ratio including the source terms for vaporization and condensation in the free, open source software OpenFOAM (Open source Field Operation And Manipulation) are also presented for the Venturi geometry. C H A P T E R 3 NUMERICAL SOLVER T wo numerical codes used for the simulation of this study are briefly described in this chapter. The Navier-Stokes Multi-Block (NSMB) solver is a numerical software developed within an european consortium solving the finite volume Navier-Stokes equations. NSMB is a multi-block structured solver and parallelized able to solve the steady or unsteady Navier-Stokes equations in their compressible or incompressible version. Open source Field Operation And Manipulation (OpenFOAM) is a free, open source software for computational fluid dynamics (CFD). OpenFOAM is a Finite Volume Method (FVM) based numerical solver for solving systems of transient transport equations. Plenty of solvers available for a wide range of domains such as incompressible, compressible, multiphase, combustion, etc . . . 3.1 NSMB Navier-Stokes Multi-blocks (NSMB) was based on a structured multi-block Euler code (EULMB) developed at Swiss Federal Institute of Technology in Lausanne (EPFL, "École Polytechnique Fédérale de Lausanne") in 1989 with the support from the European Centre for Research and Advanced Training in Scientific Computation (CERFACS, "Centre Européen de Recherche et de Formation Avancée en Calcul Scientifique") and the Royal Institute of Technology (KTH). Originally developed by Jan Vos (CFS-Engineering, Lausanne) in 1989, NSMB was developed from 1992 to the end of 2003 in the NSMB consortium, which included several universities (EPFL, Lausanne, Switzerland; SERAM ("Société d'études et de recherches de l'École nationale supérieure d'arts et métiers"), Paris, France; Institute of Fluid Mechanics of Toulouse (IMFT, "Institut de Mécanique des Fluides de Toulouse"), Toulouse, France; KTH, Stockholm, Sweden), a research institution (CERFACS, Toulouse, France), and industrial partners EADS-France (Airbus CHAPTER 3. NUMERICAL SOLVER France and EADS Space Technologies), SAAB Military Aircraft Engineering and CFS. Since 2004, NSMB is further developed by the EPF-Lausanne, ETH-Zrich, ICUBE-Strasbourg, IMFT-Toulouse, Polytechnic University of Munich, the military University of Munich, CFS Engineering and RUAG Aerospace. In addition to these groups, NSMB is still used by Airbus France, EADS-ST and KTH. The NSMB solver is a code which is parallelized in MPI (Message Passing Interface) and solves the steady or unsteady Navier-Stokes equations in their compressible or incompressible version on multi block structured grids by means of finite volume method. It provides a wide range of numerical schemes both for spatial and temporal discretisation. There are for example the central schemes (2nd and 4th order with artificial dissipation) and upwind schemes (Roe, AUSM, Van Leer, Harten . . . ) available for the spatial discretisation. Within this solver turbulence can be treated in various ways: LES (Smagorinski, structure functions Lesieur et al, Ducros et al . . . ), turbulence models from zero equations (Baldwin-Lomax, Granville . . . ), one equation (Spalart-Allmaras and several variants), two equations linear models (kε, kω, SST . . . ), nonlinear models or explicit (EARSM) and the RSM model. For all these models, their RANS-LES (DES, DDES, IDDES, WMLES) hybrid variants have been implemented. The SAS model based on Menter kω SST has also been implemented. Chimera Methodologies (overlapping meshes) and

  source Field Operation And Manipulation (OpenFOAM) is a free, open source software for computational fluid dynamics (CFD). It is owned by the OpenFOAM Foundation and licenced under the GNU General Public Licence (GPL) that gives users the freedom to modify and redistribute the software and a guarantee of continued free use within the terms of the licence. The codes are written in C++ programming language in an object-oriented manner to solve ordinary differential equations (ODEs) and partial differential equations (PDEs). The correspondence between the implementation and the original equation is clear due to the high level programming. This feature makes users straightforward to modify or mimic the existing solvers. As a result, this provides OpenFOAM with good extensibility qualities. Another distinguishing feature of OpenFOAM is that it can be used in massively parallelism through domain decomposition method, where the computational domain is split into a number of subdomains, one for each processor. Each processor receives a separate distribution of the complied code to be run on each subdomain. For the communication between processors, the Message Passing Interface (MPI) is used. In addition, it provides with plenty of pre-and post-processing utilities for users to perform. OpenFOAM is a Finite Volume Method (FVM) based numerical solver for solving systems of transient transport equations. Regarding the finite-volume discretization, a variety of discretization schemes are implemented for the temporal, convection, diffusion and source terms in the transport equations. Meanwhile, there are plenty of solvers available for a wide range of domains such as incompressible, compressible, multiphases, combustion, etc . . . The above-mentioned advantages give OpenFOAM solvers a great capabilities and extensibility.

  3.3. TURBULENCE CLOSURESThe solver interPhaseChangeFoam of OpenFOAM is selected to simulate the cavitation. It is a solver for two incompressible, isothermal immiscible fluids with phase-change (e.g. cavitation) and uses a volume-of-fluid (VOF) phase-fraction based interface capturing approach. The momentum and other fluid properties are of the "mixture" and a single momentum equation is solved. The set of phase-change models provided are designed to simulate cavitation but other mechanisms of phase-change are supported within this solver framework. The solver includes Kunz, Merkle and SchnerrSauer cavitation models. It uses the Pressure-Implicit with Splitting of Operators (PISO) algorithm to solve the Navier Stokes equations by first solving the momentum equations with pressure from the previous time step followed by solving the pressure equation for the new velocity field followed by velocity correction. Details on the numerical principles and specific implementation can be referred to the documentation of OpenFOAM (http://www.openfoam. com).

  CHAPTER 3. NUMERICAL SOLVER term on the right hand side affects only the normal stresses which equals to twice the turbulent kinetic energy. Three turbulence models, the one-equation Spalart-Allmaras model, the two-equation kε model and the Menter kω SST modle, are chosen for the study and briefly described in the following section.

  attacked this problem by first transforming the kε model into a kω type formulation before introducing a blending function dependent, among other things, on distance from the nearest wall. Compared to the original k -ω model, the differences are the values of the model constants and the presence of an additional (cross-diffusion) term (Menter discards a small additional diffusion term during the transformation). The two models are then combined by multiplying the original kω model by a function F 1 , the transformed model by (1 -F 1 ), and then adding. The result is a model which keeps the robust and accurate near-wall formulation of the original Wilcox kω model and improves freestream independence through use of the kε model in the outer part of the boundary layer. Futher to this baseline (BSL) model, Menter then added a shear-stress transport correction to form the kω model. The Menter kω SST model is known to be particularly capable of capturing flow fields featuring large separated shear layers. The equations governing the kω SST model are presented below.

  the validation of numerical method implemented in the NSMB solver for capturing the phenomenon of cavitation. Four test cases including the one-and twodimensional compressible two-phase flows with interface conditions are considered. Both the three-equation and four-equation model coupled with the HLLC scheme have been proposed

Figure 4 . 1 :

 41 Figure 4.1: Interface movement discontinuity problem. Void fraction and pressure profiles by 3-equation model (symbols) and the exact solution (solid line).

Figure 4

 4 Figure 4.2: Interface movement discontinuity problem. Void fraction and pressure profiles by 4-equation model (symbols) and the exact solution (solid line).

  Numerical solutions computed with the 3-equation and 4-equation model at 240 µs are shown in Figure 4.3 and 4.4 respectively. In this test case, strong pressure waves are propagated. The obtained result with the four-equation model are in close agreement with solutions presented

Figure 4

 4 Figure 4.4: Water-air mixture shock tube problem. Density, pressure, velocity and void fraction profiles by 4-equation model (symbols) and the exact solution (solid line).

Figure 4 . 5 :

 45 Figure 4.5: Water-air mixture expansion tube problem |u| = 2 m/s. Void fraction and pressure profiles by central scheme with 3-equation model (symbols) and 7-equation model (solid line).

Figure 4

 4 Figure 4.6: Water-air mixture expansion tube problem |u| = 100 m/s. Void fraction and pressure profiles by central scheme with 3-equation model (symbols) and 7-equation model (solid line).

Figure 4 . 7 :

 47 Figure 4.7: Initial situation for the shock bubble interaction D 0 = 0.006 m and M sh = 1.72.

Figure 4

 4 Figure 4.8: Water-air shock bubble interaction. Time evolution of the density gradient.

Figure 4

 4 Figure 4.9: Water-air shock bubble interaction. Time evolution of the pressure (in bar).

Figure 4 .

 4 Figure 4.10: Water-air shock bubble interaction. Time evolution of the axial velocity (in m/s).

Figure 4 .

 4 Figure 4.11: Water-air shock bubble interaction. Time evolution of the vertical velocity (in m/s).

  is investigated in this chapter for cavitating flow. First, the investigation is conducted on a 2D Venturi geometry with available experimental data tested in the cavitation tunnel of CREMHyG (Centre de Recherche et d'Essais de Machines Hydrauliques de Grenoble). Second, the 3D effect of the same geometry is considered. The simulation of cavitating flow is carried out by the free software OpenFOAM. The built-in solver, interPhaseChangeFoam, is used for the computation. For the Venturi case, three turbulence models, i.e. the Spalart-Allmaras, kε and kω SST models with the Reboud correction are considered together with three cavitation models of interPhaseChangeFoam.

  Figure 5.1: Schematic view of the Venturi profile.

Figure 5 . 2 :

 52 Figure 5.2: Photograph of the cavity.

Figure 5 . 3 :

 53 Figure 5.3: Enlargement of the mesh near the Venturi throat.

Figure 5 .Figure 5 . 4 :

 554 Figure 5.4 illustrates the different cavities where the time-averaged void ratio is plotted. All the cavitation and turbulence models except for Merkle and SchnerrSauer cavitation models coupled with KWReboud turbulence model show an attached cavity sheet with a large reentrant jet and with the presence of small clouds of mixture in the closure region of the sheet. The configuration of quasi-stable sheet is observed. Among them, only Kunz cavitation model coupled with the KWReboud turbulence and SchnerrSauer cavitation model coupled with the

Figure 5 .

 5 Figure 5.5 presents the evolution of the void ratio and longitudinal velocity profiles for the numerical and experimental results from station 1 to 5 with Kunz cavitation model. At stations 1 and 2, inside the attached cavity sheet, a relative strong effect of the vaporization phenomenon is clearly represented from the void ratio profile. The void ratio values of the first two stations are 0.9 and 0.95 near the wall according to experiments. At station 1, computations with the KEReboud and SAReboud turbulence model are close and with the discrepancy of about 15% compared to the experimental value. On the contrary, computation with the KWReboud turbulence model over-estimates the maximum value of the void ratio (around 10%). The numerical cavity thickness is under-estimated for all turbulence models. At station 2, the distribution is similar to that obtained for station 1, with an increase in the sheet thickness. Computations with the KEReboud and SAReboud turbulence models largely under-estimate the maximum value of the void ratio (around 20%), as observed for station 1. The computation with the KWReboud turbulence model is in better agreement with the experimental data. For the velocity profiles, at station 1 and 2, all computations over-estimate the maximum value of the longitudinal velocity. In addition, there shows the re-entrant jet phenomena with the KEReboud and SAReboud turbulence models which is not observed in the experiment.Further downstream, for stations 3, 4 and 5, the re-entrant is observed on the velocity measurement. The experiment observation indicates a recirculating behavior with a re-entrant jet extending roughly half the sheet thickness. However, the effect is not very evident especially for the computation with the KWReboud turbulence model. In addition, the thickness of the recirculating area is largely under-estimated by the KWReboud turbulence model at station 3 and 4. The intensity of the recirculating zone near the wall is under-estimated by all calculations.Regarding the void ratio profiles, an over-estimation for both the void ratio values and the thickness of the cavitation is observed by all computations. Computations with the KEReboud and SAREeboud turbulence models are similar. Both models have the similar sheet thickness and the void ratio values. For the KWReboud turbulence model, the void ratio value at the wall is in very good agreement with the experimental value although the maximum value of it and the thickness of the sheet are greatly over-estimated.The dimensionless wall pressure profiles are plotted in Figure5.6 versus the distance x-x inl et .The first five data are located inside the cavity (where the void ratio and velocity profiles are measured). Both the KEReboud and SAReboud turbulence models have a slow re-compression which results in an underestimation of the pressure downstream the cavity. The KWReboud turbulence model presents an almost constant value of pressure in the cavity, but the value keeps the same as the vaporization pressure further downstream compared to the experimental data.The Root Mean Square (RMS) wall pressure fluctuations are plotted in Figure5.7 versus the distance xx inl et . The pressure fluctuation is divided by the time-averaged pressure P av .Experimental data indicate an augmentation of pressure fluctuations at the end of the sheet cavity, with a peak located at station 5. The peak of pressure fluctuation predicted by the

Figure 5 .

 5 Figure 5.8 compares time-averaged profiles of the viscosity ratio µ t /µ, at the five station, obtained with Kunz cavitation model. The effect of the Reboud correction is not obvious for both the KEReboud and SAReboud turbulence models. The KWReboud turbulence model induces a large reduction of the ratio in the sheet at station 1 to 3.

Figure 5 . 5 :Figure 5 . 6 :Figure 5

 55565 Figure 5.5: Time-averaged void ratio (left) and velocity (right) profiles from station 1 to -Kunz model comparison

Figure 5 .Figure 5 .Figure 5 .

 555 Figure 5.9: Time-averaged void ratio (left) and velocity (right) profiles from station 1 to 5 -Merkle model comparison

The

  Root Mean Square (RMS) wall pressure fluctuations are plotted in Figure 5.15 versus the distance xx inl et . Different behaviors for the pressure fluctuations are obtained by the turbulence models. The peak position predicted by the KWReboud turbulence model is present downstream station 5 and the amplitude of the fluctuation peak is underestimated. The pressure fluctuations simulated with both the KEReboud and SAReboud turbulence models are too high.

Figure 5 .Figure 5 .Figure 5 .Figure 5 .Figure 5 .

 55555 Figure 5.16 compares time-averaged profiles of the viscosity ratio µ t /µ, at the five station, obtained with SchnerrSauer cavitation model. Computational results with the KEReboud turbu-

  station 3 to 5, computations with the Merkle and SchnerrSauer cavitation models reproduce the recirculating area although over-estimate the thickness of it. The Kunz cavitation model does not capture the recirculating area. The dimensionless wall pressure profiles are plotted in Figure 5.18 versus the distance xx inl et . The Kunz cavitation model has a slow re-compression which results in an underestimation of the pressure downstream the cavity. Computations with the Merkle and SchnerrSauer cavitation models show that the re-compression starts from station 4 and re-compress slowly afterward in comparison with the experimental data. The Root Mean Square (RMS) wall pressure fluctuations are plotted in Figure 5.19 versus the distance xx inl et . The pressure fluctuation simulated with the Kunz cavitation model is too high. The Merkle and SchnerrSauer cavitation models provide fluctuations in better agreement with experimental data inside the cavity but not in the re-compression area.

Figure 5 .Figure 5 .Figure 5 .Figure 5 .Figure 5 .Figure 5 .

 555555 Figure 5.20 compares time-averaged profiles of the viscosity ratio µ t /µ, at the five station, obtained with the KEReboud turbulence model. With both the Merkle and SchnerrSauer cavitation models, the reduction effect is similar.

Figure 5 .

 5 Figure 5.24 compares time-averaged profiles of the viscosity ratio µ t /µ, at the five station, obtained with the KWReboud turbulence model. Computational results show that the reduction effect with the SchnerrSauer cavitation model is less evident.

Figure 5 .Figure 5 .Figure 5 .Figure 5 .Figure 5 .

 55555 Figure 5.21: Time-averaged void ratio (left) and velocity (right) profiles from station 1 to 5 -kω SST model with the Reboud correction comparison

x

  inl et . The Kunz cavitation model has a slow re-compression which results in an underestimation of the pressure downstream the cavity. The Merkle and SchnerrSauer cavitation models give the similar results with an almost constant value of pressure in the cavity. The re-compression starts from station 4 and re-compress slowly afterward in comparison with the experimental data. The Root Mean Square (RMS) wall pressure fluctuations are plotted in Figure 5.27 versus the distance xx inl et . The pressure fluctuation simulated with the Kunz and Merkle cavitation models are too high. The SchnerrSauer cavitation model provide fluctuations in better agreement with experimental data inside the cavity but not in the re-compression area.

Figure 5 .

 5 Figure 5.28 compares time-averaged profiles of the viscosity ratio µ t /µ, at the five station, obtained with the SAReboud turbulence model. All the cavitation models give the similar results except at station 4 and the computed cavity is nearly steady.

Figure 5 .Figure 5 .Figure 5 .

 555 Figure 5.25: Time-averaged void ratio (left) and velocity (right) profiles from station 1 to 5 -Spalart-Allmaras model with the Reboud correction comparison

Figure 5 .

 5 Figure 5.29: View of the 3D mesh composed of 251 nodes in the flow direction and 81 nodes in each transversal direction.

Figure 5 .Figure 5 .Figure 5 .

 555 Figure 5.30: Time-averaged void ratio (left) and velocity (right) profiles from station 1 to 5 -Kunz model comparison

C

  in this thesis deals with the study and implement of the three-equation and four-equation cavitation models which are developed in LEGI into the NSMB solver.

A

  comparison of various cavitation models coupled with turbulence models by OpenFOAM on 2D and 3D Venturi geometry was proposed. The interPhaseChangeFoam solver was used to simulate the cavitation pocket by the formulation of void ratio transport equation cavitation models including Kunz, Merkle and SchnerrSauer models. Numerical results have been compared with experimental data concerning the time-averaged void ration and longitudinal velocity, wall pressure, RMS wall pressure fluctuations and turbulent eddy viscosity. For the turbulence closure, three models are considered: the one-equation Spalart-Allmaras model, the two-equation CHAPTER 6. CONCLUSIONS AND PERSPECTIVES kε model and the Menter kω SST model. The Reboud eddy-viscosity limiter is introduced

  Based on the results obtained in this work, suggestions for improvements and future work are proposed as follows. Regarding the four-equation cavitation model implemented in the NSMB solver, the mass transfer rate ṁ in the formulation of the void ration transport equation could be introduced to better modeling the cavitation phenomenon. Then coupled to the turbulence model for turbulent cavitating flow. In the meanwhile, further looking into the code of NSMB is necessary by fresh viewpoints. For OpenFOAM simulation, an investigation of the calibration of the condensation and vaporization constants appearing in the mass transfer formulation to improve the results for the Venturi case. Regarding the options of the turbulence model, the nonlinear and hybrid models could be implemented for better prediction of the structure of characterized by the incompressible region of the pure liquid and compressible region of the two-phase mixture and associated with large variations in the local Mach number, where M < 0.1 in the liquid phase and M > 1 in the mixture zone.

  is the Jacobian matrix of the convective fluxes.Expressions of the matrices are derived as follow.

  have the matrices of R and R -1 as follows,

  Figure B.1: (a) An isometric view of the BSCW (b) Cross-sectional view of the SC(2)-0414 airfoil (c) BSCW model mounted in TDT

Figure B. 2 :Figure B. 5 :

 25 Figure B.2: Case 1 (Mach 0.7, R e = 1.12 × 10 7 , AoA = 3°): Mean C p for unforced system data at 60% and 95% wing span with the SA QCR 2013 turbulence model.

Table 2

 2

	.1.

Table 2

 2

.1: Class of models for cavitating flows

Table 3

 3

	.1: Directory structure of an OpenFOAM case
	case name
	system	controlDict
		decomposePar
		fvSchemes
		fvSolution
	constant	physical properties
		polyMesh	points
			cells
			faces
			boundary
	0	BC and initial conditions
	time directories

Table 4

 4 4. VALIDATION CASES three-equation model are in good agreement with the previous investigations[START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF][START_REF] Ansari | Numerical simulation of compressible two-phase flow using a diffuse interface method[END_REF] and there is no oscillation in the solution. However, the volume fraction of gas computed by the four-equation model shows oscillation at the outlet. This problem might be the issue of the boundary condition .

	.1: Properties of air and water and initial condition for interface movement in a uniform pressure
	and velocity flow.
		0 < x < 0.5		0.5 < x < 1
		Air	Water	Air	Water
	ρ (kg/m 3 ) 10	1000	10	1000
	α	10 -6	1 -10 -6	1 -10 -6	10 -6
	p ∞ (Pa)	0	6 × 10 8	0	6 × 10 8
	γ	1.4	4.4	1.4	4.4

The numerical solution is plotted in Figure

4

.1 and 4.2 at time t = 2.79 ms and is compared to the exact one. A mesh contains 200 uniform cells is used. The results obtained with the CHAPTER

Table 4 .

 4 2: Properties of air and water and initial condition for the water-air shock tube.

		0 < x < 0.7		0.7 < x < 1
		Air	Water	Air	Water
	ρ (kg/m 3 ) 1	1000	1	1000
	P (Pa)	10 9	10 9	10 5	10 5
	α	10 -6	1 -10 -6	1 -10 -6	10 -6
	p ∞ (Pa)	0	6 × 10 8	0	6 × 10 8
	γ	1.4	4.4	1.4	4.4

Table 4 .

 4 3: Parameters of the stiffened gas EOS for water at T = 355 K.

	γ	P ∞ (Pa)	q (J/kg)	C p (J/K kg) ρ sat (kg/m 3 )
	Liquid 2.35	10 9	-0.1167 × 10 7	4267	1149.9
	Vapor 1.43	0	0.2030 × 10 7	1487	0.31

Table 4 .

 4 4: Properties of air and water and initial condition for the water-air shock tube.

		γ	ρ (kg/m 3 ) u (m/s) v (m/s) P (Pa)	P ∞ (Pa)
	Water	4.4 1000	0	0	1 ×10 5	6 ×10 8
	Air	1.4 1	0	0	1 ×10 5	0
	post-shock 4.4 1323.65	681.58 0	1.9 ×10 9 6 ×10 8

Table 5

 5

	.1.

is applied at the upstream inflow and a pressure outlet condition is used at the outlet boundary for the computational set-up. The vaporization pressure P vap is set to 2300 Pa. The time step ∆t and the maximum Courant number are set to 10 -5 s and 1.0 respectively.

Table 5 .

 5 2: Boundary conditions, flow and turbulence properties of the Venturi tested cases

Table 5

 5

	Cavitation model
	Kunz	U ∞ = 10.8(m/s)	t ∞ = 0.023(s)

.3. Here, U ∞ is set to the freestream value, t ∞ represents a relaxation time, n is the bubble number density, dnuc is the nucleation site diameter and C c and C v are the condensation rate coefficient and vapourisation rate coefficient respectively.

Case 1, BSCW Mean Cp Upper surface at 95% Wing Span

											B.3. COMPUTATIONAL RESULTS
				AePW-2 Case 1, BSCW Mean Cp					AePW-2 Case 1, BSCW Mean Cp
	-1.5		Upper surface at 60% Wing Span		-1.5		Lower surface at 60% Wing Span
								Coarse Grid								Coarse Grid
								Medium Grid								Medium Grid
		-1						Fine Grid		-1						Fine Grid
								Experiment								Experiment
	-0.5							-0.5
	CP	0							CP	0
		0.5								0.5
		1								1
		-0.2 1.5	0	0.2	0.4	0.6	0.8	1	1.2	-0.2 1.5	0	0.2	0.4	0.6	0.8	1	1.2
					X/C								X/C
				AePW-2 Case 1, BSCW Mean Cp					AePW-2
	-1.5		Upper surface at 95% Wing Span		-1.5
								Coarse Grid								Coarse Grid
								Medium Grid								Medium Grid
		-1						Fine Grid		-1						Fine Grid
	-0.5							-0.5
	CP	0							CP	0
		0.5								0.5
		1								1
		-0.2 1.5	0	0.2	0.4	0.6	0.8	1	1.2	-0.2 1.5	0	0.2	0.4	0.6	0.8	1	1.2
					X/C								X/C

2.1. MODELING OF TWO-PHASE FLOWS

4.3. WATER-AIR MIXTURE EXPANSION TUBE

ACKNOWLEDGEMENTS

F

irst of all, I would like to express my deepest gratitude and special thanks to my supervisor, Prof. Yannick Hoarau, who took time out to hear, guide and keep me on the correct path

symmetry. When this water jet impacts the stationary water at the front of the bubble (time t = 2 µs), an intense blast wave also called water hammer shock [START_REF] Hawker | Interaction of a strong shockwave with a gas bubble in a liquid medium: a numerical study[END_REF] is formed generating a high-pressure zone. The blast front, which expands continuously, is highly asymmetric due to the high-speed water jet (time t = 2.4 µs). The rightward blast wave increases as a spherical wave. Both shocks lose strength as they advance, the rightward wave more so than its leftward twin. The interaction of the blast wave with the bubble fragments lead to high pressure levels (time t = 2.7 µs). Finally, at time t = 2.9 µs, the blast wave continues its expansion and the cavity its shrinkage. These results show a good agreement with previous numerical results [START_REF] Ball | Shock-induced collapse of a cylindrical air cavity in water: a free-lagrange simulation[END_REF][START_REF] Nourgaliev | Adaptive characteristics-based matching for compressible multifluid dynamics[END_REF][START_REF] Ozlem | A numerical study of shock-induced cavity collapse[END_REF].

The pressure evolution during the collapse is illustrated in Figure 4.9. During the impact of the water jet with the stationary water at the front of the bubble, a blast wave is generated leading to the pressure increase (time t = 2 µs). As previously described, the blast fronts are highly asymmetric. The rightward wave increases as a spherical wave and expands continuously in the radial direction ((time t = 2.4 µs and after). The shock intensity decreases during the propagation, especially for the rightward front. At time t = 2.4 µs, the more intense pressure peak is generated by the leftward front on the bubble axis. At time t = 2.7 µs, the interaction of the leftward blast wave with the bubble pieces leads to a very strong pressure peak, which is the most intense reached during the collapse. At time t = 2.9 µs, the low-pressure area inside the vortices core are well illustrated. Both blast wave fronts continue to expand.

The time evolution of the axial velocity is plotted in Figure 4.10. The reflected rarefaction wave, resulting from the impact of the incident shock with the upstream bubble interaction, relaxes the pressure, which accelerates the flow and forms a high-speed water jet (time t = 1.1 µs and after). The velocity magnitude is higher than 2000 m/s. At time t = 2 µs, the water jet strikes the downstream bubble interface leading to the blast wave generation. After time t = 2.1 µs, the bubble is cut in half and forms a pair of distinct vortical structures. The developing leftward wave advances relatively slowly due to the high water velocity in the jet fluid. After time t = 2.6 µs, the front of the leftward wave can be observed (abscissa x 0.012 m).

The evolution of the vertical velocity is illustrated in Figure 4.11 at the same time as previously. When the water jet impinges with the downstream edge of the bubble, at time t = 2 µs, the bubble forms a pair of distinct vortical structures. At time t = 2.6 µs, caused by the leftward blast wave, secondary jets penetrate into the smaller bubbles and cut the initial bubble into four pieces.

CHAPTER 5. RESULTS ON THE VENTURI GEOMETRY CAVITATING FLOW peak position varies among models. With the KWReboud turbulence model, the peak is present downstream station 5, whereas the peak obtained with the SAReboud turbulence model is upstream. The KEReboud turbulence model provides fluctuations in better agreement with experimental data inside the cavity but not in the re-compression area. Figure 5.12 compares time-averaged profiles of the viscosity ratio µ t /µ, at the five station, obtained with Merkle cavitation model. At stations 2 and 3, a drastic decrease of µ t close to the wall for all models due to the Reboud correction can be observed. Therefore, a better prediction of the unsteadiness and the separation can be expected. The KWReboud turbulence model shows a large reduction of the ratio in the sheet at all stations. On the contrary, the SAReboud turbulence is not in the same case.

Entropy

and

B.1 Introduction

The aeroelastic prediction workshop aims at the assessment of the state-of-the-art in numerical methods for simulating flow-fields about wings undergoing prescribed motions or static and dynamic aeroelastic deformations at transonic flight conditions. The first AIAA Aeroelastic Prediction Workshop (AePW-1) conducted three configurations: the Rectangular Supercritical Wing (RSW), the Benchmark Supercritical Wing (BSCW) and the High Reynolds Number Aerostructural Dynamics (HIRENASD). These cases focus on the prediction of unsteady pressure distributions resulting from forced motion. The second workshop (AePW-2) extends the benchmarking effort to aeroelastic flutter solutions with flow conditions in transonic regime and focuses on a single configuration. The configuration chosen for the second workshop is the BSCW. The primary analysis condition has been chosen such that the influence of separated flow is considered to be minimal, yet a shock is still present. This is a step back in flow complexity from the BSCW cases for AePW-1.

The goal in moving to the lower transonic Mach number is to have analysis teams progress through unforced system analyses, forced oscillation solutions and flutter analyses. Revisiting the AePW-1 analysis condition is included in AePW-2 as an optional case, also extending it to include flutter solutions.

B.2 The Benchmark Supercritical Wing

The BSCW is chosen for the workshop. It is a rigid, semispan, rectangular supercritical wing with a chord of 16 inches, a span of 32 inches, and a SC(2)0414 supersonic airfoil with design normal force coefficient of 0. Tunnel (TDT). The first experiment in 1991 is a flutter test performed on the pitch and plunge apparatus (PAPA) system which provides two-degree-of-freedom flutter [START_REF] Dansberry | Experimental unsteady pressures at flutter on the supercritical wing benchmark model[END_REF].

The PAPA experimental data consists of unsteady data at flutter points and averaged data on a rigidified apparatus at the flutter condition. The second experiment in 2000 is a forced excitation test performed on the oscillating turntable (OTT), in which the wing was oscillated in pitch about an axis at the 30% chord [START_REF] Heeg | Experimental data from the benchmark supercritical wing wind tunnel test on an oscillating turntable[END_REF] 60% and 95% span stations for the unforced system by considering the effects of grid size for each selected turbulence models. The unforced system computations show good agreement with the experimental data, except that the peak value at the 10% wing span near the leading edge on the upper surface is under-predicted. In addition, form these figures the influence of the grid sizes is minor.

B.3.2 Test Case 2

The second test case focuses on flutter prediction at Mach 0.74 and 0°angle of attack. The experimental comparison results obtained from the PAPA test are available for both the inboard span station (60% wing span) and the outboard span station (95% wing span).

Comparison of the steady-state solution using different turbulence models among different grid sizes for the mean pressure coefficients at 60% and 95% span stations are shown from Figure B.9 to B.11. For the 60% span station all computational results for the unforced system converged to the same results and are in good agreement to the experimental data for both the upper and lower surface. For the 95% span station the computations show good agreement with the experimental data, except that the peak value at the 10% wing span near the leading edge on the upper surface is over-predicted. From these two sets of comparison, the computational results are in very good agreement with those of experiment. The influence of turbulence models and grid sizes is not evident for test case #2.

B.3. COMPUTATIONAL RESULTS

B.3.3 Test Case 3

This test case (Mach = 0.85, 5°angle of attack) was also analyzed in AePW-1 and the results from the different analysis teams showed that it is difficult to find a converged steady state solution. This case is chosen as an optional case for AePW-2 in order to continue working on the problem and resolving the discrepancies. For both the unforced system and forced oscillation, experimental data of mean pressure coefficients obtained from the OTT test are available at the inboard span station (60% wing span).

B.4 Conclusion

Numerical simulations were performed with the NSMB solver for the second Aeroelastic Prediction Workshop (AePW-2). For case#1, at Mach 0.7, 3°angle of attack, the effects of grid sizes and turbulence model are not obvious. The computation results are in very good agreement with the experimental data. There exists only some difference in the suction peak region in the upper surface and in some regions of the lower surface. For forced oscillation simulation, computational results are performed with the medium grid. All turbulence models predict well the pressure values after the shock region.

For case #2, at Mach 0.74, 0°angle of attack, the simulations predicted accurately the steady pressure distribution. The influnece of turbulence models and grid size is not obvious.

For case #3, at Mach 0.85, 5°angle of attack, the effect of grid sizes is not important for the computational results of the mean pressure distribution. The kω SST turbulence model captures better the upper-and lower-surface shock location, but fails to predict the pressure values at the recovery area behind the shock. For forced oscillation simulation, all turbulence models are not able to predict the mean pressure for the upper surface.