FLUIDES CRITIQUES, EFFET PISTON, MODÉLISATION Ecoulement Compressible

Champ Numérique

De Phase

Ms Mridula

Sharma " Vi

Prof Sakir

Jean Paul Caltagirone

Prof K Tapan

Iit Sengupta

Kanpur

Nacho

Elliott

Lucas

Dr Justine Lin

Shaheen, Swati, Vibhor, Suraj, Binu, Arvind, Cheshtaa, Sirisha Shivani Anandi

Keywords: COMPRESSIBLE FLOW, NEAR CRITICAL FLUID, PISTON EFFECT, NUMERICAL MODELLING COMPRESSIBLE FLOW, NEAR CRITICAL FLUID, PISTON EFFECT, NUMERICAL MODELLING

A fluid, in addition to its liquid and gas phase, is known to exist in another phase, wherein the fluid inherits some properties of both the phases. Such a fluid is called a supercritical fluid and the conditions (pressure and temperature) beyond which the fluid exists in this state is called the critical point. One of the peculiar feature of the fluids near the critical point is that the various thermo-physical properties show a singular behavior, such as diverging compressibility, vanishing thermal diffusivity etc. The flow behavior near the critical point leads to intriguing flow features ascribed to the strong thermo-mechanical coupling whose in-depth investigation can be limited by experimental constraints especially during a continuous transition from supercritical to subcritical regime. The current work focuses on analyzing the flow behavior in near-critical fluids with prime focus on supercritical fluids. This is achieved by developing a mathematical and numerical model which is followed by the validation study and error analysis of the numerical scheme wherein unusual behavior of the Courant number is observed. Subsequently, the flow behavior of supercritical fluid is studied when simultaneously subjected to thermal quench and vibration, mainly Rayleigh-vibrational and parametric instabilities, their physical mechanism and various parameters affecting them. In addition, two captivating phenomena, firstly where the temperature of the fluid region drops below the imposed boundary condition and secondly, the see-saw motion of the thermal boundary layer are observed and physical explanations are provided. In order to investigate the flow dynamics in subcritical regime, phase-field modelling approach is explored for isothermal conditions. The model is examined for elementary test cases illustrating the feasibility to extend the model for a continuous transition from supercritical to subcritical regime.

Les simulations numériques nous ont permis de relever deux phénomènes particulièrement surprenants : (i) la température du fluide à l'intérieur du domaine devient inférieure à la trempe de température imposée à la frontière et (ii) une oscillation des doigts d'instabilité apparaît dans la couche limite thermique dans la direction de la vibration. Dans le cas des fluides sous le point critique (cas diphasique), le modèle compressible développé est couplé à un de champ de phase ("phase field") dans les conditions isothermes. Des cas tests élémentaires ont été considérés avec succès. Une discussion est proposée afin d'étendre le modèle dans le cas d'une transition continue du régime supercritique au régime sous-critique et vice-versa.

x RÉSUMÉ Un fluide est connu généralement pour exister principalement sous forme de liquide ou gaz.

Cependant, au-delà du point critique thermodynamique, défini en termes de pression et de température du fluide, une distinction claire entre les phases liquide et vapeur cesse d'exister. Le fluide est alors appelé fluide supercritique. Lorsque l'on s'approche du point critique les propriétés thermophysiques deviennent singulières, par exemple la compressibilité diverge et la diffusivité thermique s'annule, ce qui rend les fluides proches du point critique particulièrement intéressants dans de nombreuses applications allant de l'espace à la micro technologie. D'un point de vue plus fondamental, le comportement singulier de ces propriétés conduit à des caractéristiques d'écoulement surprenantes et inattendues, principalement attribuées au couplage thermomécanique du champ d'écoulement. Alors que les contraintes expérimentales peuvent limiter une analyse en profondeur de ces systèmes, un modèle mathématique approprié peut aider à mieux comprendre les caractéristiques d'écoulement dans les fluides proches du point critique, notamment la transition continue de l'état sous critique à l'état supercritique et vice-versa dans l'hypothèse d'un milieu continu. Les travaux de recherche dans cette thèse marquent un premier pas vers cet objectif dans lequel un modèle mathématique compressible est d'abord développé pour analyser les écoulements dans les fluides supercritiques suivi d'une modélisation par champ de phase du cas diphasique sous le point critique. Le modèle développé calcule la masse volumique directement à partir de l'équation de la conservation de la masse sans la nécessité d'un algorithme itératif de couplage vitesse-pression. Le modèle est résolu numériquement par la méthode des volumes finis avec le code du laboratoire, Thetis. La validation du modèle a été faite dans le cas de la propagation de l'onde acoustique dans l'eau. De plus, on a observé que des solutions précises avec des schémas implicites pour des systèmes non linéaires sont possibles avec des nombres de Courant élevés.

Ceci s'explique par les contributions à l'erreur provenant non seulement du taux de croissance des erreurs conventionnelles mais aussi du retard de phase et de la différence des vitesses de groupe des ondes numériques et physiques. Le modèle est ensuite utilisé pour étudier les phénomènes d'écoulement dans des fluides supercritiques lorsqu'ils sont assujettis à une trempe thermique et à une vibration simultanées qui mènent aux instabilités thermo-vibrationnelles, en particulier celles de Rayleigh-vibrationnelles et paramétriques. La masse volumique étant évaluée directement à xi partir de l'équation de conservation de la masse, le modèle n'est pas alors limité par des petites perturbations thermiques imposées par l'hypothèse de l'équation d'état linéarisée dans la littérature.

Les résultats des simulations numériques et des expériences concordent bien. Une description détaillée des mécanismes physiques à l'origine de ces instabilités est menée par l'analyse de l'effet des différents paramètres (telles que l'amplitude de la trempe de température, la proximité du point critique, l'amplitude et la fréquence de la vibration) sur l'amplitude critique de vibration pour le démarrage de ces instabilités et le comportement des motifs ondulatoires observés dans les couches limites thermiques. Un phénomène intéressant a été observé : la température du fluide à l'intérieur de la cellule complètement fermée descend en dessous de la température imposée à la limite. Cette anomalie a été expliquée en détail pour différentes conditions aux limites pour la température. Par ailleurs, les doigts d'instabilité apparaissant proche de la paroi subissant la trempe de température, ont montré un mouvement de va-et-vient dans la direction de vibration et ceci a été expliqué en termes d'effet piston forcé (vibration) par analogie à l'effet piston classique dû au changement de pression induit par le changement de température. L'étude du cas diphasique (sous le point critique) par une approche par champ de phase associée au modèle compressible décrit précédemment, a été considérée. Des cas tests élémentaires sous le point critique et relativement proche du point critique ont été considérés. Une discussion de l'extension de ce modèle pour différentes applications dans le domaine diphasique sous-critique est également proposée.

CHAPTER 1. INTRODUCTION

A fluid in general is known to exist primarily in two phases, either a liquid or a gas.

However, a French scientist, Baron Charles Cagniard de Latour [START_REF] De | Experimente unter hohem drucke[END_REF], while analyzing a phase equilibrium between equal volumes of liquid and gaseous CO2 in a sealed gun, found that the expected splashing noise of a liquid phase disappeared above certain conditions of temperature.

This led to the discovery of a new phase of a fluid, wherein the phase boundary between liquid and gas phase ceases to exist, and the fluid is called a supercritical fluid (SCF). Following this discovery, the subject of supercritical or near-critical fluids (fluid near the critical point) has piqued interests of several researchers and till date ascribed to experimental and numerical challenges, several expects remain unexplored.

In the quest to further understand the fluid behavior near the critical point, the current thesis focuses on the study of near critical fluids, primarily supercritical fluid by developing compressible mathematical and numerical models.

In this chapter, a brief introduction to the near-critical fluids is presented followed by the description of the thesis objectives and outline.

Fluids near the critical point

The fluid flow and thermal behavior of a fluid is governed by its various thermo-physical properties which depend on whether the fluid is in solid, liquid or gas phase. The existence of a fluid in either of this phase is governed by pressure and temperature conditions (𝑃, 𝑇). In each phase, the fluid possesses different properties which make these phases distinguishable from each other. The existence of a fluid in different phases as a function of pressure, temperature and volume is expressed using a phase-diagram as shown in Fig. 1-1 for any pure substance (a substance which is spatially uniform in chemical composition and homogenous). The projection of this 3D illustration on 𝑃 -𝑇 plane provides further insight into the existence of various phases under different conditions, as marked in the Fig. 1-1(b). Of particular interests are the curves known as the liquid-vapor co-existence curves, which represent the loci of (𝑃, 𝑇) such that the two phases are in equilibrium with each other, unlike other regions wherein the fluid is known to exist only in a single-phase. An interesting behavior can be observed on moving along the liquid-vapor co-existence curve in 𝑃 -𝑇 plane where the co-existence curve seems to end at a temperature denoted by 𝑇 𝑐 i.e. there is phase transition from two-phase region to singlephase region. This point is called the critical point of a fluid (denoted by 𝑃 𝑐 , 𝑇 𝑐 ) beyond which a clear demarcation between the liquid and the gas phase ceases to exist and a single-phase fluid is attained. Instead, the fluid behaves as an intermediate between the liquid and the gas phase exhibiting some properties of the liquid, such as high density and solubility like a liquid and low viscosity and high compressibility like a gas. The fluid beyond this critical point is termed as a supercritical fluid.

More recently, a striking observation was reported by Simeoni et al. [START_REF] Simeoni | The Widom line as the crossover between liquid-like and gas-like behaviour in supercritical fluids[END_REF] and Gorelli et al.

[3] wherein they have shown that the homogeneous regime considered above the critical point can be divided into liquid-like and gas-like regions, similar to the subcritical state. Using inelastic Xray scattering to study the dispersion, 𝑖. 𝑒. the dependence of the sound-speed on the frequency (a phenomenon peculiar to liquid and not gases), they were able to attribute this property to supercritical fluids. The line dividing the two phases in the supercritical regime, was described as an extension of the co-existence curve below the critical point whose properties were investigated by Sciortino et al. [START_REF] Sciortino | Line of compressibility maxima in the phase diagram of supercooled water[END_REF] using molecular dynamics computations and this line was later termed as the 'Widom' line (as shown in Fig. 12). A similar line, defined as percolation line, has also been identified by Campi et al. [START_REF] Campi | Percolation line of self-bound clusters in supercritical fluids[END_REF] and was later investigated in detail by Skor and Nezbeda [START_REF] Škvor | Percolation line and response functions in simple supercritical fluids[END_REF]. 𝜌 𝑐 is the critical density, 𝐶 𝑃,𝐼𝐺 is the ideal gas specific heat at constant pressure, R is the gas constant, and γ the isentropic exponent (reproduced from [START_REF] Banuti | Similarity law for Widom lines and coexistence lines[END_REF] ).

Banuti [START_REF] Banuti | Crossing the Widom-line -Supercritical pseudo-boiling[END_REF] investigated the same problem by mixed analytical and empirical approaches and showed that unlike the phase transition in subcritical regime, the phase transition in supercritical state occurs over a finite temperature interval and thus categorized it as a non-equilibrium process.

In addition to the energy required to overcome the usual intermolecular attraction, the transition in supercritical state requires an additional energy to increase the temperature. Banuti et al. [START_REF] Banuti | Similarity law for Widom lines and coexistence lines[END_REF] have further proposed a Widom line functional based on the Clapeyron equation. In their recent article, Banuti et al. [START_REF] Banuti | Seven questions about supercritical fluids -towards a new fluid state diagram[END_REF] summarized the fundamental aspects in supercritical regime, highlighting the relevance of liquid-like supercritical and gas-like supercritical fluid in combustion and injection processes. However, in the current work, the supercritical region has been considered as a singlephase region.

Looking at these unusual behaviors, a very important question then arises, how does this transition affect the thermo-physical properties?

Behavior of the properties and critical exponents

The behavior of various thermo-physical properties can be analyzed if the relation between various thermodynamic variables is known, for example using an equation of state defining the relation between pressure, density and temperature. Intuitively, such an equation can also be used to analyze the phase transition and was first described using Van der Waals (vdW) equation of state near the critical point of a fluid. An elementary analysis using this equation of state shows an unusual behavior of various properties such as vanishing density difference, diverging isothermal compressibility etc. However, very close to the critical point, this has been found to be incapable to capture an accurate behavior of the thermo-physical properties [START_REF] Zappoli | Heat Transfers and Related Effects in Supercritical Fluids[END_REF]. This is primarily attributed to the behavior of fluctuations on approaching the critical point. It is known that locally, the various properties exhibit fluctuations from their mean values, which are more statistical in nature. With majority of the macroscopic thermo-physical properties being related to the statistical distribution of these fluctuations, it is evident that these fluctuations will have a significant effect on various properties. It has been established that on approaching the critical point, the behavior of macroscopic properties is no longer ascribed to the local behavior of molecules but due to the collective behavior arising from the long-range interactions. These long-range interactions, which exist along short-range interactions (defining the background contribution) thus govern the behavior of various properties. Consequently, the specificity of a system gets hidden leading to the universal behavior which has been termed as critical universality. This implies, that on approaching the critical point, the various properties will behave in a similar manner independently of the fluid being considered, and thus have been defined based on a single parameter with critical exponent of each property [START_REF] Zappoli | Heat Transfers and Related Effects in Supercritical Fluids[END_REF][START_REF] Stanley | Introduction to Phase Transitions and Critical Point Phenomena[END_REF] . These critical exponents have been derived by renormalizing the fluctuations (Renormalization Group Theory) and are described in Table 1-1 [START_REF] Zappoli | Heat Transfers and Related Effects in Supercritical Fluids[END_REF] . Thus, by virtue of the dominance of these long-range interactions, the various thermodynamic properties tend to show a singular behavior, 𝑖. 𝑒. some properties such as, isothermal compressibility(𝜒 𝑇 ), thermal expansion coefficient ( 𝛽 𝑃 ), thermal conductivity (𝑘), specific heat at constant volume and constant pressure (𝐶 𝑉 , 𝐶 𝑃 ) diverge while thermal diffusivity (𝐷 𝑇 ) and sound speed (𝑐 𝑠 ) tend to zero as can be observed from the behavior of critical exponents in Table 1-1. Thermal conductivity (𝑘) 𝜀 -𝜑 𝜑 = 0.58 kinematic viscosity (𝜐) 𝜀 -𝜂 𝜂 = 0.04

Correlation length 𝜀 -𝜈 𝜈 = 0.63

The general behavior of these thermodynamic properties is obtained by fitting experimental results to a law with a constant background term (corresponding to the short-range interactions or mean contribution) and a temperature dependent critical contribution term as a function of reduced proximity to the critical point, 𝜀 = 𝑇-𝑇 𝑐 𝑇 𝑐

. These two contributions are generally separated from each other by assuming an additive relation, wherein the critical part is added to the regular part [START_REF] Zappoli | Heat Transfers and Related Effects in Supercritical Fluids[END_REF].

For any given quantity, 𝑌, this can be expressed as, 𝑌(𝜀) = 𝑌 𝑐 (𝜀) + 𝑌 𝑟𝑒𝑔 (𝜀), where 𝑌 𝑐 (𝜀) and 𝑌 𝑟𝑒𝑔 (𝜀) denote the critical and regular contribution. The background is generally expressed by a polynomial function, as [START_REF] Zappoli | Heat Transfers and Related Effects in Supercritical Fluids[END_REF],

𝑌 𝑟𝑒𝑔 (𝜀) = 𝐴 0 + 𝐴 1 𝑇 + 𝐴 2 𝑇 2 + . . . . . .. (1.2.1)

while the critical contribution is represented by non-analytical functions as a power law multiplied by crossover functions (to account for deviations away from the critical point) as [START_REF] Zappoli | Heat Transfers and Related Effects in Supercritical Fluids[END_REF],

𝑌 𝑐 (𝜀) = 𝜓 𝑌0 𝜀 ±𝜓 (1 + 𝑎 𝑌 (1) 𝜀 𝛥 + 𝑎 𝑌 (2) 𝜀 2𝛥 + ⋯ ) (1.2.2)

Here, 𝜓 𝑌0 is the critical amplitude, while 𝜓 is a universal exponent and (1 + 𝑎 𝑌 (1) 𝜀 𝛥 + 𝑎 𝑌 (2) 𝜀 2𝛥 + ⋯ ) denotes the crossover function with 𝛥(≈ 0.502) representing the universal exponent while 𝑎 𝑌 (𝑖)

represent coefficients which depend on the type of fluid. However, in most of the cases, the behavior can be well defined by neglecting the higher order terms in cross-over functions and the properties of a fluid in general can be represented by [START_REF] Zappoli | Heat Transfers and Related Effects in Supercritical Fluids[END_REF],

𝑌(𝜀) = 𝐴 0 + 𝜓 𝑌0 𝜀 ±𝜓 (1.2.

3)

It can therefore be deduced that on approaching the critical point, the contribution from the critical term becomes predominant. When these distinctive properties of the fluid become distinguishable from mean contributions for certain values of pressure and temperature, the fluid is termed as near-critical fluid.

Piston Effect: Heat transfer mechanism in near-critical fluids

The heat transfer in any substance or system is primarily governed by its thermal diffusivity, which reflects how fast the thermal perturbation will homogenize the system. On approaching the critical point, the thermal diffusivity tends to vanish and thus, intuitively it was believed that this will have as an impeding effect on the thermal homogenization of a supercritical (near-critical) fluid when subjected to thermal heating. However, Nitsche and Straub [START_REF] Nitsche | The critical''HUMP"of Cv under microgravity: results from D1-Spacelab experiment"Wärmekapazität[END_REF] in their experimental studies with supercritical SF6 observed that the bulk temperature followed the wall temperature very fast thereby causing thermal relaxation in seconds as opposed to the predicted scale of days based on thermal diffusion. This phenomenon was explained in 1990s by three independent teams (Onuki et al. [START_REF] Onuki | Adiabatic heating effect near the gas-liquid critical point[END_REF] , Boukari et al. [START_REF] Boukari | Critical speeding up in pure fluids[END_REF], and Zappoli et al. [START_REF] Zappoli | Anomalous heat transport by the piston effect in supercritical fluids under zero gravity[END_REF]) and was termed as piston effect and can be explained as follows.

The phenomenon can be understood by considering a cell filled entirely with a supercritical fluid subjected to a heat flux at one of its boundary. A very thin thermal boundary (TBL) layer is formed due to the vanishing thermal diffusivity. The fluid in the TBL expands due to its high thermal expansion at constant pressure (𝛽 𝑃 ) causing a propagation of an acoustic wave in the bulk which compresses the bulk due to its very high isothermal compressibility (𝜒 𝑇 ). With each passing compression wave, the fluid converts some of the kinetic energy into thermal energy, heating the bulk adiabatically. Further, the compression wave is associated with change in pressure and thus of the same nature as acoustic wave, which results in thermal homogenization to occur on very short time scales termed as piston effect time scale. Similar argument can be applied when the boundary wall is cooled wherein the fluid in the TBL contracts and adiabatic cooling thermalizes the bulk fluid. Since this resembles to the action of a piston which compresses a gas in a closed cylinder, this phenomenon has been termed as the piston effect and forms the basis of understanding of most of the phenomena in supercritical (near-critical) fluids.

It was shown by Onuki et al. [START_REF] Onuki | Fast adiabatic equilibration in a single-component fluid near the liquid-vapor critical point[END_REF] that the piston effect can be present in any compressible fluid which can be summarized as follows. The temperature evolution of a fluid, in the absence of any convection can be described by the classical heat diffusion equation. However, for a compressible fluid, it is necessary to add the effect of thermal expansion in the form of work done by the pressure force. This work done varies with the average temperature of the fluid and the ratio of heat capacities, (𝛾 = 𝐶 𝑃

𝐶 𝑉

). What makes this effect dominant and observable in supercritical (nearcritical) fluids is the diverging behavior of 𝛾 (𝐶 𝑃 diverges with a critical exponent much higher than 𝐶 𝑉 , see Table 1-1). Thus, even though the piston effect can be considered ubiquitous in any compressible medium, it becomes significant only near the critical point where the effects of diffusion are subdued while those of compressibility increase rapidly. The time scale corresponding to the homogenization caused by the piston effect is defined as a piston effect time scale and is defined by [START_REF] Onuki | Fast adiabatic equilibration in a single-component fluid near the liquid-vapor critical point[END_REF],

𝑡 𝑃𝐸 ≈ 𝑡 𝐷 (𝛾-1) 2 (1.3.1)
where 𝑡 𝐷 represents the diffusion time scale (

𝐿 2 𝐷 𝑇
), 𝐿 being the characteristic length scale.

Table 1-2 compares the various time scales for hydrogen with varying proximities to the critical point.

Applications of near-critical fluids

The need of scientific research has been largely driven by challenging applications in various fields, varying from space to micro technology. Over the past several decades, there has been a

growing demand for supercritical fluids in various industrial and scientific applications such as varied as alternative eco-friendly refrigerants [START_REF] Lorentzen | A new, efficient and environmentally benign system for car airconditioning[END_REF][START_REF] Kim | Fundamental process and system design issues in CO2 vapor compression systems[END_REF], chemical extraction/separation processes, supercritical chromatography, drying and catalysis [START_REF] Knez | Industrial applications of supercritical fluids: A review[END_REF]. In addition, these also find their applications in particle formation ranging from nano to macro dimensions, supercritical water oxidation for destruction of aqueous based organic waste [START_REF] Marrone | Supercritical water oxidation-Current status of full-scale commercial activity for waste destruction[END_REF] , electronic chip manufacturing [START_REF] Weibel | An overview of supercritical CO2 applications in microelectronics processing[END_REF],

drilling technologies [START_REF] Rothenfluh | Penetration length studies of supercritical water jets submerged in a subcritical water environment using a novel optical Schlieren method[END_REF][START_REF] Augustine | Hydrothermal flames: From phenomenological experimental demonstrations to quantitative understanding[END_REF] and rocket fuels [START_REF] Pizzarelli | CFD analysis of transcritical methane in rocket engine cooling channels[END_REF][START_REF] Neill | Practical uses of liquid methane in rocket engine applications[END_REF] Table 1-2: Typical physical timescales of dynamic processes in supercritical 𝐻 2 along its critical isochore confined in a one-dimensional slab of 10 cm In space applications, the storage of cryogenic fluids (Oxygen, Hydrogen) is done under supercritical conditions as the management of these fluids is a real problem due to the uncertain localization of gas and liquid phases. Among these several applications, let us look at some applications pertinent to supercritical (near-critical) fluid and see what additional attributes lie in these.

ε(= T-T c T c ) x 10

➢

In rocket engines, methane is being considered as a denser and cheaper replacement of hydrogen in launch vehicles. The interest in such flows for liquid rocket engines is driven by the fact that the thrust chamber is cooled by one of the available propellants, which flows in a suitable narrow channel. A peculiar behavior is observed if methane is used as coolant: it enters channels at supercritical pressure and subcritical temperature and then, under heating from the hot gas, its temperature increases and can reach the supercritical regime. A further interesting behavior can be ascribed when just at the start of rocket engine, the pressure will not be high and thus the methane can transition from supercritical to sub-critical regime.

➢

In CO2 air to-water heat pump or in supercritical CO2 based refrigeration systems, the evaporation is performed at low temperatures in the subcritical regime whereas the condensation occurs at the supercritical state. The overall efficiency is increased which is beneficiary in a certain number of applications such as the production of domestic hot water or industrial drying processes. This is attributed to the fact that heat rejection in supercritical state is primarily a single-phase phenomenon unlike conventional system wherein phase transition reduces the thermodynamic efficiency. An important path that completes the refrigeration/heating cycle can be identified with transition from supercritical to sub-critical regime.

➢

In super-critical water reactors, supercritical water oxidation (SCWO) (a high-pressure and a high-temperature process for the destruction of toxic, hazardous or non-bio-degradable aqueous organic waste) is carried out above the critical point of water. Corrosion of reaction vessels and plugging due to salt precipitation are one of the major challenges of the process. In order to prevent these issues, it is desirable to have the inlet flow entering the reactor under subcritical conditions and thus the analysis of correct physics involves transition from subcritical to supercritical regime.

➢

Recently, hydrothermal spallation drilling has been considered as an alternative for deep drilling, depths > 10 𝑘𝑚, for geothermal applications. In order to induce thermal spallation, the necessary high temperatures and heat fluxes can be achieved using hydrothermal flame technologies. The hydrothermal flames in such conditions are usually formed using water which at such depth exceeds its critical conditions and unlike in its liquid state, supercritical water is a good solvent for compounds such as oxygen, nitrogen or carbon dioxide. Further, in the absence of any mass transfer limiting phase boundaries, combustion reactions between organic fuels (e.g. methane)

and oxidizing agents (e.g. oxygen, air) are easily facilitated in supercritical water. Therefore, while initial conditions on the ground pertain to the sub-critical state, below the ground water is in supercritical state and hence a transition from one state to other.

It can be seen that while in some of the applications, only the supercritical pressures may be attained with sub-critical temperature, in others there is a continuous transition from single phase to two phase regions.

Motivation and Scope of the current work

With the description of piston effect, explaining the speeding up of thermal homogenization, it can be well ascertained that the singular behavior of the thermo-physical properties can cause various intriguing phenomena. While experimental constraints can limit an in-depth insight into the flow behavior due to complex conditions near the critical point, numerical modelling can thus prove to be a practical and flexible tool in understanding the fundamentals of flow and transport characteristics in this region, especially a continuous transition from supercritical to sub-critical state and vice-versa in the realm of continuum mechanics. The current thesis thus marks a first step towards this goal and the important objectives of the thesis are three folds,

• First, to develop a mathematical and numerical model which can capture the physics of highly compressible supercritical fluid. The importance of such an objective is motivated by the use of a linear state equation in the existing models which circumvents the analysis to a higher thermal heating/cooling and the closest approach to the critical point.

• Secondly, to perform numerical investigations using the developed model in order to gain insights into the various fundamental mechanisms that arise in supercritical fluids in weightlessness conditions when subjected to simultaneous thermal quench and longitudinal vibrations. While this has been primarily motivated by applications in space technology, the absence of gravity evades the effect of convection which may otherwise subdue the effect arising solely from the thermo-mechanical coupling in supercritical fluid.

• Lastly, to explore the possibility to extend this model to analyze the flow behavior in subcritical state which may then be extended to study a continuous transition. One of the primary challenge attributed to this, unlike usual two-phase liquid-vapor system wherein the liquid phase is considered to be incompressible, both the phases, viz. liquid and gas, are highly compressible.

In order to reach these objectives, the current work is primarily oriented towards the Chapter 4 presents the results pertaining to the investigation of thermo-vibrational instabilities, mainly Rayleigh-vibrational and parametric. The results are validated with experimental observations followed by the description of the physical mechanism causing these instabilities. In addition, the effect of various parameters on the critical amplitude for onset of these instabilities, effect on wavelength and a stability plot is described.

Chapter 5 highlights two intriguing phenomena which were observed in conjunction with Rayleighvibrational instabilities. The first one being the drop of fluid temperature below the imposed temperature at the boundary whilst the second one explains the observed see-saw motion of the thermal boundary layer. Both these have been ascribed to the high compressibility leading to strong thermo-mechanical coupling in supercritical fluids.

Chapter 6 presents elementary results in sub-critical state investigated using phase-field modelling approach under isothermal conditions, such as stability of stagnant bubble, coalescence of two drops, phase-separation etc. The primary difference with regard to other models is the use of massfraction as the phase-field parameter.

Chapter 7 summarizes the work presenting the concluding remarks, highlighting the challenges and future perspectives derived from the current work.

CHAPTER 2. MATHEMATICAL MODELLING

The mathematical analysis of a fluid flow comprises of translating the physical laws of conservation, namely, mass, momentum and energy into the mathematical relations. These relations which are generally in the form of partial differential equations are termed as mathematical model of the fluid system. In case of a compressible fluid, such as near-critical fluids in the context of present work, an additional relation between pressure (𝑃), temperature (𝑇) and density (𝜌), known as the equation of state, is required to close the set of conservation equations. In general, it is inevitable to use a pressure-velocity coupling algorithm when solving the compressible model which renders the solution more intricate in addition to limiting the accuracy of density not reaching the machine precision. As in near-critical fluids, even a small variation in density can significantly affect the flow characteristics, it is thus highly desirable to calculate density with high precision.

In order to address these issues, a comprehensive description of a mathematical model to analyze single-phase highly compressible supercritical fluid is initially presented. The salient feature of the model lies in its ability to calculate density directly from the mass conservation (continuity equation) without the need of any pressure-velocity coupling algorithm. This is achieved by directly incorporating the dependence of pressure on density and temperature in the momentum equation. The development of such a model is motivated by the eventual aim of modelling the continuous transition from supercritical to subcritical state and vice-versa. For the analysis of two phase flow in sub-critical state, the phase-field modelling approach is introduced at the end of the chapter, discussing the theoretical background and contemporary work in this field.

One of the primary advantage of using phase field modelling approach is that it evades any geometrical reconstruction of the interfaces. Subsequently, a compressible phase-field model, which is in conjunction with prior developed compressible model, is described for isothermal systems.

Single-phase compressible model

Let us consider a fluid confined in a domain 𝛺 and bounded by a surface 𝛤, described by its The last part of ( It is to be noted here that the above set of equations are going to be solved numerically. This implies, the integral involving the divergence of velocity field can be simplified using several possible choices of numerical integration such as rectangle rule (constant value), trapezoidal rule (averaging at two-time steps) or Simpsons rule and so on. Here considering the simplification based on rectangle rule, for example, the integral ∫ including all the important physics essential to investigate any thermo-fluidic system. Some of the prominent features of the model can be summarized as,

• The momentum equation (2.1.14) is completely autonomous as it does not contain any unknown pressure. This circumvents the need of any pressure-velocity coupling algorithm as required in usual solution methodologies.

• The density is calculated directly from the continuity equation which ensures mass conservation resolvable to the machine precision.

• The dependence of pressure on density and temperature is explicitly incorporated in the momentum equation using continuity and energy equation as described in (2.1.7) and (2.1.8).

The thermodynamic properties in the current work are evaluated using relations obtained from Renormalization Group Theory as described in [START_REF] Zappoli | Heat Transfers and Related Effects in Supercritical Fluids[END_REF]. However, in general these may also be obtained from property data base such as NIST [START_REF]Thermophysical properties of pure fluids database[END_REF]. Further, in the current work the viscous dissipation is neglected while no volumetric source term is considered. Thus, the final set of governing equations as used in the current work can be summarized as,

𝜌 (

𝜕𝑽 𝜕𝑡 + 𝑽. 𝜵𝑽) = -𝜵 (𝑃 0 -( These equations are solved numerically using Finite Volume Method as will be described in CHAPTER 3.

Phase field model

The fluid in sub-critical regime is a two-phase system with singular behavior of various thermo-physical properties when approaching the critical point. This implies that unlike conventional two-phase fluid system, where the liquid phase is incompressible, in sub-critical regime both phases are highly compressible. Hence in order to analyze the behavior in subcritical state, a phase-field model is presented. Owing to preliminary analysis in sub-critical regime, attention has been restricted only to isothermal conditions 𝑖. 𝑒. without the energy equation.

Basics of Phase-field approach

The problem of multiphase flow has been intensively studied in literature owing to its critical importance in diverse systems, such as heat transfer systems (evaporators, condensers, electronic cooling systems etc.), power systems (Rankine cycle, internal-combustion engines etc.), process systems (chemical reactors, porous media etc.), environmental control (air conditioners etc.), biological systems (cardiovascular system etc.) and so on. The theoretical and numerical investigation of these systems is quite complex primarily attributed to the effect of surface tension and the existence of an interface. The mathematical approach to study two or multiphase system can be divided mainly into two categories, interface tracking and interface capturing. In interface tracking approach, such as VOF (volume of fluid method), each phase is identified using a continuum variable (such as color function in VOF). The effect of interface dynamics (surface tension) is directly incorporated into the governing equations by evaluating local gradient of the continuum variable. One of the main argument of this method, which is based on surface tension approach, is that material properties, such as density, exhibit a sharp discontinuity across the interface. However, as described by Gibbs [START_REF] Gibbs | On the equilibrium of heterogeneous substances[END_REF], the interface region is actually a rapid but smooth transition of physical quantities between the bulk fluid values. This idea was further developed by Rayleigh [START_REF] Rayleigh | XX. On the theory of surface forces.-II. Compressible fluids[END_REF] and Van der Waals [START_REF] Rowlinson | The thermodynamik theory of capillarity under the hypothesis of a continuous variation of density[END_REF] wherein thermodynamic principles were used to develop gradient theories for the interface. This concept was extended by Landau [START_REF] Landau | On the theory of phase transitions part i[END_REF][START_REF] Landau | On the theory of phase transitions part ii[END_REF] where the spatial description of material properties was described based on an order parameter (parameter which distinguishes one phase from the other), which laid an informal foundation of phase-field modelling. The interface thickness, however, is not arbitrary and is determined by the molecular force balance at the interface and its value is closely related to the range of molecular interactions [START_REF] Rowlinson | Molecular Theory of Capillarity[END_REF]. This is schematically illustrated in Fig. 2-1. Thermodynamically, the effect of these finite interactions can be represented in free energy by including the effect of the gradient of the order parameter. This was initially proposed by Cahn and Hilliard [START_REF] Cahn | Free Energy of a Nonuniform System. I. Interfacial Free Energy[END_REF][START_REF] Cahn | Free Energy of a Nonuniform System. III. Nucleation in a Two-Component Incompressible Fluid[END_REF], who derived an expression of free energy density in a non-homogenous system based on the concentration gradient and further used it to the study the problem of spinodal decomposition. Their model is now well known in literature as the Cahn-Hilliard model or Cahn-Hilliard equation. While density or composition may seem to be the intuitive choice of the order parameter (also termed as phase-field parameter), any other variable can also be used to define the two distinct phases [START_REF] Cahn | Free Energy of a Nonuniform System. I. Interfacial Free Energy[END_REF]. The bulk phase is thereby represented by a distinct constant value of the phase-field parameter while it varies smoothly in the interfacial region. Furthermore, it is known that any isolated system will tend to maximize its entropy or minimize its energy (for isothermal system) which implies that the free energy defined on the basis of phase-field parameters, must obey this thermodynamic constraint. This forms the basic principle in phase-modeling approach wherein the evolution of the phase parameter is obtained by satisfying the aforementioned thermodynamic constraint and since the interface can be uniquely located by smooth change in the phase-field variable, this interface is implicitly captured by studying the dynamics of the free energy function. As the interface is captured rather than tracked, phase-field method forms a part of the family of interface capturing method. It is to be mentioned that since these methods assume a finite interface thickness, which is ascribed to the molecular interaction or diffusion in the interface region, these have also been termed as diffuse interface models.

Contemporary work in phase-field modelling for fluid systems

The phase-field method was initially developed to study the behavior of the solidification of materials. Korteweg [START_REF] Korteweg | Sur la forme que prennent les ´equations du mouvements des fluides si l'on tient compte des forces capillaires causees par des variations de densite considerables mais continues et sur la theorie de la capillarite dans l'hypoth`ese d'une variation continue de la densit´e[END_REF] proposed that stresses and convection can be induced in a fluid system due to non-uniform density (or concentration or temperature) distribution. In order to include these effects, it was proposed to incorporate a specific stress tensor in the Navier-Stokes equations which is now well known in the literature as the Korteweg stress tensor. Since in phase-field modelling, non-homogeneity in the system is attributed to the gradient of the phase-field parameter, it is therefore possible to couple the phase-field and fluid equations using the Korteweg stress tensor based on the phase-field parameter. This idea made it possible to extend the phase-field modelling to two or multiphase flow problems following which phase-field modelling has been extensively used to study the various fluid flow phenomena varying from flow instabilities, bubble dynamics, thermo-capillary effect, phase change and so on. Here, some of the pioneering works covering various aspects are summarized.

The hydrodynamic coupling of phase-field models was initially studied by Gurtin et al. [START_REF] Gurtin | Two-phase binary fluids and immiscible fluids described by an order parameter[END_REF] where this methodology was used to analyze the coarsening effects in binary fluids. The problem of solidification with convection was addressed by Anderson et al. [START_REF] Anderson | A phase-field model of solidification with convection[END_REF]. Jacqmin [START_REF] Jacqmin | An energy approach to the continuum surface tension method[END_REF][START_REF] Jacqmin | Calculation of Two-Phase Navier-Stokes Flows Using Phase-Field Modeling[END_REF] presented how the phase-field approach can be extended to solve the full set of governing equations of fluid dynamics, 𝑖. 𝑒. Navier-Stokes equations, to study the various fluid phenomena such as break-up of bubbles, sloshing, wave breaking, contact lines etc. Unlike the Korteweg stress form, the coupling was introduced using the potential form of the free energy. Faraday instabilities were studied using phase-field modelling by Takagi and Matsumoto [START_REF] Takagi | Numerical simulation of two-dimensional Faraday waves with phase-field modelling[END_REF] and the results were verified with the linear theory. A similar analysis pertaining to the Rayleigh-Taylor instability was performed by Celani et al. [START_REF] Celani | Phase-field model for the Rayleigh-Taylor instability of immiscible fluids[END_REF] in case of two immiscible fluids in the limit of small Atwood numbers. Khatavkar et al.

[44] studied the behavior of the impact of a micron size drop on a smooth and flat solid surface using this approach. They investigated the effect of various material parameters such as impact velocity, droplet diameter etc. on the impact behavior and extended the analysis to non-90-degree contact angles. Borcia et al. [START_REF] Borcia | Static and dynamic contact angles -A phase field modelling[END_REF] used phase-field model to analyze the static and dynamic contact angles wherein the contact angle was controlled by boundary conditions for the density field (phasefield parameter) at the solid walls. Using the same approach, Borcia et al. [START_REF] Borcia | Drops on an arbitrarily wetting substrate: A phase field description[END_REF] further investigated the effect of variable wettability on the stability of a thin liquid film on a flat homogeneous solid and reported transition of a liquid film from a flat-film to a drop on a hydrophobic surface. In another work, Borcia et al. [START_REF] Borcia | Can Vibrations Control Drop Motion?[END_REF] investigated the effect of lateral and vertical harmonic oscillations on a liquid droplet resting on a solid plate and thus described a mechanism to control the motion of drops.

The phase-field model has also been largely used to study the phenomena of thermocapillary in fluid systems. By including the effect of advective transport in Cahn-Hilliard equation, Jasnow and Vinals [START_REF] Jasnow | Coarse-grained description of thermo-capillary flow[END_REF] studied the thermo-capillary induced flows, mainly the motion of droplets and phase separation via spinodal decomposition in a temperature gradient. A quasiincompressible phase-field model was derived by Antanovskii [START_REF] Antanovskii | A phase field model of capillarity[END_REF] for binary mixtures to study thermocapillary phenomena, where it was assumed that immiscible liquids can mutually penetrate into each other in the interface region in such a way that the sum of the mass diffusive flow rates of the two fluids vanishes. However, for incompressible flow, this led to a non-solenoidal velocity field in the interface region following which the model was termed as quasi-incompressible. Borcia

and Bestehorn [START_REF] Borcia | Phase-field model for Marangoni convection in liquid-gas systems with a deformable interface[END_REF] studied the Marangoni convection in a liquid-gas system with a deformable interface heated from below by describing the Korteweg stress term to be dependent on the temperature field. Similar to spinodal decomposition in solids, Lamorgese and Mauri [START_REF] Lamorgese | Diffuse-interface modeling of liquid-vapor phase separation in a van der Waals fluid[END_REF] studied the liquid-vapor phase separation in a van der Waals fluid subjected into the unstable range of its phase diagram. Owing to fundamentals of phase-field method arising from thermodynamic principles, Abels et al. [START_REF] Abels | Thermodynamically consistent, frame indifferent diffuse interface models for incompressible two-phase flows with different densities[END_REF] and Shen and Yang [START_REF] Shen | A Phase-Field Model and Its Numerical Approximation for Two-Phase Incompressible Flows with Different Densities and Viscosities[END_REF] derived the phase-field models for different density ratios from first principle approach. Ding et al. [START_REF] Ding | Diffuse interface model for incompressible two-phase flows with large density ratios[END_REF] described the model for two-phase incompressible fluid flows with large viscosity and density contrasts wherein the condition of solenoidal velocity field was derived from the conservation of mass of binary mixtures. The model differed from the previously derived quasi-incompressible models owing to the assumptions made regarding the relationship between the diffuse fluxes of the two species. Guo and Lin [START_REF] Guo | A thermodynamically consistent phase-field model for two-phase flows with thermocapillary effects[END_REF] derived a thermodynamically consistent model, i.e., based on the principles of maximization of the entropy both for compressible as well as incompressible fluid, though illustration was made only for incompressible fluids. While a review of developments pertaining to phase-field models can be found in [START_REF] Anderson | Diffuse-interface methods in fluid mechanics[END_REF][START_REF] Plapp | Phase-Field Models[END_REF] , other pioneering works regarding phase modelling have been addressed primarily from mathematical point of view in [START_REF] Fabrizio | A thermodynamic approach to non-isothermal phase-field evolution in continuum physics[END_REF][START_REF] Morro | A phase-field approach to non-isothermal transitions[END_REF][START_REF] Huang | A blow-up criterion for classical solutions to the compressible Navier-Stokes equations[END_REF][START_REF] Heida | On compressible Korteweg fluid-like materials[END_REF].

Relevance to the current work

Phase-field models hold a special significance when analyzing fluids near the critical point which is largely due to the assumption of finite thickness of the interface. This is in coherence with the physical behavior of interface thickness which shows a diverging behavior on approaching the critical point [START_REF] Wu | Critical Liquid-Vapor Interface in SF6. II. Thermal Excitations, Surface Tension, and Viscosity[END_REF][START_REF] Wu | Critical Liquid-Vapor Interface in SF6. I. Thickness of the Diffuse Transition Layer[END_REF][START_REF] Royall | Bridging length scales in colloidal liquids and interfaces from near-critical divergence to single particles[END_REF][START_REF] Beysens | Thickness of fluid interfaces near the critical point from optical reflectivity measurements[END_REF]. Thus, motivated by the physical nature of phase-field modelling in the context of near-critical fluids in sub-critical regime, a compressible model to analyze the flow behavior in sub-critical state in isothermal conditions is proposed in the current work. Appropriate changes are incorporated in the previous developed compressible model in §2.1, for a single-phase flow to couple the evolution of phase-field parameter and Navier-Stokes equations. The model can be further extended to non-isothermal systems as will be discussed in §7.2. In the subsequent section, the mathematical formulation of elementary phase-field equation is presented along with the description as to how the current model becomes different from the ones considered in literature, primarily attributed to the choice of phase-field parameter.

Mathematical description

The beginning point of the phase-field method corresponds to defining the free energy of the whole system which includes the effect of interfacial dynamics. This was initially presented by Cahn and Hilliard [START_REF] Cahn | Free Energy of a Nonuniform System. I. Interfacial Free Energy[END_REF] for a non-homogenous system based on the concentration though it can be written in terms of any intensive scalar property (called the order or phase-field parameter), denoted by 𝜙 in the current work. The order parameter thus helps to differentiate one phase from the other.

For example, considering density as a parameter, then each phase can be identified based on its density. The free energy (𝐹 𝑉 ) per unit volume in a generic form for a non-homogenous system 𝑤. 𝑟. 𝑡 𝜙 can be written in terms of 𝜙 as, (see Appendix E),

𝐹 𝑉 = 𝛽 ̅ 𝑓 0 (𝜙) + 𝛼 ̅ 2 |𝜵𝜙| 2 (2.2.1)
Here 𝑓 0 (𝜙) is the contribution to the free energy from the bulk while the second term incorporates the energy contribution due to change in 𝜙 across the interface. Further, 𝛼 ̅ and 𝛽 ̅ are constants related to physical parameters, surface tension (𝜎) and interface thickness (𝜉) (Appendices E.1 and E.2). Thus, it can be observed that there is an increase in the free energy of the system with the introduction of the interface which is implicitly represented by the gradient of 𝜙. However, governed by the thermodynamic principle in order to attain equilibrium, the system will adjust so that the overall free energy is always minimized for an isothermal system. In order to minimize the increase in free energy caused by the gradient term, an intuitive approach will be to consider a more diffused interface. However, this can only be achieved if there exists a higher concentration or contribution of 𝜙 in the bulk (near the interface region) thereby resulting in increase of the contribution of bulk free energy, eventually increasing the free energy of the system. Thus, these two effects are opposing to each other and equilibrium will exist as a compromise between these two effects. Mathematically, this can be described as a condition when the variation of free energy 𝑤. 𝑟. 𝑡 𝜙 will vanish. In the original work of Cahn and Hilliard [START_REF] Cahn | Free Energy of a Nonuniform System. I. Interfacial Free Energy[END_REF], as the order parameter under consideration was the composition, the variation of the free energy with composition, by definition, was termed as chemical potential. In coherence with same analogy, the chemical potential in the present case is defined as,

𝜂 = 𝛿𝐹 𝑉 𝛿𝜙 (2.2.2)
The potential can be understood as, how the free energy will change as the parameter 𝜙 changes in the system. Thus, the equilibrium conditions can be mathematically represented by 𝜂 = 0. Using (2.2.2), the equilibrium conditions can be thus represented as, (for a constant 𝛼 ̅) In general, the bulk free energy is represented by a double well potential such that each value at its minima uniquely represents each phase. Fig. 2-2 shows the profile of bulk free energy function, 𝑓 𝑜 (𝜙) = 𝜙 2 (1 -𝜙) 2 , considered in the present work, where minima corresponds to 𝜙 = 0 and 1.

𝜂 = 𝛽 ̅ 𝑓 𝑜 ′(𝜙) -𝛼 ̅ 𝛻 2 𝜙 (2.2.3)
The evolution equation of the phase-field parameter is governed by whether the phase-field parameter is conserved or not. In case of conserved phase-field parameter, which is the case in most fluidic systems, this is obtained from Cahn-Hilliard equation given by,

Coupling to Navier-Stokes equation

In surface tension methods, the effect of interface is directly forced into the momentum equation in terms of surface tension force. In order to introduce a similar effect in phase-field modelling, several approaches have been described in literature such as from the direct approach of Korteweg [START_REF] Korteweg | Sur la forme que prennent les ´equations du mouvements des fluides si l'on tient compte des forces capillaires causees par des variations de densite considerables mais continues et sur la theorie de la capillarite dans l'hypoth`ese d'une variation continue de la densit´e[END_REF], wherein a pre-assumed stress relation is directly introduced into the momentum equation.

With recent developments in the phase-field modelling, several authors have developed constitutive relations of stress in terms of phase-field parameter thereby inducing coupling between the phasefield and Navier-Stokes equations. Jacqmin [START_REF] Jacqmin | An energy approach to the continuum surface tension method[END_REF][START_REF] Jacqmin | Calculation of Two-Phase Navier-Stokes Flows Using Phase-Field Modeling[END_REF] derived the necessary changes in the momentum equation in terms of potential rather than stress form. The underlying idea behind this was based on the following facts,

• the movement of interface can cause it to become thick/thin thereby causing a change in free energy.

• this effect of the increase in energy is countered by a diffuse-interface force resulting in change of kinetic energy.

Based on this approach, it was found that the hydrodynamic coupling to include the effects of interfacial dynamics can be incorporated in terms of a force, which for compressible fluid can be defined as, (Appendix F),

𝑓 𝑆𝑇 = -𝜙𝛁𝜂 (2.2.6)
The above potential form can be converted to the stress form and the coupling is achieved by adding the following term to the momentum equation (see Appendix F),

-𝛼 ̅ 𝜵. (𝜵𝜙 ⊗ 𝜵𝜙) + 𝜵(𝛼 ̅|𝜵𝜙| 2 ) (2.2.7)
where the first term is known as the Korteweg stress tensor. A similar relation has been obtained by various authors based on different approaches. Gurtin [START_REF] Gurtin | Generalized Ginzburg-Landau and Cahn-Hilliard equations based on a microforce balance[END_REF] developed constitutive relations based on the physical reasoning that changes in interface region can be attributed to micro-forces acting in the interfacial region. A more general phase-field model based on thermodynamic principles was developed by Guo and Lin [START_REF] Guo | A thermodynamically consistent phase-field model for two-phase flows with thermocapillary effects[END_REF]. Among the stress and potential form, as described by Jacqmin , the stress form is more suited for dynamical situations as it ensures conservation of momentum [START_REF] Jacqmin | An energy approach to the continuum surface tension method[END_REF][START_REF] Jacqmin | Calculation of Two-Phase Navier-Stokes Flows Using Phase-Field Modeling[END_REF]. In the present work, the stress form will be used for the necessary coupling.

Choice of phase-field parameter

The phase-field parameter is used to identify the phase in two or multiphase system. A very intuitive choice in this context can be related to the density of the fluid under consideration. This may seem a good approach in case of incompressible fluids wherein density does not change and thus the stress term (Eq. (2.2.7)) will act only at the interface. However, some authors have also considered density as phase-field parameter in compressible fluids [START_REF] Lamorgese | Diffuse-interface modeling of liquid-vapor phase separation in a van der Waals fluid[END_REF][START_REF] Pecenko | A diffuse-interface approach to two-phase isothermal flow of a Van der Waals fluid near the critical point[END_REF][START_REF] Pecenko | Non-isothermal twophase flow with a diffuse-interface model[END_REF]. This is primarily facilitated by the fact that in compressible flow, the application of equation of state helps to solve for the pressure field while density can be calculated from continuity equation. This helps to evade the need to solve the phase-field equation. However, as in near critical fluids, the fluid is highly compressible, choosing density as a parameter will cause density gradients within the same phase and thus will not only affect the interface region but also the bulk phase which is not physical in nature. It is also to be mentioned here that unlike surface tension methods (such as VOF), where the forcing term acts only at interface (by virtue of evaluation of curvature at the interface), in phase-field method this force is indirectly incorporated by means of a stress tensor and acts based on the gradient of the phase-field parameter. This necessitates the need to have an appropriate phase-field parameter so that its gradient will persist only at the interface. Thus, in order to prevent the evoked issue, in the current work, the behavior of the compressible fluid (near the critical point) is being investigated by means of phase-field parameter which is considered to be the mass fraction of each phase. The motivation to choose this as the phase-field parameter is attributed to the following reasons,

• The gradient of phase-field parameter will be present only in the interface region as in the bulk, this term will be zero thereby precluding any effect that may arise due to density gradient in the bulk.

• In a compressible flow, it is necessary to ensure mass conservation. Thus, selecting mass fraction as the parameter ensures the mass conservation of the system in the bulk and interface region thereby providing a more physical meaning to the phase-field parameter.

Phase-field model with mass-fraction as order parameter

The phase-field model with parameter other than density has been proposed by various authors in literature such as volume fraction by Ding et al. [START_REF] Ding | Diffuse interface model for incompressible two-phase flows with large density ratios[END_REF], mass fraction by Lowengrub and

Truskinovsky [START_REF] Lowengrub | Quasi-incompressible Cahn-Hilliard fluids and topological transitions[END_REF], Abels et al. [START_REF] Abels | Thermodynamically consistent, frame indifferent diffuse interface models for incompressible two-phase flows with different densities[END_REF], Morro [START_REF] Morro | A phase-field approach to non-isothermal transitions[END_REF] and more recently by Guo and Lin [START_REF] Guo | A thermodynamically consistent phase-field model for two-phase flows with thermocapillary effects[END_REF]. Most of these have been used to analyze incompressible flows. Guo and Lin [START_REF] Guo | A thermodynamically consistent phase-field model for two-phase flows with thermocapillary effects[END_REF] and Morro [START_REF] Morro | A phase-field approach to non-isothermal transitions[END_REF] proposed a phase-field model based on mass fraction as one of the possible choice but these models have not been analyzed with any test cases. In the present work, we develop a phase-field model (i.e. Cahn-Hilliard equation) based on mass-fraction as discussed by Guo and Lin [START_REF] Guo | A thermodynamically consistent phase-field model for two-phase flows with thermocapillary effects[END_REF]. The model is coupled with the momentum equation using stress tensor as discussed in §2.2.5. However, unlike usual compressible models, where the density is explicitly calculated from the equation of state, the density here is calculated from mass conservation of each phase followed by the closure relation in terms of phase-field parameter.

Derivation of the phase-field model

The . Based on the above terminologies, the following relations can be described,

𝜌 ̃1 = 𝑀 1 𝑉 , 𝜌 ̃2 = 𝑀 2 𝑉
(2.2.8)

𝛾 1 + 𝛾 2 = 1 (2.2.9) 𝐶 1 + 𝐶 2 = 1 (2.2.10) 𝜌 ̃𝑖 = 𝑀 𝑖 𝑉 = 𝑀 𝑖 𝑀 𝑉 𝑀 = 𝛾 𝑖 𝜌 (2.2.11) 𝜌 = 𝜌 ̃1 + 𝜌 ̃2 (2.2.12)
Using Eq. (2.2.10), we have

𝑉 1 𝑉 + 𝑉 2 𝑉 = 1
which after simple calculations yields to,

𝛾 1 𝜌 1 + 𝛾 2 𝜌 2 = 1 𝜌 (2.2.13)
Now assuming that the two fluids move with different velocities 𝑽 𝒊 , the equation of mass balance for each material volume (in the entire domain) in terms of apparent density can be written as,

𝜕 𝜕𝑡 (𝜌 ̃1) + 𝜵. (𝜌 ̃1𝑽 𝟏 ) = 0 (2.2.14) 𝜕 𝜕𝑡 (𝜌 ̃2) + 𝜵. (𝜌 ̃2𝑽 𝟐 ) = 0 (2.2.15)
In the bulk region, the mass flux is ascribed only due to the advection while in the interfacial region between the two fluids, there also exists a contribution from the diffusion flux. In a physical sense this is driven by the gradient of the parameter (𝛾 1 ) across the interface and results in its smooth transition from one phase to another. Denoting the diffusive mass flux component of fluid 1 by 𝐽 𝑑1

(whose mathematical description will be provided shortly) , we can thus write, The diffusive flux here is attributed to the gradient of potential (𝜂) which exists near the interface region for a smooth transition of the parameter 𝜙 from one phase to another. Thus, we have 𝐽 𝑑1 = 𝑀 0 𝛻𝜂, which results in the same form as the Cahn-Hilliard equation (2.2.5) with advection and can be written as,

𝜌 ̃1𝑽 𝟏 = 𝜌 ̃1𝑽 -𝐽 𝑑1 (
𝜕𝜙 𝜕𝑡 + 𝑽. 𝛻𝜙 = 1 𝜌 𝛻. (𝑀 0 𝛻𝜂 ) (2.2.24)

Final set of governing equations for isothermal sub-critical flow

Based on mass fraction as the phase-field parameter, 𝜙, the final set of governing equations for the analysis of two-phase system in isothermal conditions is described now. It is important to mention here that in the current work, 𝜙 = 0 and 1 denote the two different phases, wherein 𝜙 = 1 denotes the liquid or heavier fluid (fluid 1) which in physical terms this implies the region is completely occupied by the liquid phase while 𝜙 = 0 denotes the gas or lighter fluid (fluid 2).

The free energy functional in terms of phase field parameter as described in (2.2.1) can thus be written as,

𝐹 𝑉 = 𝜌𝐹 = 𝜌𝛽𝑓 0 (𝜙) + 𝜌𝛼 2 |𝜵𝜙| 2 (2.2.25)
Here, 𝐹 is the free energy per unit mass, following which the potential, 𝜂 = 𝛿(𝜌𝐹) 𝛿𝜙 yields,

𝜂 = 𝜌𝛽𝑓 𝑜 ′ (𝜙) -𝛼 𝜵. (𝜌𝜵𝜙) (2.2.26)
The Cahn-Hilliard equation for the evolution of 𝜙 can thus be written as (𝑓 𝑜 ′ (𝜙) is written as 𝑓 𝑜 ′ for the sake of brevity),

𝜕𝜙 𝜕𝑡 + 𝑽. 𝜵𝜙 = 1 𝜌 𝜵. (𝑀 0 𝜵(𝜌𝛽𝑓 𝑜 ′ -𝛼 𝜵. (𝜌𝜵𝜙)) ) (2.2.27)
For the double-well bulk free-energy functional, 𝑓 𝑜 (𝜙) = 𝜙 2 (1 -𝜙) 2 , considered in the current work, the expression for potential can be written as,

𝜂 = 𝜌𝛽(4𝜙 3 -6𝜙 2 + 2𝜙) -𝛼 𝜵. (𝜌𝜵𝜙) (2.2.28)
following which the phase-field equation reduces to,

𝜕𝜙 𝜕𝑡 + 𝑽. 𝜵𝜙 = 1 𝜌 𝜵. (𝑀 0 𝜵(𝜌𝛽(4𝜙 3 -6𝜙 2 + 2𝜙) -𝛼 𝜵. (𝜌𝜵𝜙)) ) (2.2.29)

Coupling with the momentum equation

As described in §2.2.5, in order to couple the phase-field parameter with Navier-Stokes equation, the stress form of the interfacial force is added in the momentum equation (2.1.17) as described by,

𝐹 𝑆𝑇 = -𝛼 𝜵. (𝜌𝜵𝜙 ⊗ 𝜵𝜙) + 𝜵(𝛼𝜌|𝜵𝜙| 2 ) (2.2.30)

Closure relations

In order to ensure coherence between the evolution of phase-field variable and fluid properties, it is essential to write the density and other properties as a function of phase-field variable. This was described for the evolution of density in (2.2.13) which in terms of 𝜙 can be written as,

1 𝜌 = 𝜙 𝜌 1 + 1-𝜙 𝜌 2 (2.2.31)
Here 𝜌 1 and 𝜌 2 correspond to the densities of fluid 1 (𝜙 = 1) and fluid 2 (𝜙 = 0), respectively. It is to be noted here that, unlike an incompressible flow, the densities of both the fluids will change.

The density of each phase is evaluated based on its mass conservation in the entire domain as in case of single-phase flow (2.1.20) which yields,

𝜌 1 = 𝜌 1 0 𝑒𝑥𝑝 (-𝛿𝑡 𝜵 • 𝑽) (2.2.32) 𝜌 2 = 𝜌 2 0 𝑒𝑥𝑝 (-𝛿𝑡 𝜵 • 𝑽) (2.2.33)
following which they are advected as,

𝜕𝜌 1 𝜕𝑡 = 𝑑𝜌 1 𝑑𝑡 -𝑽 • 𝜵𝜌 1 (2.2.34) 𝜕𝜌 2 𝜕𝑡 = 𝑑𝜌 2 𝑑𝑡 -𝑽 • 𝜵𝜌 2 (2.2.35)
thereby leading to the density field given by (2.2.31). The other thermodynamic properties, namely isothermal compressibility and viscosity are closed in a similar manner as the density and can be written as,

1 𝜇 = 𝜙 𝜇 1 + 1-𝜙 𝜇 2
(2.2.36)

1 𝜒 𝑇 = 𝜙 𝜒 𝑇,1 + 1-𝜙 𝜒 𝑇,2 (2

Chapter summary

The mathematical description for the analysis of highly compressible near-critical fluids is developed in this chapter. With prime focus of this work on supercritical fluids, the model is initially developed for a single-phase fluid system. The model is described by the calculation of density directly from the mass conservation without the need of any pressure-velocity coupling algorithm. This is made possible by directly incorporating the dependence of pressure on density and temperature in the momentum equation. In order to investigate the behavior in sub-critical regime, a compressible model using phase-field modelling is described for isothermal systems.

In the next chapter, the numerical methodology to solve the developed mathematical model is described followed by the validation studies and the error analysis.

CHAPTER 3. NUMERICAL MODELLING AND VALIDATION TEST CASES

The Furthermore, the error analysis of the numerical scheme used in the current work (implicit secondorder central difference in space and forward first order Euler in time) is presented and it is shown that accurate results are possible with higher Courant numbers in non-linear system while low Courant numbers do not always imply accuracy.

Basics of Finite Volume Method

The mathematical representations of the fluid system in the form of partial differential equations can be discretized by various approaches such as finite difference method (FDM), finite volume method (FVM), finite element method (FEM) etc. Among these, FVM is profoundly used in problems related to the fluid flow as it ensures conservation of fluxes across the control volume.

The solution methodology is generally carried in the following steps:

➢ The starting point of FVM is the integral form of the governing equations. The physical domain is divided into number of small control volumes and conservation laws are applied to each control volume thereby ensuring that all fluxes are conserved. This helps to establish a precise relationship between the numerical methodology and conservation laws in a physical sense.

In the above procedure, a very important aspect pertains to defining an appropriate grid points for pressure (or scalar) and velocity fields. An intuitive approach is to define all the flow variables, 𝑖. 𝑒. pressure and velocity, on the same grid which can lead to uniform pressure field even though an oscillatory pressure may exist in the system [START_REF] Henk | An Introduction to Computational Fluid Dynamics: The Finite Volume Method[END_REF]. One possible solution is achieved by solving the pressure (and other scalar variables) at ordinary nodal points and velocity field on grid points which are staggered 𝑤. 𝑟. 𝑡 scalar nodes and centered on the cell faces. This is illustrated in Fig. 3-1 for 2D configuration, which shows the control volume for a scalar (such as pressure (𝑃), temperature (𝑇) etc.), 𝑥 -velocity and 𝑦 -velocity. The ordinary or scalar node under consideration is generally denoted by P while W, E, N and S represent scalar nodes on east, west, north and south of P, respectively. The velocity or flux nodes are marked on the west, east, south and north faces by w, e, s and n, respectively. Thus, the discretized equations for each variable is written for its respective cell as shown in Fig. 3-1.

It can thus be understood that the definition of a sub-volume or control volume in the domain is subjective in nature and depends on the flow variable under consideration, though these are usually identified based on scalar nodes.

➢ In the next step, the surface and volume integrals are transformed into discrete algebraic relations over each sub-volume (or cells). The variation of the flow variables is approximated by appropriate interpolation profiles relating the surface values to the cell values. The interpolation profiles or numerical schemes, as they are commonly referred to, are therefore prone to errors induced by approximations. In order to effectively capture the physical behavior of the flow phenomena, it is essential to make an appropriate choice of the numerical scheme. This is primarily governed by the type of flow being investigated.

Similarly, for unsteady problems, integration over a finite time step is performed by means of appropriate temporal discretization scheme. Depending on how the solution at the previous time step is passed to the current time step, these are classified as explicit, semiimplicit or implicit. The choice of temporal and spatial step is however not independent of each other and is governed by the propagation of error dynamics in the spatial-temporal plane. In order to ensure a minimal error, it is always desirable to have a neutrally stable scheme 𝑖. 𝑒. a scheme where the error growth rate (ratio of the error at the current and previous time step) is unity.

The discretized equations are then subjected to appropriate boundary conditions which represent physical conditions (such as heat flow, solid boundary wall 𝑒𝑡𝑐.). In physical terms, these conditions represent the constraint or the disturbances from the boundary propagating into the domain.

➢ The final step consists of solving the system of algebraic equations. The discretization process results in algebraic equations which may be linear or non-linear depending on the type of partial differential equations and numerical schemes. The size and set of the equations (which are solved in matrix form of type 𝐴𝑋 = 𝑏), is largely governed by the dimensionality and grid points used in the domain, though the matrices obtained from this discretization process are usually sparse 𝑖. 𝑒. most of their elements are zero. Even though, theoretically any valid method may be used to solve this system of equations, computational power often poses constraints on the solution process. In general, the resulting matrices can be solved by any of the two families of methods:

▪ Direct methods: These methods have a definite formula to evaluate the unknown variables 𝑖. 𝑒. they compute the solution in finite number of steps. For example, Gaussian Elimination, QR factorization etc.

▪ Iterative methods: In iterative methods, an initial guess to solution is made following which the solution is obtained once it converges below the preset criteria. Thus, the number of steps are not predefined in iterative methods. These methods are generally preferred in case of non-linear systems. Examples of these methods include Gauss-Seidel, Jacobi etc.

The above described methodology is incorporated in a programming language and used to solve the governing equations. In the current work, these are solved using the home-made code, Thetis, developed at Laboratory I2M which is described in the following section

Numerical code Thetis

The code Thetis is a FVM based code written in FORTRAN 90 for solving the fluid-flow problems in parallel computing. Some of its elementary features can be described as:

▪ Staggered grid: The scalar and velocity grids in Thetis follow a staggered arrangement. The scalar nodes are defined at the boundary of physical domain while the velocity mesh is shifted in the appropriate direction (for example in 𝑥 direction for velocity in 𝑥 -direction).

Thus, the computational domain corresponding to velocity field extends out of the physical domain. Fig. 3-2 shows the grid architecture as used in Thetis.

▪ Mesh: In addition to a simple constant mesh, Thetis has the capability to develop a nonuniform regular mesh of exponential nature. This type of mesh is particularly useful in regions such as near the walls to capture the effect of thin thermal boundary layer (as in the case of supercritical fluids in the current work) without the need to discretize the bulk domain for the same size. When initially defining the mesh, the grid points are created for the scalar nodes following which velocity nodes are interpolated from them.

▪ Vector form of the solution:

The momentum equation is solved in the vector form 𝑖. 𝑒. the components of the velocity field in all the directions (based on dimensionality of the problem) are evaluated simultaneously. In order to achieve this, the neighboring nodes of any node under consideration are located using connection matrices. These matrices help to build the coefficient matrix from the values of the flow variables, grid spacing 𝑒𝑡𝑐. based on the type of discretization scheme used in the analysis. The same methodology is applied for other variables for construction of their coefficient matrices. ▪ Solvers: Thetis provides parallel capabilities to solve the system of linear equations obtained after assembling the discretized governing equations. While both the direct and iterative solvers with preconditioner are available, in this work MUMPS (Multifrontal Massively Parallel Sparse direct Solver [START_REF] Amestoy | Multifrontal parallel distributed symmetric and unsymmetric solvers[END_REF]), a direct solver based on the Gaussian Elimination method for large sparse system of equations, is used for the simulations.

Thetis does not include the developed mathematical model by default. Thus, appropriate changes in the various subroutines are implemented to solve the desired model. The numerical model along with the numerical schemes used in the current work will be described now.

Numerical model

The numerical algorithm to solve the mathematical model is described herein. Initially, the single-phase model is presented followed by numerical methodology of phase-field model (Cahn-Hilliard equation and relevant changes in the momentum equation).

Single-phase flow

One of the salient feature of the model presented in CHAPTER 2 is that it does not involve any unknown pressure in the momentum equation which circumvents the need of any pressure velocity coupling and thus calculates pressure directly from the relation given by Eq. (2.1.10).

These equations are solved in a sequential manner. The second order space discretization is used for diffusion and transport term primarily motivated by its higher accuracy. Further, most of the analysis in the current work pertains to flow analysis in a closed system (cavity) and thus no high velocity fields or Péclet numbers (ratio of convective and diffusive flux), which govern the stability criteria for second order scheme are expected. The temporal behavior is discretized using an implicit first-order Euler scheme. A step by step description of the numerical algorithm is described below.

Step 1. Initial conditions: The initial conditions for the flow variables (𝑉 , 𝑃, 𝜌, 𝑇) at 𝑡 = 0 are known as given in the problem formulation.

Step 2. Discretization of Energy equation: The energy equation solves the temperature field (𝑇) implicitly as:

𝜌 (𝑛) 𝐶 𝑉 (𝑛) ( 𝑇 (𝑛+1) -𝑇 (𝑛) 𝛿𝑡 + 𝑽 (𝑛) . 𝜵𝑇 (𝑛+1) ) = - 𝛽 𝑃 (𝑛) 𝑇 (𝑛) 𝜒 𝑇 (𝑛) 𝜵 • 𝑽 (𝑛) + 𝜵 • (𝑘 (𝑛) 𝜵𝑇 (𝑛+1) ) (3.3.1)
Step 3. Discretization of Momentum Equation :The momentum equation given by (2.1.17) is solved implicitly for the velocity field in vector form as, 𝜌 (𝑛) ( 𝑉 (𝑛+1) -𝑉 (𝑛) 𝛿𝑡 + 𝑽 (𝑛) . 𝜵𝑽 (𝑛+1) ) = -𝜵 (𝑃 0(n) -𝛿𝑡 (

1 𝜒 𝑇 (𝑛) + 𝛽 𝑃 2(𝑛) 𝑇 (𝑛+1) 𝜌 (𝑛) 𝐶 𝑉 (𝑛) 𝜒 𝑇 2(𝑛) ) 𝜵 • 𝑽 𝑛+1 -𝛿𝑡 ( 𝛽 𝑃 (𝑛) 𝜌 (𝑛) 𝐶 𝑉 (𝑛) 𝜒 𝑇 (𝑛) ) 𝜵. ((-𝑘 (𝑛) )𝜵𝑇 (𝑛+1) )) + 𝜵 • [𝜇 (𝑛) (𝜵𝑽 𝑛+1 + 𝜵 𝒕 𝑽 𝑛+1 - 2 3 𝜵 • 𝑽 𝑛+1 𝐼)] (3.3.2)
Here, the updated temperature field, 𝑇 (𝑛+1) , obtained by solving (3.3.1) is used.

Step 4. Calculation of Pressure: The equilibrium Lagrangian pressure 𝑃 𝐿 for time step (𝑛 + 1) can thus be calculated from this new velocity field 𝑉 𝑛+1 and temperature field 𝑇 𝑛+1 as,

𝑃 𝐿 0(𝑛+1) = 𝑃 0(n) -( 1 𝜒 𝑇 (𝑛) + 𝛽 𝑃 2(𝑛) 𝑇 (𝑛+1) 𝜌 (𝑛) 𝐶 𝑉 (𝑛) 𝜒 𝑇 2(𝑛) ) 𝜵 • 𝑽 (𝑛+1) 𝛿𝑡 - 𝛿𝑡 ( 𝛽 𝑃 (𝑛) 𝜌 (𝑛) 𝐶 𝑉 (𝑛) 𝜒 𝑇 (𝑛) ) 𝜵. ((-𝑘 (𝑛) )𝜵𝑇 (𝑛+1) ) (3.3.3)
The above pressure relation was derived from its total derivative (or Lagrangian form) and hence the subscript is used to differentiate it from 𝑃 0(𝑛+1) which is obtained as explained below.

Step 5. Calculation of density: The updated velocity field 𝑽 at (𝑛 + 1) is further used to evaluate the Lagrangian form of the density (𝜌 𝐿 ) at time step (𝑛 + 1) using the following relation (as in

(2.1.20)), 𝜌 𝐿 (𝑛+1) = 𝜌 (𝑛) 𝑒𝑥𝑝 (-𝛿𝑡 𝜵. 𝑽 (𝑛+1) ) (3.3.4)
Here similar analogy as pressure is used to define 𝜌 𝐿 .

Step 5. Advection of flow variables: In order to ensure coherence between behavior of flow variables in terms of Eulerian and Lagrangian form, the density and pressure fields are then advected from their respective total derivatives as in (2.1.21),

( 𝜌 (𝑛+1) -𝜌 (𝑛) 𝛿𝑡 ) = ( 𝜌 𝐿 (𝑛+1) -𝜌 𝐿 (𝑛) 𝛿𝑡 ) -𝑽 (𝑛+1) . 𝜵𝜌 𝐿 (𝑛+1) (3.3.5) ( 𝑃 (𝑛+1) -𝑃 0(𝑛) 𝛿𝑡 ) = ( 𝑃 𝐿 (𝑛+1) -𝑃 𝐿 (𝑛) 𝛿𝑡 ) -𝑽 (𝑛+1) . 𝜵𝑃 𝐿 (𝑛+1) (3.3.6)
Step 6. Updating the properties: The thermodynamic properties as a function of temperature and other variables are updated using the flow variables obtained at 𝑛 + 1.

The above process is repeated in time by setting the values at time (𝑛 + 1) to the values at time (𝑛). In solving the above equations, the matrix system is solved using direct solver MUMPS [START_REF] Amestoy | Multifrontal parallel distributed symmetric and unsymmetric solvers[END_REF].

Phase field Model

The phase field model to analyze two-phase flow in subcritical state was described in §2.2.

It comprises primarily of Cahn-Hilliard form of equation for the evolution of phase-field parameter 𝜙. The numerical discretization of the set of equations described in §2.2.8 is presented herein.

In step-by-step algorithm described in §3. The phase-field equation is solved implicitly with first-order Euler in time and while QUICK scheme is used for spatial discretization. The choice for higher order spatial discretization scheme is motivated by a fourth order derivative in the phase-field equation (2.2.29) and thus include the influence from more neighboring mesh points to reduce the discretization errors by bringing them in a wider influence. Further, due to the double well form of bulk-free energy function, the phase-field equation includes a non-linear term. In order to solve this in implicit manner, the non-linear terms are split as shown below. The phase-field equation, is thus discretized as, ( 𝜙 (𝑛+1) -𝜙 (𝑛) 𝛿𝑡 + 𝑽 (𝑛+1) . 𝜵𝜙 (𝑛+1) ) = The density of each fluid (or phase 1 and 0) is calculated in similar manner as for a single-phase flow,

𝜌 𝐿,1 (𝑛+1) = 𝜌 𝐿,1 (𝑛) 𝑒𝑥𝑝 (-𝛿𝑡 𝜵. 𝑽 (𝑛+1) ) (3.3.9)
𝜌 𝐿,2 (𝑛+1) = 𝜌 𝐿,2 (𝑛) 𝑒𝑥𝑝 (-𝛿𝑡 𝜵. 𝑽 (𝑛+1) ) (3.3.10)

Here, the subscript 𝐿 denotes the Lagrangian form of the density as was described in §2.1 for a single-phase fluid. These are then advected from their total derivatives as,

( 𝜌 1 (𝑛+1) -𝜌 1 (𝑛) 𝛿𝑡 ) = ( 𝜌 𝐿,1 (𝑛+1) -𝜌 𝐿,1 (𝑛) 𝛿𝑡 ) -𝑽 (𝑛+1) . 𝛁𝜌 𝐿,1 (𝑛+1) (3.3.11) ( 𝜌 2 (𝑛+1) -𝜌 2 (𝑛) 𝛿𝑡 ) = ( 𝜌 𝐿,2 (𝑛+1) -𝜌 𝐿,2 (𝑛) 𝛿𝑡 ) -𝑽 (𝑛+1) . 𝛁𝜌 𝐿,2 (𝑛+1) (3.3.12)
The density and other thermo-physical properties for the next iteration are calculated using the closure relation given by (2.2.31), (2.2.36) and (2.2.37) as,

1 𝜌 (𝑛+1) = 𝜙 (𝑛+1) 𝜌 1 (𝑛+1) + 1-𝜙 (𝑛+1) 𝜌 2 (𝑛+1) (3.3.13) 1 𝜇 (𝑛+1) = 𝜙 (𝑛+1) 𝜇 1 + 1-𝜙 (𝑛+1) 𝜇 2 (3.3.14) 1 𝜒 𝑇 (𝑛+1) = 𝜙 (𝑛+1) 𝜒 𝑇,1 + 1-𝜙 (𝑛+1)
𝜒 𝑇,2

(3.3.15)

Validation of the Numerical Model

Single-phase model

The mathematical model developed in §2.1 and solved using numerical algorithms described in §3.3.1 is validated by analyzing the propagation of acoustic wave in water (isothermal conditions) and thus only momentum and continuity equations are solved. The results obtained using this were further extended to explain how propagation of acoustic wave in water helps in walking mechanism of basilisk lizard in Sharma et al. [START_REF] Sharma | Numerical modeling of the impact pressure in a compressible liquid medium: application to the slap phase of the locomotion of a basilisk lizard[END_REF]. The validation of the complete model has already been described in [START_REF] Amiroudine | A Lagrangian-Eulerian compressible model for the trans-critical path of near-critical fluids[END_REF] wherein density change caused by piston effect is compared with experimental observations of Miura et al. [START_REF] Miura | High-speed observation of the piston effect near the gas-liquid critical point[END_REF].

Problem description: propagation of acoustic wave in water

Consider a 1D container filled completely with water and open only from the top to permit inflow into the system as shown in Fig. 3-3. The flow of water into the system from the top results in compression of the water present in the system to accommodate the mass entering in the fixed volume of the container. This leads to the generation and propagation of pressure waves in the system. For a pressure wave front moving at the speed of sound 𝑐 𝑠 , and fluid velocity 𝑉 0 , the pressure change across the control volume can be expressed as ∆𝑃 = 𝜌𝑐 𝑠 𝛥𝑉 [START_REF] White | Fluid Mechanics[END_REF] (Appendix A for details). This expression can also be re-written as:

𝑃 = 𝜌𝑐 𝑠 𝑉 0 (3.4.1)
while describing changes 𝑤. 𝑟. 𝑡 some stationary end of the control volume. The time taken by the pressure wave to reach back to the inlet is given by (2 𝐻 0 /𝑐 𝑠 ) where 𝐻 0 is the depth of the chamber.

For the depth considered here (𝐻 0 = 0.5𝑚) the analytical time for the pressure wave to return to the inlet is thus 0.67 × 10 -3 𝑠. Hence, to numerically capture the propagation of these acoustic waves, a timestep smaller than 10 -4 𝑠 is selected. Table 3-1 describes the physical parameters required in this study.

Analysis with different Courant Numbers

Since the numerical schemes are prone to errors, it is essential to analyze their effect on the accuracy of the solution. In an unsteady problem, the choice of time-step and grid size is dependent on each other for the stability and accuracy of the solution. Thus, it is essential to analyze the behavior of the numerical scheme before validation of the numerical model. In order to achieve this, the propagation of acoustic wave in water is analyzed for different grid sizes and Courant numbers. In terms of time steps, these cases correspond to 10 𝜇s, 1 𝜇s and 0.1 𝜇s, respectively. Thus, for the highest Courant number used here, the time scales resolved is the finest. One notices oscillations in several cases, particularly for 𝑁 𝑐 < 1 while non-oscillatory stable solutions are observed for 𝑁 𝑐 = 10. The observed behavior seems to contradict the conventional argument (based on the analysis of [START_REF] Charney | Numerical Integration of the Barotropic Vorticity Equation[END_REF]) wherein stability of the solution is expected to decline with increasing values of Courant number [START_REF] Livne | An implicit method for two-dimensional hydrodynamics[END_REF][START_REF] Briley | Solution of the multidimensional compressible Navier-Stokes equations by a generalized implicit method[END_REF]. This unusual observation needs to be investigated and is discussed now. However, before proceeding further, a brief description to error analysis is presented herein.

Basics of error analysis: A brief introduction

The numerical solution to any PDE representing a physical system consists of discretizing the PDE over the physical domain using appropriate spatial differencing scheme such as central difference, upwind etc. and temporal scheme such as first order Euler scheme. However, in the process of discretization, different types of errors creep into the solution which may lead to erroneous conclusion from the numerical results. On the other hand, certain schemes may become unstable under certain combinations of spatial grid spacing and time step. This may eventually lead to growth of error in time, oscillation in solutions etc. Hence, it becomes necessary to study the properties and behavior of the numerical scheme and errors induced due to the same.

When the physical system is solved numerically, the solution can be thought of as numerical waves travelling in the spatial and temporal domain. This can be obtained by transforming the solution in Cartesian domain to spectral domain of wavenumber and frequency. In order to describe the behavior of the numerical scheme, the following terms can be defined as:

• Phase Velocity Phase velocity (𝑐 𝑁 ) refers to velocity of individual wave.

• Group Velocity

Group velocity ( 𝑉 𝑔𝑁 ) refers to the velocity of the entire group of waves as a whole, here denoting numerical waves. The medium through which these waves travel may be characterized by the relation between 𝑐 𝑁 and 𝑉 𝑔𝑁 .

▪ 𝑐 𝑁 = 𝑉 𝑔𝑁 : non-dispersive medium

▪

Error propagation equation

In order to obtain accurate numerical solution in transient analysis, implicit schemes are known to be unconditionally stable for linear systems while non-linearity in advection systems often limits the stability. Nevertheless, the implicit schemes for non-linear systems allow larger Courant numbers than explicit schemes. In the methodology of von Neumann [START_REF] Charney | Numerical Integration of the Barotropic Vorticity Equation[END_REF], it is assumed that the propagation of signal (𝑖. 𝑒. how a variable evolves in space and time) and the corresponding error follows the same dynamics. This assumption may intuitively seem to be satisfying for linear as well as quasi-linear form of non-linear equations. However, in recent works of Sengupta et al.

[77], it has been shown that this does not hold true owing to the dispersion or phase error or when the numerical method is not strictly neutrally stable for the model 1D convection equation. In order to illustrate this, consider the 1D wave equation which is the classical linear system used for error analysis,

𝜕𝑢 𝜕𝑡 + 𝑐 𝜕𝑢 𝜕𝑥 = 0, 𝑐 > 0 (3.4.2)
𝑐, is the wave speed. The error is defined as: 𝑒 = 𝑢(𝑥, 𝑡) -𝑢 𝑁 (𝑥, 𝑡), with 𝑢(𝑥, 𝑡) and 𝑢 𝑁 (𝑥, 𝑡) representing the exact and numerical solutions, respectively. The error equation can thus be written as (see Appendix C for details),

𝜕𝑒 𝜕𝑡 + 𝑐 𝜕𝑒 𝜕𝑥 = -𝑐 [1 - 𝑐 𝑁 𝑐 ] 𝜕𝑢 𝑁 𝜕𝑥 -∫ 𝑑𝑐 𝑁 𝑑𝑘 [∫ 𝑖𝑘 ′ 𝑈 0 [|𝐺|] 𝑡 ∆𝑡 𝑒 𝑖𝑘 ′ (𝑥-𝑐 𝑁 𝑡) 𝑑𝑘 ′ ] 𝑑𝑘 - ∫ 𝐿𝑛|𝐺| ∆𝑡 𝑈 0 [|𝐺|] 𝑡 ∆𝑡 𝑒 𝑖𝑘(𝑥-𝑐 𝑁 𝑡) 𝑑𝑘 (3.4.3)
It can be seen that unlike the conventional error dynamics, which comprises only of the error growth rate |𝐺| (ratio of the error at the current and previous time step), the error propagation comprises of several terms. Thus, the characteristics of the numerical scheme for the lowest error can be primarily evaluated in terms of the following parameters [START_REF] Sengupta | Error dynamics: Beyond von Neumann analysis[END_REF][START_REF] Sengupta | High Accuracy Compact Schemes and Gibbs' Phenomenon[END_REF][START_REF] Sengupta | High Accuracy Computing Methods: Fluid Flows and Wave Phenomena[END_REF]: ) and numerical group velocity identical to the phase speed. An important aspect worth mentioning here is the difference between the physical wave speed (𝑐) and the numerical wave speed (𝑐 𝑁 ). The physical wave speed is the actual wave speed which is attempted to be captured in computational domain. Thus, in the numerical domain, the waves travel with speed 𝑐 𝑁 and in ideal scenario 𝑐 = 𝑐 𝑁 , an assumption in classical error analysis [START_REF] Charney | Numerical Integration of the Barotropic Vorticity Equation[END_REF]. However, based on the numerical scheme and the resulting dispersion relation, these

•
are not equal and this contributes to the error as described by Eq. (3.4.3).

Error characteristics of the current numerical scheme

In order to investigate the behavior of the numerical scheme used in the current work (second order space discretization and first order Euler in time), it is important to analyze its characteristics which is done using the 1D wave equation (3.4.2). This is further used to explain the unusual behavior pertaining to stable solutions at high Courant numbers and inaccurate solutions at lower ones (Fig. 34).

The numerical amplification rate for Eq. It therefore provides a possibility to analyze the propagation of error based on the chosen time step and grid size for the numerical analysis. It can be inferred that the selection of an appropriate operating point (𝑖. 𝑒. choice of spatial and temporal resolution) is not only dependent on the Courant Number but is also governed by the appropriate combination of the three parameters represented in Fig. 345. This means an erroneous and oscillating solution is possible even for Courant number less than or equal to unity owing to the contribution of various terms in Eq. (3.4.3) It is to be mentioned here that the error analysis (characteristics) presented herein for this wave equation does not include the advective and viscous terms. Though it may seem reasonable to question the applicability of the results obtained to explain the observations for the complete set of Navier-Stokes equations, it can be argued that the simple wave equation (in the form given by Eq. (3.4.2)), could represent the whole problem in the following manner: in the case of inclusion of viscous terms in the error analysis of wave equation, the numerical oscillations would be dampened out while an opposite effect can be expected from the nonlinear advection term. The net effect of both these terms will either lead to an increase or decrease of the stability zone but the qualitative representation of the numerical characteristics can be asserted to remain the same. The changes that would be encountered, can thus be stated to affect the characteristics more in a quantitative way.

Moreover, here the analysis pertains to investigating how effective is the current numerical scheme (as described above) to capture the dynamics governed by the Navier-Stokes equations. Thus, if the numerical scheme is unable to give precise and accurate results for a simple wave equation, it is highly unlikely to get precise and accurate results for the complete set of Navier-Stokes equations.

Hence, for the sake of clarity, the characteristics obtained by correct error dynamics for simple wave equation are well suited to explain the observed results in Fig. 34.

Discussion of observations for different Courant Numbers

It is observed in Fig. 3-4 that highly accurate and stable results are obtained for high Courant number while their low values showed some oscillatory behavior. In order to explain this on the basis of numerical characteristics of the current scheme, Fig. 3456 ) respectively, in the (𝑁 𝑐 , 𝑘Δ𝑥)-plane. It is to be mentioned that the largest to smallest grid size corresponds to a time step of 10 𝜇𝑠, 1 𝜇𝑠 and 0.1 𝜇𝑠 for each Courant number (similar plot is also depicted in Fig. 34567) for 𝑁 𝑐 =10). In order to mark the considered operating points in (𝑁 𝑐 , 𝑘∆𝑥)-plane in the contour plots, the fundamental spatial frequency of the signal (

𝑃 𝑃 𝑎𝑡𝑚

) for each combination of 𝑁 𝑐 and ∆𝑥 is evaluated for three different time instances (0.3 𝑚𝑠, 0.6 𝑚𝑠 and 0.9 𝑚𝑠). Using Fourier transform of the output signal (

𝑃 𝑃 𝑎𝑡𝑚

), the relevant value of peak wavenumber is found to be 𝑘 2𝜋 ~ 1 𝑚 -1 . However, the pressure pulse propagating will behave like a Heaviside function, 𝑖. 𝑒., the wavenumber spectrum will be very wide and corresponding wavenumber will require fine spacing to capture discontinuous behavior (as shown in Fig. 34), with the variation now considered in space. In that respect, the case of 𝑁 𝑐 = 0.11 for ∆𝑥 = 125𝑚𝑚 will be very difficult. Comparatively, for the nominal fundamental case of 𝑘 2𝜋 ~ 1 𝑚 -1 , one seeks the properties for 𝑘∆𝑥 = 0.7854. The solutions are found to be quite satisfactory, except for larger grid size (∆𝑥 = 125𝑚𝑚) where the solution deviates significantly from the exact solution (see Fig. 34, 𝑁 𝑐 = 0.11). An explanation for such a deviation is presented now.

Consider firstly the case with 𝑁 𝑐 = 0.11, where oscillations at the initial time are observed for ∆𝑥~125𝑚𝑚, for which the contour plots are zoomed in Fig. 3456). The computed oscillatory solution is due to what is known as the Gibbs' phenomenon [START_REF] Sengupta | High Accuracy Compact Schemes and Gibbs' Phenomenon[END_REF] and is related to peculiar behavior of the numerical group velocity for high value of 𝑘∆𝑥. This is attributed to the upstream propagating q-waves to be noted for high value of 𝑘∆𝑥. (In general, in order to capture the physical wave numerically, it is essential that both physical and numerical waves travel in the same direction.

However, due to the characteristics of the numerical scheme, some of the numerical waves can start moving in the opposition direction thereby causing an unphysical solution and such waves been termed as 𝑞 -waves [START_REF] Sengupta | High Accuracy Compact Schemes and Gibbs' Phenomenon[END_REF][START_REF] Sengupta | High Accuracy Computing Methods: Fluid Flows and Wave Phenomena[END_REF]). For other grid sizes (Fig. 3456)) such high values of 𝑘∆𝑥 do not arise.

For the numerical stability, implicit schemes are not constrained with respect to the Courant number and none of the computed cases display error growth. However, a comparison of the results, as shown in Fig. 34(a) and Fig. 34), illustrates that an increasing Courant number leads to more accurate solutions and capture a sharp step-like profile (Heaviside function) more efficiently for the same time step. It is observed from Fig. 34567(a) that the error growth rate for all grid sizes is very close to unity for 𝑁 𝑐 = 10, which explains the accuracy of the solution obtained for higher Courant number. Similarly, a nearly vanishing value of phase error in Fig. 34567(c), reveals that numerical waves travel with nearly the same speed as the physical waves and observance of accurate solution.

The fact that numerical waves follow the physical waves nearly at the same speed is further evident from the ratio of 𝑉 𝑔𝑁 𝑐 being close to unity. This explains a stable and accurate solution observed in Fig. 34for higher Courant numbers. ). The marked points correspond to the fundamental wavenumber

𝑘 2𝜋 ~ 1 𝑚 -1 .
It is to be mentioned that even though it may seem obvious to obtain an oscillation-free solution for higher Courant numbers (for the same time step) owing to a better spatial resolution, a deviation from a sharp step-like profile is still observed at time step Δ𝑡 = 10 𝜇𝑠 (largest grid spacing ∆𝑥 = 1.49 𝑚𝑚 for 𝑁 𝑐 = 10). The observed lack of accuracy is due to a lower grid resolution, while the step-like function requires resolving larger range of wavenumbers.

As illustrated in §3. These extreme cases have been purposely selected so as to elucidate the results effectively in a quantitative manner, though a similar explanation will hold for all the cases showing such discrepancies. The former case corresponds to a deviation from the step like profile while the latter represents the desired sharp peak as shown in Fig. 34). Fig. 345678illustrates the point corresponding to the second dominant wavenumber in (𝑁 𝑐 , 𝑘Δ𝑥)-plane for 𝑁 𝑐 = 0.11. It can be seen that while the value of error growth rate (|𝐺|) is nearly equal to 1, the negative value of 𝑉 𝑔𝑁 /𝑐 and high value of phase lag implies that the constituting numerical waves are not in coherence with the physical wave. Specifically, when the solution is discontinuous, a large range of 𝑘∆𝑥 is excited.

Along with this comes the problem of negative group velocity of the components with 𝑘∆𝑥 exceeding 𝜋/2 (as shown in Fig. 3-5 and a zoomed view shown in Fig. 3-8 for the case of ∆𝑥 = 125𝑚𝑚 and the corresponding 𝑘∆𝑥 ≈ 𝜋). For the highest wavenumbers, the components of solution will move upstream from the point of discontinuity, as shown in Fig. 34. This is the so-called Gibbs' phenomenon, which is described in greater detail in [START_REF] Sengupta | High Accuracy Compact Schemes and Gibbs' Phenomenon[END_REF]. Thus, the net signal (

𝑃 𝑃 𝑎𝑡𝑚

) is not well reconstructed leading to a deviation from the desired value as observed in Fig. 34. On the contrary, it can be observed from Fig. 3456789(for 𝑁 𝑐 = 10 and ∆𝑥 ≈ 0.0149 mm) that for various higher dominant wavenumbers, the values of characteristics are close to the desired ideal values, especially the phase lag error, 1 -𝑐 𝑁 𝑐 unlike for 𝑁 𝑐 = 0.11 where it is significantly high (Fig. 345678). Even though higher wavenumbers get excited 57 like the previous case, finer spatial resolution impedes group velocity from moving upstream. This explains an accurate solution observed in Fig. 34) and thus a sharp step like profile which satisfies the analytical solution.

It can also be observed that all the three parameters for 𝑁 𝑐 = 0.11 and ∆𝑡 = 0.1 𝜇𝑠 (∆𝑥 ≈ 1.47 𝑚𝑚) are closer to the desired ideal values as compared to 𝑁 𝑐 = 10 for the same time step and ∆𝑥 ≈ 0.0149 𝑚𝑚. However, from the results shown in Fig. 34, it can be seen that solutions for the latter are more stable 𝑖. 𝑒. no oscillations for time instances very close to 𝑡 = 0, which may give an impression of an anomalous behavior. As mentioned before, the marked operating points correspond to time instances beginning from 0.3 𝑚𝑠 which thereby preclude the operating points denoting initial instances of the solution.

58

In order to explain this observation, for the sake of completeness, the signal at time instances 𝑡 = 0.5, 1, 1.5 and 2𝜇𝑠 is investigated for the fundamental wavenumber for 𝑁 𝑐 = 0.11 and 10 with ∆𝑡 = 0.1 𝜇𝑠. It is found that in both the cases, all the possible wavenumbers were nearly equally plausible. It is to be mentioned here that for the same time step, a higher Courant number entails a wider spectrum of wavenumbers. This implies that for 𝑁 𝑐 = 10, the wavenumbers available to reconstruct/capture the signal (i.e. propagation of numerical wave) are considerably higher.

Consequently, the output thus observed is more coherent with the exact signal (propagation of the physical wave) leading to no-oscillations in the solution.

On the contrary, the possible wavenumbers for 𝑁 𝑐 = 0.11 (for same time step as 𝑁 𝑐 = 10) are considerably lower resulting in higher values of wavenumber (𝑘) and thus q-waves and relatively poor reconstruction of the signal. Thus, oscillations are observed in the solution as can be observed in Fig. 34. The above explanation therefore suggests that the lower values of Courant Number do not necessary entail stable solution while higher values can lead to stable solutions. However, obtaining stable solution for higher Courant numbers as compared to lower values is unusual and has been highlighted in the above discussion. This is attributed to the correct and complete error dynamics which includes effects of numerical error growth such as dispersion and phase lag errors described above.

It is to be noted here that Courant Number has been defined based on the actual wave speed 𝑖. 𝑒. the physical speed of the wave whose propagation is being captured numerically. Intuitively, it seems reasonable to ascertain that for an accurate numerical analysis of wave propagation, it is desirable that the smallest grid size in computational domain be greater than the distance traversed by the wave in a given time step i.e. 𝑐∆𝑡 < ∆𝑥. However, in the computational domain, the physical waves are actually represented by the propagation of numerical waves which move with the phase speed 𝑐 𝑁 and depend on wavenumber for dispersive schemes. Hence, the intuitive logic of grid size being larger than the distance moved by the wave needs to be written as 𝑐 𝑁 ∆𝑡 < ∆𝑥 which will be dependent on the wavenumber. Moreover, as the solution at a point is the summation of solution due to all wavenumbers, (when 𝜔 = 𝑐(𝑘)), therefore defining constraint of this sort can cause confusion. Instead, an analysis of the kind presented here advocates appropriate grid size and time step. Further, the results presented herein bring out an important fact that the criterion of numerical stability (Courant numbers lesser than one) is dependent on the numerical scheme, as well as, the type of PDE being investigated.

Validation results for the propagation of acoustic wave in water

The problem to investigate propagation of acoustic wave in water is explained in § 3.4.1.1 while results for various Courant numbers are shown in Fig. 34. The knee-like jump after a certain time (two acoustic times) represents the addition effect of the pressure waves which reflect from the bottom wall and interfere with the pressure wave at the entrance. This is further illustrated in Thus, the distance traversed by pressure wave (as highlighted by point 𝑀 1 ) shows the speed of pressure-wave propagation is the same as the speed of sound in water. It is observed that with the increase in resolution in space and appropriate choice of the time step, the error between the analytical and numerical time (2 𝐻 0 /𝑐) becomes small (around 3-4% for 𝛿𝑡 = 0.1 𝜇𝑠 and 𝑁 𝑐 =1 and 0.1). Furthermore, for 𝛿𝑡 = 0.1 μs, the pressure error at time 0.2𝑚𝑠 was about 1.5%. This shows that the above mathematical and numerical models can be used to study the propagation of pressure waves even in media with a very low compressibility (water in this case).

Phase-field (Cahn-Hilliard) equation

In this section, validation results for the Cahn-Hilliard equation (2.2.29) are described based on mass-fraction as the order parameter. Cahn and Hilliard [START_REF] Cahn | Free Energy of a Nonuniform System. I. Interfacial Free Energy[END_REF][START_REF] Cahn | Free Energy of a Nonuniform System. III. Nucleation in a Two-Component Incompressible Fluid[END_REF] proposed this equation to describe the spontaneous phase separation in binary fluid mixtures (binary alloy). A brief description of the governing mechanism is initially presented followed by the description of the results. Before proceeding further, it is important to mention two important parameters which play a crucial role in numerical aspects of phase-field modelling. These are,

• Cahn-number (𝐶𝑛): Cahn number is the ratio of the interface thickness (𝜉) to the characteristic length scale (𝐿). It is to be noted that in order to effectively track the interface dynamics, it is essential to have a certain number of grid cells in the interface region. Jacqmin [START_REF] Jacqmin | An energy approach to the continuum surface tension method[END_REF][START_REF] Jacqmin | Calculation of Two-Phase Navier-Stokes Flows Using Phase-Field Modeling[END_REF] recommended to consider a minimum of three nodes in the interface region. • Péclet number (𝑃𝑒): Peclet number in general represents the ratio of convective to diffusive fluxes. In case of phase field modelling, the diffusion is controlled by the mobility (𝑀 0 ) following which 𝑃𝑒 is defined by (for free energy expression given by (2.2.25)),

𝑃𝑒 = 𝑈𝐿 𝑀 0 𝛽 ̅ (3.4.8)
where, 𝑈 is the characteristic velocity while 𝛽 ̅ is the coefficient of bulk free energy. Now, in order to understand the separation mechanism, let us consider a mixture which is non-homogeneous 𝑤. 𝑟. 𝑡 a parameter (𝜙) following which the double well free energy function is defined by (2.2.25). The equilibrium composition corresponds to potential (𝜂) equal to zero which thereby yields an equilibrium composition, 𝜙 = 0 and 𝜙 = 1. It is also well known that the negative curvature of the free energy, 𝜕 2 𝑓 0 𝜕𝜙 2 < 0, corresponds to an unstable region. For the bulk-free energy considered in the present case, this is schematically shown in Fig. 3-11. is considered to be filled with a binary mixture of two phases represented by 𝜙 = 0 and 𝜙 = 1. At 𝑡 = 0, the initial configuration corresponds to an unstable region, 𝜙 0 = 0.5 with small random perturbation.

Thus, the initial condition is represented by, 𝜙((𝑥, 𝑦),0) = 𝜙 0 + 𝜙 𝑟𝑎𝑛𝑑 (3.4.9) where 𝜙 𝑟𝑎𝑛𝑑 is a random distribution between -0.5 × 10 -2 and 0.5 × 10 -2 . Fig. 3-12 shows the evolution and separation of the mixture into binary phases at various time instances.

Chapter Summary

The numerical methodology to solve the mathematical model developed in CHAPTER 2 is described in this chapter followed by validation studies. For a single-phase fluid system, the results are validated by investigation of propagation of acoustic wave in water while for phase-field model, a classical problem of spinodal decomposition is investigated. In addition, the numerical scheme used in the current work is investigated for error analysis and based on the error propagation equation, unusual behavior of low and high Courant numbers affecting the solution is explained.

In the next chapter, the model is subsequently used to analyze the thermo-vibrational instabilities and some new findings are highlighted.

supported by the numerical studies of Jounet et al. [START_REF] Jounet | Low-frequency vibrations in a near-critical fluid[END_REF] wherein the fluid was observed to oscillate homogenously at low frequencies with greater amplitude than in a normal gas. This was explained by virtue of strong thermo-mechanical coupling occurring in the thermal boundary layer. A further extension to the two-dimensional system revealed that the low frequency regime was affected only in the presence isothermal longitudinal boundaries. This was explained by virtue of the formation of longitudinal TBL following which a strong thermo-mechanical coupling provokes a curvature in the streamlines near the walls. When subjected to simultaneous thermal quench and mechanical vibration, evidences of instabilities have been cited by various experimental studies [START_REF] Amiroudine | Thermovibrational instability in supercritical fluids under weightlessness[END_REF][START_REF] Wunenburger | Frozen wave induced by high frequency horizontal vibrations on a CO2 liquid-gas interface near the critical point[END_REF][START_REF] Garrabos | Thermoconvectional phenomena induced by vibrations in supercritical SF6 under weightlessness[END_REF][START_REF] Daniel A Beysens | Effect of oscillatory accelerations on two-phase fluids[END_REF].

One such work has been reported by Amiroudine and Beysens [START_REF] Amiroudine | Thermovibrational instability in supercritical fluids under weightlessness[END_REF] where finger like instabilities were observed when a cell with a sample of hydrogen near its critical density was quenched from 𝑇 𝑐 + 3𝑚𝐾 to 𝑇 𝑐 + 2𝑚𝐾 (𝑇 𝑐 being the critical temperature) and simultaneously subjected to vibrations in weightlessness conditions, as shown in Fig. 4-1. The range of vibrational amplitudes and frequencies considered were higher than the previous work of Amiroudine and Beysens [START_REF] Amiroudine | Thermovibrational instability in supercritical fluids under weightlessness[END_REF] (amplitude varied from 0.05 to 5 times the cell size while frequencies varied between 2.78 𝐻𝑧 and 25 𝐻𝑧). Finally, they presented the stability plot which described the type of instability to expect under given thermo-vibrational parameters.

67

In addition to above described instabilities in zero-gravity conditions, instabilities were also reported for the case of 2 miscible fluids when subjected to vibration parallel to the interface [START_REF] Gaponenko | Dynamics of the interface between miscible liquids subjected to horizontal vibration[END_REF][START_REF] Shevtsova | Two-scale wave patterns on a periodically excited miscible liquid-liquid interface[END_REF]. Shevtsova et al. [START_REF] Shevtsova | Two-scale wave patterns on a periodically excited miscible liquid-liquid interface[END_REF] showed that for amplitudes greater than cell dimensions, secondary instabilities (Faraday instabilities) were observed over frozen wave instabilities. One of the primary assumption in the aforementioned numerical studies [START_REF] Amiroudine | Thermovibrational instability in supercritical fluids under weightlessness[END_REF][START_REF] Gandikota | Rayleigh and parametric thermo-vibrational instabilities in supercritical fluids under weightlessness[END_REF] in SCFs is the linear equation of state for the calculation of density. This assumption not only limits a higher temperature perturbation but also circumvents an accurate analysis very close to the critical point.

The latter is due to the fact that the thermo-vibrational phenomenon is primarily governed by large density gradients in the TBL. Thus, a precise density calculation is essential to ensure that all the physical phenomena are accurately captured by the model. Since the model developed in §2.1 directly evaluates density from mass conservation, it thereby permits a more rigorous investigation of these instabilities, in terms of higher quench rates and decreasing proximity to the critical point (results as close as 𝑇 𝑖 -𝑇 𝑐 = 5𝑚𝐾 have been obtained). In the following section, thermovibrational instabilities are presented as investigated using the developed model in model §2. 𝑚𝐾 and 10% to 50%, respectively.

Problem Description

When subjecting SCFs to simultaneous quench and vibration, the type of instabilities will be governed by the relative direction between the vibrational acceleration and the temperature gradient. These can be categorized as Rayleigh-vibrational (RV) or parametric when the relative direction is normal or parallel, respectively. The nomenclature of these instabilities is attributed to their resemblance to Rayleigh-Bénard convection and parametric (Faraday) instabilities [START_REF] Gandikota | Rayleigh and parametric thermo-vibrational instabilities in supercritical fluids under weightlessness[END_REF],

respectively, as will be discussed in the following sections. In order to gain insight into the physical mechanism leading to each instability, it is essential to analyze each of them separately. This is attained by imposing appropriate boundary conditions SCFs. The appearance of these finger like structures correspond to Rayleigh-vibrational instabilities which will be discussed in detail in §4.4. A further verification is performed by comparing the wavelength (distance between adjacent fingers, 𝜆 , as marked in Fig. 4-4) of the Rayleigh-vibrational instabilities for various conditions.

However, unlike the previous case, the percentage error in experimental and numerical wavelengths is evaluated here. Fig. 4-5 compares wavelengths for different proximities to the critical point (for same temperature quench 𝛿𝑇 = 400 𝑚𝐾 ). As can be observed from Fig. 4-4, the wavelengths are not uniform along the length. The wavelength therefore obtained from experimental data is the average wavelength (dividing cell dimension by the number of fingers) and a similar procedure is adopted in case of numerical investigations. It can be observed that the wavelengths evaluated from the current numerical study show a similar trend and the error decreases from 29 % to 10 % on moving away from the critical point. This is attributed to significant variations in properties on approaching the critical point which can affect the behavior both in experiments and numerical simulation causing a higher discrepancy. Further, unlike the experimental conditions, where the wall temperature is changed slowly, a sudden quench in numerical study also adds up to the discrepancy. A noticeable observation pertains to the decrease in the wavelength on approaching the critical point which is explained in §4.4.3. velocity difference (∆𝑽) due to the difference in density, ∆𝜌, between any two layers (bulk and TBL region in the present case) can be described as ∆𝑽 = ( ∆𝜌 𝜌 𝑀

)𝐴𝜔, where 𝜌 𝑀 is the mean density.

This results in a pressure difference ∆𝑃 = 𝜌 𝑀 ( ∆𝜌 𝜌 𝑀 ) 2 𝐴 2 𝜔 2 . The destabilization due to this Bernoulli-like pressure difference between different fluid layers is stabilized by the effect of diffusive forces (viscous and thermal). However, when these weakening forces overcome the stabilizing effect, the TBL is destabilized resulting in the appearance of finger-like structures as shown in Fig. 456. Thus, the appearance of these finger-like structures corresponds to a minimum of energy, considering both potential and kinetic energy. Though the primary effect leading to these instabilities is attributed to the shear between fluid layers, which therefore has Kelvin-Helmholtz type origin (the term "Kelvin-Helmholtz instability" is usually associated with two immiscible fluids whereas in the present context only a single fluid is considered), as suggested in literature [START_REF] Gershuni | Thermal Vibrational Convection[END_REF], these are termed as "Rayleigh-vibrational instability."

A further motivation behind this terminology is ascribed to the similarity in the mechanism which leads to instability in the Rayleigh-Bénard cell configuration. Consider a fluid contained in a cell which is maintained between two horizontal plates separated by a distance 𝐿 and a constant temperature difference 𝛿𝑇, the lower plate being at a higher temperature. A fluid particle which lies between two plates separated by L will rise once the convective time across it is lower as compared to the diffusion time, 𝑖. 𝑒. 2𝜐𝐷 𝑇 [START_REF] Gershuni | Thermal Vibrational Convection[END_REF],

where symbols have their usual meaning as described before. Here, 𝐿 is the characteristic length scale across which diffusion takes place. Since in supercritical fluids, the diffusion length scale is very small as compared to geometrical length scale (𝐿), it is reasonable to consider the thickness of the TBL (𝛿 𝑇𝐵𝐿 ) as the characteristic length scale. Thus, the stability criterion can be defined in terms of a dimensionless number defined as, 𝑅𝑎 𝑣 = (𝛽 𝑃 𝛿𝑇 𝐴 𝜔 𝛿 𝑇𝐵𝐿 ) 2

2𝜐𝐷 𝑇 [START_REF] Amiroudine | Thermovibrational instability in supercritical fluids under weightlessness[END_REF]. This includes the effect of various thermo-physical properties attributed to (i) the distance from the critical point, (ii) quench conditions (𝛿𝑇) and (iii) vibrational parameters (frequency and amplitude of vibration).

The effect of these parameters on Rayleigh-vibrational instabilities is presented below.

Critical Amplitude for Rayleigh vibrational instabilities

As described in the previous section, since there exists a competition between shear and diffusive forces, it is therefore possible to define a critical amplitude (𝐴 𝑐𝑟 ) of vibration for a fixed frequency beyond which these instabilities may be observed. In order to ascertain the critical amplitude, the simulation was run till 𝑡 = 50𝑠 and the smallest amplitude at which instability was detected (waviness in the TBL) was chosen as the critical amplitude for that frequency and quench rate. It can be observed that the critical amplitude decreases with increase in frequency and quench rate. This can be explained as follows. A higher frequency entails a higher acceleration and thus a larger difference in the inertial forces between layers of different densities when compared to a lower frequency. Similarly, a higher quench rate will result in a higher temperature gradient thereby causing a higher density variation leading to the same conclusion, for the same vibration parameters (𝑖. 𝑒. frequency and amplitude). This implies that for the same stabilizing action provided by the diffusive forces, the desired threshold to overcome this will be attained at a lower amplitude for a higher frequency or quench rate as illustrated in Fig. 4-8(a). A significant remark not to be missed here is that physical meaningful results have been obtained as close as 𝑇 𝑖 -𝑇 𝑐 = 5 𝑚𝐾 with quench rate as high as 50% illustrating the strength and capability of the model to capture the behavior of a highly compressible system with large property variations.

Another important aspect pertaining to the critical amplitude is the proximity to the critical point. Furthermore, a higher value of 𝐴 𝑐𝑟 is observed for a lower quench temperature which can be ascribed to a lower density gradient resulting from it. Consequently, a lower difference in inertial forces for same vibrational parameters necessitates a higher amplitude to overcome the stabilizing action of diffusive forces leading to the observed behavior in Fig. 4-8(b).

Analysis of wavelengths as a function of various parameters

The Rayleigh-vibrational instabilities described above lead to the appearance of finger like structures, as reported in literature [START_REF] Amiroudine | Thermovibrational instability in supercritical fluids under weightlessness[END_REF][START_REF] Gandikota | Rayleigh and parametric thermo-vibrational instabilities in supercritical fluids under weightlessness[END_REF] and are shown in Fig. 456. These fingers can be characterized by the distance between them 𝑖. 𝑒. wavelength of the instabilities. It has been mentioned in [START_REF] Amiroudine | Thermovibrational instability in supercritical fluids under weightlessness[END_REF][START_REF] Gandikota | Rayleigh and parametric thermo-vibrational instabilities in supercritical fluids under weightlessness[END_REF] that wavelengths do not change over a period of time. However, it is observed that after the onset of instabilities, the wavelength remains constant for a certain period It can be seen that initially the distance between fingers remains constant for certain duration of time (Fig. 4-9(a-b)). However, after the onset, these not only grow into the bulk but also become thicker in size due to diffusion enhanced by vibrational acceleration. When these are close enough to each other, the adjacent fingers merge with each other and hence results in the reduction of the number of fingers ( 𝑖. 𝑒. increase in wavelength) as highlighted in Fig. 45678910). The merging process can be further explained by virtue of interaction of vortices which are formed when fingers emerge at the onset of these instabilities. This is illustrated in Fig. 4-10(a-d) wherein velocity vectors are plotted for the merging process of 2 fingers highlighted in Fig. 456789). For the sake of clarity in illustration, the velocity vectors near the wall have not been shown. It can be observed that the velocity field (vortices) corresponding to two individual and distinct fingers (Fig. 4-10(a))

begin to interact with each other as shown in Fig. 4-10(b). In addition to diffusion, this is also attributed to vibrational acceleration in horizontal direction which causes an increase in the horizontal velocity as fingers grow into the bulk. Subsequently, the velocity vectors tend to overlap and finally merge to form a single finger (Fig. 45678910)). The merging of fingers, in general, will depend on amplitude of vibration when other thermal and mechanical parameters are held constant.

However, in the current work, the wavelength at the onset of instabilities have been measured

at the middle of the cell in order to circumvent the influence of boundaries. A higher quench rate implies a higher temperature gradient and thus a larger density difference. Thus, for the same acceleration, a higher Bernoulli-like pressure difference will exist between the layers for a higher quench rate. This implies that the stabilizing viscous forces will be weakened at more local sites along the TBL in order to minimize the energy (potential and kinetic).

• Effect of acceleration

Subsequently, more fingers (or a lower wavelength) with a higher quench rate is observed for the same acceleration. An anomalous behavior in this context can be seen when for 𝑇 𝑖 -𝑇 𝑐 = 5 𝑚𝐾, wherein at the same higher acceleration, the wavelength increases with the quench rate. While this is clearly visible for 𝑞 = 30%, it can be visualized that when the line for 𝑞 = 50% is extrapolated, a similar behavior will be observed. A further investigation reveals that a similar behavior can also be attained with higher acceleration (greater than the ones considered in Fig. 4-11(a)) for 𝑇 𝑖 -𝑇 𝑐 = 500 𝑚𝐾, for 𝑞 = 10% and 𝑞 = 50% as demonstrated in Fig. 456789101112.

As mentioned previously, a higher Bernoulli-like pressure difference attributed to either a higher acceleration or a larger density variation will result in a higher number of plausible local sites where the TBL will be destabilized. In addition, when these fingers grow into the bulk, they have been found to merge as depicted in Fig. 456789and Fig. 45678910. Now, considering the same analogy, as the acceleration is increased, the wavelength at the onset becomes so small (due to the increase in probable sites of onset of instabilities) that as the fingers physically appear and protrude into the bulk, they merge into one another. As a result, a decrease in the effective number of fingers and thus a higher wavelength is observed. In order to further support this argument, wavelengths for 𝑞 = 10% and 𝑇 𝑖 -𝑇 𝑐 = 500 𝑚𝐾 for cell size ℎ = 7𝑚𝑚 and ℎ = 14 𝑚𝑚 are compared in Fig. 4-13. It can be seen that the change from decreasing trend of wavelength is attained at a higher acceleration for cell size of ℎ = 14 𝑚𝑚. This is coherent with the above explanation because for a longer cell, a wider span is available for the appearance of fingers before these merge into each other leading to a higher wavelength. Thus, for a longer cell, a higher acceleration will be required so that plausible number of local sites for onset of instabilities become nearly the same as the span of TBL thereby ensuing coalescence of fingers when they appear as instabilities. A similar explanation holds for the observations in Fig. 4-11(b) wherein such an anomalous behavior is observed even at a low acceleration. This can be related to the diverging behavior of thermal expansion which eventually leads to large pressure difference even at low acceleration and hence the observed behavior as in Fig. 4-11(b). It is to be remarked here that a higher wavelength (which implies lower number of fingers) is also observed for a cell with larger dimensions (Fig. 45678910111213). This is due to the fact that for a constant vibrational acceleration, the local sites along the TBL prone to instabilities will be nearly the same. Hence, for a longer cell, the fingers will be separated by a larger distance resulting in a higher wavelength as observed in Fig. 45678910111213. • Effect of same quench temperature The experimental observations in Fig. 4-5 illustrated the effect of proximity to the critical point for the same vibration parameters and quench value (𝛿𝑇). An additional analysis to confirm this behavior with proximity to the critical point is shown in Fig. 4-14 where wavelengths for two different quench temperatures (𝛿𝑇 = 100 𝑚𝐾 and 5 𝑚𝐾) are compared for different proximities to the critical point value for various accelerations. It can be clearly seen that for a constant acceleration, the wavelength of instability decreases when approaching the critical point. This behavior is attributed to higher Bernoulli-like pressure difference when approaching the critical point (due to larger density variations) which results in more fingers (and hence lower wavelengths) as explained previously.

Rayleigh vibrational number as a distance from critical point

The Rayleigh vibrational number was introduced in §4.4.1 and has been known to describe the onset criteria for Rayleigh-vibrational instabilities as a function of various thermo-physical and vibrational parameters [START_REF] Amiroudine | Thermovibrational instability in supercritical fluids under weightlessness[END_REF]. One of the important parameter in 𝑅𝑎 𝑣 is the thickness of the thermal boundary layer (𝛿 𝑇𝐵𝐿 ) which changes as a function of time. Thus, in order to evaluate the critical 𝑅𝑎 𝑣,𝑐 , the time at which the TBL becomes unstable, described by 𝑡 𝑐 , is noted from the numerical simulations following which the critical thickness of the TBL is evaluated using 𝛿 𝑇𝐵𝐿,𝑐 = 2𝜋√𝐷 𝑇 𝑡 𝑐 . Subsequently, the critical Rayleigh vibrational number 𝑅𝑎 𝑣,𝑐 can be calculated using the expression, 𝑅𝑎 𝑣,𝑐 = (𝛽 𝑃 𝛿𝑇 𝐴 𝜔 𝛿 𝑇𝐵𝐿,𝑐 ) 2 2𝜐𝐷 𝑇 as defined in §4.4.1. Amiroudine and Beysens [START_REF] Amiroudine | Thermovibrational instability in supercritical fluids under weightlessness[END_REF] evaluated 𝑅𝑎 𝑣,𝑐 for a fixed frequency and amplitude as a distance from the critical point and found it to increase on approaching the critical point with power exponents of -0.86 and -0.31 when close and far from the critical point, respectively. Similarly, Gandikota et al. [START_REF] Gandikota | Rayleigh and parametric thermo-vibrational instabilities in supercritical fluids under weightlessness[END_REF], analyzed 𝑅𝑎 𝑣,𝑐 for a fixed frequency and amplitude over a wider range of amplitudes than Amiroudine and Beysens [START_REF] Amiroudine | Thermovibrational instability in supercritical fluids under weightlessness[END_REF] and observed a similar behavior.

In the present work, the 𝑅𝑎 𝑣,𝑐 is analyzed only as a function of a single parameter, vibration frequency. The amplitude of vibration was varied from 0.15ℎ to 1.5ℎ . In addition, the quench temperature was also varied for each case with quench rate was varying from 𝑞 = 5 % to 30%. This is motivated by the fact that it is more reasonable to draw out the behavior of 𝑅𝑎 𝑣,𝑐 based on a single parameter rather than fixing more parameters. Thus, based on different combinations of amplitude, proximity to the critical point and quench rate (and thus quench temperature), the critical Rayleigh-vibrational number was calculated and is shown in Fig. 4-15 for two different frequencies, 𝑓 = 10 𝐻𝑧 and 30 𝐻𝑧. It can be seen that 𝑅𝑎 𝑣,𝑐 increases on approaching the critical point. This is mainly attributed to two factors, the vanishing thermal diffusivity (𝐷 𝑇 ) on approaching the critical point which increases the diffusion time scale and secondly the diverging thermal expansion (𝛽 𝑃 ) which causes a higher Bernoulli-like pressure difference and thus a lower convective time scale. Furthermore, it can be seen that, similar to the studies of Amiroudine and Beysens [START_REF] Amiroudine | Thermovibrational instability in supercritical fluids under weightlessness[END_REF], 𝑅𝑎 𝑣,𝑐 can be described by a power law fit as 𝑅𝑎 𝑣,𝑐 = 𝑅𝑎 1,2 𝜀 𝑚 1,2 with two different slopes when close and far (𝑇 -𝑇 𝑐 > 100 𝑚𝐾) to the critical point. The increase in 𝑅𝑎 𝑣,𝑐 on approaching the critical point illustrates the same behavior, 𝑖. 𝑒. a higher exponent (𝑚 1 ≈ -1.7) when near and a lower (𝑚 1 ≈ -1) when far from the critical point, respectively (exponents obtained by power law fit 𝑅 2 ~ 0.92-0.99). A lower value of exponents in the work of Amiroudine and Beysens [START_REF] Amiroudine | Thermovibrational instability in supercritical fluids under weightlessness[END_REF] can be attributed to lower quench rates considered due to the assumption of a linear equation of state for density calculation. However, a higher quench rate, as favored by the proposed mathematical model, makes it possible to account for higher variations especially when approaching the critical point. As a result, it is reasonable to expect a faster increase in the critical Rayleigh vibrational number as described above.

Parametric instabilities

It is well established that a two-phase fluid (miscible or immiscible) exhibits Faraday instabilities when subjected to vibrations perpendicular to its interface [START_REF] Kumar | Parametric instability of the interface between two fluids[END_REF] . On a similar basis, when the SCF close to its critical point is simultaneously quenched and subjected to vibrations of sufficiently high amplitudes in the direction parallel to the temperature gradient, the TBL becomes unstable leading to the appearance of waves or finger like structures. In order to understand the physical phenomenon inducing these instabilities, consider the schematic as shown in Fig. 4-17(a) where the SCF, initially at temperature 𝑇 𝑖 , is quenched on the left wall (𝐴𝐷). The density of the fluid in the TBL will be significantly higher than that of the bulk fluid due to lower temperature therein. Thus, if considering the fluids in the bulk and the TBL as two different fluids (owing to the significant difference in their properties), an equivalent configuration can be represented as shown in Fig. 4-17(b). This resembles to an arrangement similar to that of parametric (Faraday) instability wherein two immiscible fluids are acted upon by vibration under zero-gravity. Thus, when the natural frequency of the system matches with the external frequency, wavelike patterns as shown in Fig. 4-16 are observed which means that in one cycle of vibration, the TBL near each quenched wall is destabilized only once thereby exhibiting a harmonic behavior similar to that of parametric instabilities. This is can be observed in Fig. 4-16 where for one time-period of vibration (𝑓 = 10 𝐻𝑧, 𝑇 = 0.1 𝑠), the instabilities are observed once on each side. The appearance is hence found to depend on the vibration parameter, frequency, following which these are termed as parametric instabilities. Gandikota et al. [START_REF] Gandikota | Rayleigh and parametric thermo-vibrational instabilities in supercritical fluids under weightlessness[END_REF] described these instabilities as sub-harmonic in nature. This is due to the fact that their characterization pertains to the observance of these instabilities once every half time-period of vibration. However, it would be more appropriate to ascertain the periodicity based on their behavior at a particular TBL (or interface like conditions) and not considering the cell as a whole.

In continuation with the above explanation, a further insight into the mechanism leading to the onset of parametric instabilities is described herein. Gandikota et al. [START_REF] Gandikota | Rayleigh and parametric thermo-vibrational instabilities in supercritical fluids under weightlessness[END_REF] ascribed these instabilities to formation of vortices though the exact mechanism leading to these vortices is not explained. Furthermore, it has also been mentioned in [START_REF] Gandikota | Rayleigh and parametric thermo-vibrational instabilities in supercritical fluids under weightlessness[END_REF] that these instabilities are observed when the direction of vibration is anti-parallel to the temperature gradient and thus leads to Rayleigh-Taylor like configuration. While the results in the present work were found to be in coherence with former remark, the observations differed regarding the Rayleigh-Taylor like configuration. For instance, the instabilities are observed near the wall 𝐴𝐷 (temperature gradient towards positive 𝑥 -direction) when the cell moves in negative 𝑥 -direction. This is supported by Fig. 4 ). An equivalent configuration is shown in Fig. 4-17(c) wherein a lighter fluid rests atop a heavier fluid with downwards acting acceleration which is a stable configuration 𝑤. 𝑟. 𝑡 to Rayleigh-Taylor instability. Thus, these instabilities are described to resemble to parametric-like instabilities. In the following section, an in-depth explanation is provided leading to the appearance of these instabilities. In addition, it is also described that why the vortices formed near the corner lead to non-uniform onset of parametric instabilities (as in Fig. 4-16), 𝑖. 𝑒. waves initially appear near the horizontal walls and then propagate towards the middle of the TBL.

Mechanism of Parametric instabilities

As illustrated in Fig. 4-17(a-b), our configuration is equivalent to parametric instabilities in case of two immiscible fluids under zero-gravity. However, in the present case, as there is only a single fluid, no interface exists between fluids of different densities. For the sake of clarity and ease in understanding, let us call the boundary separating the bulk and the TBL as a 'virtual-interface' (the bulk and the TBL region may be considered as two different fluids as described in the previous section). Thus, the instabilities will be observed when the virtual-interface is adequately perturbed.

The important question which arises here is, how is this perturbation obtained? In order to explain this, let us consider a case with left wall (𝐴𝐷) quenched and the cell vibrating along 𝑥 -direction.

Under the same vibrational acceleration, the inertia associated with lighter fluid (bulk region) will be lower as compared to the heavier fluid (TBL region). When a fluid element from the bulk impacts the element of the TBL (or vice-versa), the transfer of momentum therefore yields different velocities in both the regions near the virtual-interface. Since the overall momentum is conserved, a small velocity field in 𝑦 -direction will be induced in the fluid-elements near the virtual-interface.

Physically, the conditions near the virtual-interface can be visualized as several objects of same mass adjacent to each other (fluid elements of the bulk region) impacting objects in same arrangement but of different mass (fluid elements in the TBL region). As the objects are free to slide among themselves, upon impacting objects of different mass, a velocity will also be induced in the direction normal (𝑦 -direction) to their initial motion (𝑥 -direction) due to momentum conservation. Thus, in the present context, the net velocity field of the fluid elements near the virtual interface will be inclined 𝑤. 𝑟. 𝑡 the bulk velocity field (which is primarily in 𝑥 -direction). This will thereby produce rotational or vortex like characteristics which act as a perturbation but they are stabilized by shear action of viscous forces. With the growth in the thickness of the TBL, the vortex strength increases due to higher inertia associated with thicker TBL resulting in stronger perturbation. When the perturbations are able to overcome the stabilizing action of viscous forces, instabilities are observed as shown in Fig. 45678910111213141516.

It is to be mentioned that even though viscous effects may be very small, as mentioned in Gandikota et al. [START_REF] Gandikota | Rayleigh and parametric thermo-vibrational instabilities in supercritical fluids under weightlessness[END_REF], their action is sufficient enough to stabilize the small perturbations induced by the velocity field. It is the formation of these vortices as described in [START_REF] Gandikota | Rayleigh and parametric thermo-vibrational instabilities in supercritical fluids under weightlessness[END_REF] that are responsible for parametric instabilities. The above explanation therefore suggests that parametric instabilities are possible even between infinite parallel plates. In order to verify this, the system was investigated with symmetric horizontal wall conditions imposed at walls 𝐴𝐵 and 𝐶𝐷 (infinite parallel plates like conditions) for 𝑇 𝑖 -𝑇 𝑐 = 500 𝑚𝐾, 𝛿𝑇 = 50 𝑚𝐾, 𝑓 = 20 𝐻𝑧 and 𝐴 = 6 𝑚𝑚. The results are shown in Fig. 45678910111213141516171819) for symmetric conditions and in Fig. 45678910111213141516171819) for wall conditions at various time instances. The illustrations have been made only near the left wall (𝐴𝐷) for the sake of brevity.

While instabilities observed for symmetrical conditions substantiates our explanation, it is worth noticing that unlike wall conditions, the instabilities appear uniformly across the entire span of TBL for symmetrical conditions. This difference can be ascribed to the vortices formed near the corners which can be explained as follows. The thermal quench on the vertical walls leads to the formation of the TBL while vibration causes the fluid to move from the TBL towards the bulk and vice-versa.

Critical Amplitude for parametric instabilities

A judicious investigation into the phenomenon leading to the onset of parametric instabilities has been described in the previous section. It can therefore be concluded that for a fixed proximity to the critical point and a given frequency, there exists a critical amplitude (which represents the critical acceleration) for which the perturbation will grow to overcome the stabilizing effect of viscous forces. Fig. 4-21(a) illustrates the critical amplitude for different proximities to the critical point (𝑇 𝑖 -𝑇 𝑐 = 2000 𝑚𝐾, 500 𝑚𝐾 and 50 𝑚𝐾 ) as a function of frequency and quench rate. It can be observed that critical amplitude decreases with increase in frequency and quench rate. While a higher frequency entails a higher inertial force, a higher quench rate implies a heavier fluid (alongwith larger density gradient) and thus a larger difference in associated momentum. Thus, the threshold to overcome the stabilizing viscous forces is attained with lower amplitudes for the both aforementioned cases. In order to compare the effect of proximity to the critical point, as explained in case of RV instabilities, it is more reasonable to compare the behavior for same quench conditions. This is shown in Fig. 4

Effect of cell size on critical amplitude

The onset of parametric instabilities as previously described is primarily due to the perturbation in the velocity field near the transition region of the bulk and the TBL. The perturbation will be higher if larger momentum transfer occurs from the bulk fluid to the TBL and vice-versa. This implies that the change in cell dimensions along the direction of vibration can have effect on the critical amplitude. The argument is supported from Fig. 4-22 which illustrates the critical amplitude for 𝑇 𝑖 -𝑇 𝑐 = 500 𝑚𝐾 and 𝑇 𝑖 -𝑇 𝑐 = 20 𝑚𝐾 for different cell sizes for quench rate of 𝑞 = 10% and 𝑓 = 20 𝐻𝑧. It can be seen that with all other parameters held constant, the critical amplitude decreases with the increase in cell dimension which can be explicated as follows. The induced perturbation will be of higher strength by virtue of higher momentum associated with higher volume (or mass) in a cell of higher dimensions thereby leading to a lower critical amplitude for onset of instabilities. It is also worth noticing that the exponent of the decrease of this critical amplitude is lower (-0.22 for 𝑇 𝑖 -𝑇 𝑐 = 20 𝑚𝐾 and -0.45 𝑇 𝑖 -𝑇 𝑐 = 20 𝑚𝐾 ) (based on a power law fit) when approaching the critical point. This physically implies that the critical amplitude is not significantly affected by the change in cell size when close to the critical point. This can be explained based on the diverging isothermal compressibility (𝜒 𝑇 ) when approaching the critical point. When the fluid from the center of the cell moves towards the walls, the transfer of momentum is highly affected by the compressibility of the fluid. Thus, if it was an incompressible medium, it will act as a solid body and transfer the same momentum. However, with an increasing isothermal compressibility (𝜒 𝑇 ), a significant portion of the momentum is dampened out thereby transmitting a lower momentum. This implies that the overall effect of cell volume (or mass/momentum) would be lower when fluid is more compressible (the aspect of damping will be described in more detail in §CHAPTER 5).

Consequently, the difference between the critical amplitude for different cell sizes does not vary too much on approaching critical point and thus a lower exponent as mentioned above is observed.

Wavelength in parametric instabilities

Gandikota et al. [START_REF] Gandikota | Rayleigh and parametric thermo-vibrational instabilities in supercritical fluids under weightlessness[END_REF] calculated the wavelengths in parametric instabilities and tried to verify the dispersion relation for immiscible fluids given by [START_REF] Kumar | Parametric instability of the interface between two fluids[END_REF] ,

𝜔 2 = (𝜌 1 -𝜌 2 )𝑔𝑘+𝑘 3 𝜎 𝜌 1 +𝜌 2 (4.5.1)
Here 𝜌 1 and 𝜌 2 correspond to density of the two fluids while 𝜎 and 𝑘 represent surface tension and wavenumber, respectively.

95

A similar analysis is presented here with higher quench rates for a given proximity to the critical point (𝑇 𝑖 -𝑇 𝑐 = 500 𝑚𝐾). In the work of Gandikota et al. [START_REF] Gandikota | Rayleigh and parametric thermo-vibrational instabilities in supercritical fluids under weightlessness[END_REF], wavelength was observed to decrease with 𝐴𝜔 2 with an exponent of -0.35 based on a power law fit. A nearly similar exponent in observed in the current work for 𝑞 = 10% (𝜆 ~(𝐴𝜔 2 ) -0.36 ) while a higher exponent is observed (𝜆 ~(𝐴𝜔 2 ) -0.47 ) for 𝑞 = 50%. This can be ascribed to a larger variations in density causing dominance of non-linear effects thereby leading to a higher exponent (-0.47) than given by (4.5.1), which is obtained using linear stability analysis.

Stability Analysis

One of the primary motivations to analyze SCF under thermo-vibrational instabilities has been to gain insights into the type of instabilities one may expect at different frequencies and amplitudes of vibration. While these may be used in certain processes to our advantage, for example mixing, in others, it may be desired to isolate the system to these instabilities. For example, in case of a heterogeneous chemical reaction taking place on the walls of a container, large density perturbations near the boundaries can cause an enhancement of mass transfer rate. In the presence of such a reaction, in case of corrosive SCF stored in a closed container and acted upon by vibrations, this can accelerate the corrosion rate. Gandikota et al. [START_REF] Gandikota | Rayleigh and parametric thermo-vibrational instabilities in supercritical fluids under weightlessness[END_REF] presented a plot limited by lower quench rate conditions (10% due to linear state equation). As the current model makes it possible to analyze system with a higher quench rate, a similar plot, but in 3D as a function of different quench rates and frequencies, for two different proximities to the critical point 𝑇 𝑖 -𝑇 𝑐 = 500 𝑚𝐾 and 𝑇 𝑖 -𝑇 𝑐 = 50 𝑚𝐾 is illustrated in Fig. 4-24.

The plot thus shows that for amplitudes below the blue plane, the system will be stable 𝑤. 𝑟. 𝑡 to both the described instabilities. Similarly, for amplitude bounded by blue and yellow planes, Rayleigh-vibrational instabilities will be first observed while above the yellow plane, both these instabilities will occur simultaneously. Thus, it can be observed that higher amplitudes are required for the onset of parametric instabilities in both the cases. While it can be intuitively ascertained that it is easy to destabilize the fluid layers in parallel due to shear, it can be explained more judiciously based on the mechanism causing both these instabilities. The Bernoulli-like pressure in Rayleighvibrational instabilities was observed to grow with square of the velocity difference while in parametric, owing to momentum transfer it is likely to follow a linear relation. Furthermore, since the momentum transfer occurs between two highly compressible fluids, some of the momentum is also expended in compressing the fluid whereas no such loss can be ascertained in case of Rayleighvibrational instabilities thereby leading to the observed behavior. It can also be seen that when close to the critical point, the difference between the critical amplitude of both these instabilities increase with quench rate. This can be explained by virtue of the fact that Rayleigh-vibrational instabilities are more sensitive to the difference in inertial forces, 𝑖. 𝑒. density variations and thus quench rates, due to the instability arising from fluid layers in shear/parallel to each other.

Chapter Summary

The two types of thermo-vibrational instabilities, defined by the relative direction between the temperature gradient and vibration (Rayleigh-vibrational when perpendicular and parametric when parallel) are described. The physical mechanism causing these instabilities is discussed in detail. It is found that for both of these instabilities, it is possible to define a critical amplitude for their onset which is observed to decrease with the increase in quench rate and frequency of vibration. Further, it also decreases on approaching the critical point. These observations are attributed to a higher difference in the inertial forces (higher perturbation causing vortex formation) at a higher value of both these factors, both in Rayleigh-vibrational and parametric instabilities. An interesting observation pertains to the increase in wavelength at higher acceleration in case of Rayleigh-vibrational instabilities while in case of parametric instabilities, the critical amplitude is found to depend on size of cell. Finally, a stability plot highlighting critical regions in terms of thermo-vibrational parameters (quench rate, frequency and critical amplitude) is presented.

In the next chapter, we discuss two intriguing phenomena observed along with Rayleighvibrational instabilities, firstly, the drop of the temperature in the fluid domain below the imposed temperature at the boundary and the second being the see-saw motion of the thermal boundary layer.

CHAPTER 5. HEAT-SINK AND SEE-SAW EFFECTS IN RAYLEIGH-VIBRATIONAL INSTABILITIES

The simultaneous action of thermal quench and linear vibrations in supercritical fluids lead to thermo-vibrational instabilities as described in CHAPTER 4. In this chapter, the following intriguing phenomena observed along-with Rayleigh-vibrational instabilities are described and explained:

• When a fluid is subjected to a temperature imposed at the boundary, it is intuitive to ascertain that the temperature in the fluid domain will be constrained by that temperature limit. However, it is observed that upon being quenched, the temperature of the SCF can fall below the imposed condition. These anomalous regions within the bulk are defined as sink-zones as they act like sink for heat (Sharma et al. [START_REF] Sharma | Cooling beyond the boundary value in supercritical fluids under vibration[END_REF])

• The thermal homogenization in SCF is known to be caused by the action of piston effect, 𝑖. 𝑒.

propagation of the pressure wave due to the high compressibility of the SCF resulting in a temperature change. But how will the temperature field be affected if the pressure change is caused by external means, for example vibration? Will thermo-mechanical coupling lead to alteration in thermal behavior? These aspects are described in the second half of the chapter and the details can also be found in Sharma et al. [START_REF] Sharma | See-saw motion of thermal boundary layer under vibrations: An implication of forced piston effect[END_REF]. It is observed that coupling between temperature variations in the bulk due to vibration (termed as forced piston effect (FPE) in vibration for reasons presented therein) and temperature in the TBL causes see-saw motion of the TBL.

Sink Zones: Regions with temperature below the boundary value

Contemporary prior work

The anomalous behavior of thermo-physical properties in SCFs have led to several unique observations such as drop in the temperature of the bulk below the initial temperature on heating.

This unconventional cooling was observed by Beysens et al. [START_REF] Beysens | Heat can cool near-critical fluids[END_REF] in their experimental work. A cell filled with SCF was subjected to a heat pulse with boundary temperatures held constant and at the end of heat pulse, the bulk fluid was observed to cool below the initial temperature after a long thermal transient. This effect was further found to escalate with increasing acceleration of free fall (or gravitational field). Another captivating observation in SCF was reported by Boukari et al. [START_REF] Boukari | Calculation of the dynamics of gravity-induced density profiles near a liquid-vapor critical point[END_REF] wherein they numerically investigated a 1D problem to compare the time scales of density and temperature equilibration in the presence and absence of gravity. A temperature gradient was found to be formed rather quickly in the bulk and which led to fall in fluid temperature below the boundary temperature. A nearly similar observation was also reported by Zappoli et al. [START_REF] Zappoli | Thermoacoustic and buoyancy-driven transport in a square side-heated cavity filled with a near-critical fluid[END_REF] wherein they investigated the mechanism of heat and mass transport in a side-heated square cavity in the presence of gravity. In all these previous works, this surprising phenomenon has been observed but providing a better and more clearer explanation constitutes one of the objectives of the present study.

In the current work, a near similar observation is made when under the action of simultaneous thermal perturbations and mechanical vibration, the temperature of the bulk fluid drops below the imposed boundary temperature when the direction of acceleration is normal to the temperature gradient. This is a primary difference from the study of Boukari et al. [START_REF] Boukari | Calculation of the dynamics of gravity-induced density profiles near a liquid-vapor critical point[END_REF] wherein the direction of acceleration and temperature gradient are parallel to each other due to 1D study.

Another notable difference is that the direction of acceleration (gravity plays an equivalent role of acceleration) remains constant while it changes in the present investigation due to vibration.

Preliminary observations

Rayleigh-vibrational instabilities were introduced and explained in §4. This implies that the temperature in the domain drops below the imposed boundary conditions which is also illustrated in insets of Fig. 5-1(a-b). In order to ensure that this behavior is not attributed to the adiabatic boundary conditions imposed at the walls 𝐵𝐶 and 𝐴𝐷, a similar analysis with all walls (𝐴𝐵 , 𝐵𝐶, 𝐶𝐷 and 𝐴𝐷) (Fig. 4-3) quenched is performed and same observations are made. Further, analysis of the system with finer mesh sizes yield similar results which disregards these observations being attributed to numerical artifacts. However, the results pertaining to these have not been shown for the sake of brevity.

Features of sink-zones

Before moving further some peculiar features of sink-zones as observed from results are highlighted which will help in better elucidating the physical phenomena underlying their formation.

• Sink-zones are also observed when all the walls except top are adiabatic. This is depicted in in Fig. 5-3(a) which shows the contour plots of temperature (highlighting the sink-zones with negative values) when only the top wall is quenched whilst other parameters are same as in Fig. 5-1. However, when the bottom wall is maintained at 𝑇 𝑖 𝑖. 𝑒. isothermal conditions, no such zones are observed as can be seen in Fig. 5-3(b). Thus, it can be remarked that the formation of sink-zones is attributed to the absence of any possible mechanism of heat inflow into the system. This aspect can be explained for isothermal configuration (Fig. 5-3(b)) as follows: the decrease in temperature of the bulk caused by the piston effect leads to the formation of a hot thermal boundary layer near the bottom wall (𝐶𝐷) which leads to heat flow into the system at wall 𝐶𝐷. Thus, no local heat sink-zones are formed in the fluid domain. Thus, the aforementioned conditions form one of the necessary condition for the formation of sink-zones in the bulk. While this condition may seem quite intuitive as any heat flow into the system will preclude the feasibility of temperature falling below the boundary value, it will be shown in subsequent section that this is not the sufficient condition for the appearance of sink-zones. • With increase in time, the region (fluid domain) where the temperature is below the boundary value spans a larger section of the bulk. Thus, the domain of occupancy of sink-zones grows in size with time. In addition, there is a slight variation in their size depending on the direction of vibration 𝑖. 𝑒. sink-zones tend to decrease in size when the fluid moves towards them.

These two aspects can be observed in Fig. 5-4.

• The formation of sink-zone is observed even at far proximities to the critical point, though very late in time (depending on the proximity and for high acceleration). This further supports their formation being attributed to some physical behavior peculiar to SCF which is further illustrated Fig. 5-5(a-b) for various proximities to the critical point. It is to be mentioned that even though there are significant differences in the configuration analyzed by Boukari et al. [START_REF] Boukari | Calculation of the dynamics of gravity-induced density profiles near a liquid-vapor critical point[END_REF], the same case is initially investigated and described for two reasons.

Firstly, a similar qualitative observation will serve as an additional validation of the mathematical model and observed phenomena. Secondly, as will be shown in subsequent section, the physical phenomena leading to the formation of sink-zones in 2D is closely related to the case when the direction of vibration (acceleration) is along temperature gradient, a clear explanation of the same will aid in better understanding of its formation and the results presented herein.

One-dimensional analysis with both walls quenched

We consider a 1D case similar to the work of Boukari et al. [START_REF] Boukari | Calculation of the dynamics of gravity-induced density profiles near a liquid-vapor critical point[END_REF] which is schematically presented in Fig. 5-6(a). The fluid is initially at 20 𝑚𝐾 from the critical point (same as in case of Boukari et al. [START_REF] Boukari | Calculation of the dynamics of gravity-induced density profiles near a liquid-vapor critical point[END_REF] ) and a quench of 𝛿𝑇 = 10𝑚𝐾 is imposed on both the sides. In addition, the fluid is subjected to an acceleration of 100𝑔 in the negative 𝑦 -direction, where 𝑔 is the magnitude of acceleration due to gravity. It is to be mentioned that the motive for selecting this large value of acceleration was to impede the long duration of simulations. A smaller value of acceleration will affect the behavior of the system only quantitatively (in time) and increase unnecessary computation time. Thus, a higher value of acceleration is selected for the current analysis and explanation.

However, in later sections, some results for lower values of acceleration are also presented. In addition to analyzing the case of 100𝑔, a similar case under weightlessness is also analyzed in order to make necessary comparisons essential to explain the results. with the studies of Boukari et al. [START_REF] Boukari | Calculation of the dynamics of gravity-induced density profiles near a liquid-vapor critical point[END_REF] wherein temperature drops below the temperature imposed at the boundary after a certain period of time. The principle leading to this peculiar behavior which will be further used to explain the formation and behavior of sink-zones is now explained. The quench at the top and bottom wall induces a very thin thermal boundary layer (TBL) and the high value of thermal expansion (𝛽 𝑃 ) (contraction in the present case) leads to large density gradients within the TBL (property of near-critical fluids). The large value of thermal expansion at constant pressure (𝛽 𝑃 ) in conjunction with very high compressibility (𝜒 𝑇 ) cause the fluid in the bulk to expand and thereby its temperature to decrease. This is the classical piston effect as can be observed in Fig. 5-6(b) (dotted lines) in zero-g conditions. Now consider the case when acceleration is acting in negative 𝑦 direction (case (i)). The heavier fluid in the top TBL and fluid in the bulk, exerts a force (due to its self-weight) on the fluid beneath it. Consequently, the fluid is compressed with maximum compression being experienced by the fluid near the bottom, while near the top the fluid expands in addition to usual expansion by piston effect. Thus, a high compressibility in the bulk consequently affects the temperature field. Therefore, the overall temperature profile is then a superimposed effect of expansion of fluid due to the piston effect and compression/expansion due to self-weight depending on the location of fluid element from the bottom. This results in temperature near the bottom being higher than in weightlessness condition as can be clearly seen in Fig. 56.

The combined action of compression of bulk fluid by this heavy fluid and action of selfweight causes a faster decrease in temperature near the top as compared to the bottom.

Consequently, while the temperature of the fluid near the top reaches the wall temperature, the fluid near the bottom is yet to attain the imposed wall temperature and thus keeps on expanding. This higher temperature persisting near the bottom has been mentioned as a sharp temperature gradient by Boukari et al. [START_REF] Boukari | Calculation of the dynamics of gravity-induced density profiles near a liquid-vapor critical point[END_REF] in their study. This temperature gradient near the bottom thus causes the fluid to expand thereby decreasing the temperature everywhere. As a result, the temperature near the top wall falls below the boundary value and extends into the bulk fluid with time. This can also be observed in the Fig. 5-6(b) for 𝑡 = 5 𝑠 and 𝑡 = 25 𝑠. As mentioned in [START_REF] Boukari | Calculation of the dynamics of gravity-induced density profiles near a liquid-vapor critical point[END_REF] , the temperature of the fluid will continue to fall below the wall temperature until the fluid in the top TBL starts to expand due to heating from the top wall and thus exerting an opposite effect 𝑖. 𝑒. compressing the bulk fluid.

When both these opposite effects (the expansion and compression of the bulk due to bottom and top boundary layer respectively) become equal, the temperature would stop to fall below the boundary value. Subsequently, it will try to return to the imposed values as imposed by boundary conditions. It is to be mentioned that owing to 1D consideration, convection has been ignored in the present analysis.

The above description thus illustrates the possibility that the fluid temperature can fall below the boundary value. While in the above described case, both sides are quenched, a more relevant analysis in the context of the present work is with one side adiabatic. This assertion will become clearer in the subsequent section. ), the temperature of the fluid at the bottom is higher than zero-g conditions and it does not fall below the top wall temperature. While the explanation for the higher temperature at the bottom wall (in case of top wall quenched) is nearly the same as discussed for both the walls quenched, of peculiar interest is the reasoning pertaining to the former case (bottom wall quenched).

One-dimensional analysis with one side adiabatic

For the schematic shown in Fig. 5-8(a), when the top wall is quenched, it induces a thin TBL near the top wall. Due to the action of self-weight and movement of the heavier fluid downwards, the fluid in the bulk is compressed/expanded near the bottom/top, respectively, thereby causing the temperature at the top to decrease more rapidly. However, since no TBL is formed at the bottom due to adiabatic condition, the temperature does not fall below the boundary value, as can be seen in Fig. 5-8(b) even at 𝑡 = 200𝑠. This is attributed to the absence of expansion of the bulk fluid from the bottom as in case of both walls quenched. Further, a higher temperature near the bottom wall can be explained by virtue of increase in temperature due to the compression under self-weight. In case when the bottom wall is quenched while the top wall is adiabatic (Fig. 5-7(a)), the density gradients near the bottom are formed owing to the temperature change in the TBL near the bottom. It is noticeable that in this case, the heavier fluid is already at the bottom and thus it may seem to be a stable configuration. However, it is observed the temperature drops below the wall temperature on the top adiabatic wall as can be seen in Fig. 5 case, it is logical to ascertain that same suite will be followed for higher dimensions.

See-saw motion of the thermal boundary layer

Rayleigh-vibrational instabilities were introduced in §4.4. In addition to the appearance of finger like structures, it is also observed that coupling between the temperature variations in the bulk due to vibration and temperature in the TBL leads to the see-saw motion of the TBL. This is illustrated in Fig. 5-13 wherein TBL seems to perform the said motion about an axis normal to the 2D plane. It is to be mentioned here that the difference in the temperature profiles at these 𝑦 positions is because finger like structures due to the Rayleigh-vibrational instability have protruded at 𝑦 = 0.3 𝑚𝑚 but are yet to reach 𝑦 = 1 𝑚𝑚. A closer look of the contour plots shows that the defined see-saw motion is the variation in the relative thickness of the TBL along the direction of vibration. This is illustrated more clearly in Fig. 5-14(a-b) where the thickness of the TBL is marked (approximatively) at two arbitrary 𝑥 positions by 𝛿 𝑇𝐵𝐿,1 and 𝛿 𝑇𝐵𝐿,2 and is found to change along the direction of vibration (thickness at 𝛿 𝑇𝐵𝐿,1 being greater than 𝛿 𝑇𝐵𝐿,2 when moving towards positive 𝑥 -direction and vice versa). Though a first glance at these findings may present their observance as naive and intuitive, it will be shown in subsequent section that it is not the change in the actual thickness of the TBL but the relative thickness which thereby presents a see-saw motion of the TBL (hence the use of relative thickness in caption of Fig. 567891011121314). Further, the see-saw motion of the TBL is defined separately from the oscillatory behavior of temperature field purely due to its clear visibility as compared to the oscillations in the bulk. In order to explain this behavior of TBL, initially the phenomenon of forced piston effect is described by virtue of which changes in the temperature field occur in the bulk and coupling of these changes with the TBL eventually leads to see-saw motion of the TBL.

Forced Piston Effect

In case of heat transport by piston effect, the temperature change in the bulk is caused by the conversion of mechanical energy of the pressure waves (propagating in the bulk due to high expansion or compression of SCF in the TBL) into thermal energy. On a similar basis, the change in the temperature field of the bulk/TBL due to pressure change caused by an external force is defined here as mechanical or forced piston effect (FPE). In the present work, this effect is caused by vibrational forces. Further justification for the use of this terminology can be attributed to the similar argument which holds for the piston effect, 𝑖. 𝑒. changes in the temperature field due to compression/expansion in a compressible fluid. The only difference is that while piston effect is caused by the propagation of acoustic or pressure wave, FPE arises due to local pressure variations created by external forces. In addition, the visibility and the importance of FPE is attributed to the high compressibility of SCFs. The effect of FPE will depend on the type of external forces acting on the system. For example, in case of angular rotation, the FPE will cause pressure to vary along the radial direction and thus temperature variations will be observed in the aforesaid direction. Thus, in the present case, the forced piston effect (FPE) refers to the piston effect by virtue of the external force, |𝑭| = 𝜌𝐴𝜔 2 𝑠𝑖𝑛(𝜔𝑡). In order to explain FPE in vibrations more clearly, the change in the temperature when the supercritical fluid is acted upon by vibrational forces only is analyzed. In such scenario, the top wall (𝐴𝐵) is maintained at 𝑇 𝑖 while others are adiabatic. The choice of these boundary conditions is motivated by the fact that it is intended to have changes in the temperature field only due to vibration and circumvent any change that may occur due to the main piston effect. the fact that the temperature increases at 𝑃2 while it shows opposite trend at 𝑃1. The increase in the maximum value of this temperature difference (denoted by ∆𝑇 𝑝𝑒𝑎𝑘 ) with amplitude can be explained as follows. A higher amplitude of vibration implies a higher acceleration and therefore a higher momentum associated with the fluid. When a fluid with higher momentum impacts the wall, it experiences a higher reaction force and thus a stronger compression. A similar reasoning can be asserted for the expansion side. A higher compression thereby leads to a higher temperature change (or difference 𝑇 𝑃1 -𝑇 𝑃2 ) as can be seen in Fig. 5-15(a). An analogous argument follows for higher frequencies, though results pertaining to the same have not been shown for the sake of brevity. The observations in Fig. 5-15(b) are attributed to the fact that on moving towards the center of the cell, the compressive forces from the wall experience a damping like action due to the high compressibility of the supercritical fluid.

The damping action mentioned here can be explained in a more physical sense as follows.

In case of an incompressible fluid, the fluid element can be considered as a solid, where the same magnitude of the force acting on one end will be transmitted to another end. However, in case of a compressible fluid, some part of the force will be utilized in compressing the fluid element (𝑖. 𝑒.

performing the compressible work on the fluid element) while the remaining fraction will be transmitted to the subsequent fluid element. A larger fraction of the force will be utilized in doing the work on the fluid element when the compressibility is higher. This reduction in force transmitted from one fluid element to another due to compressibility is referred to as damping-like action in the present case as it tends to dampen the action of applied force. Therefore, it implies that for a fluid with higher compressibility, a lower fraction of the force will be transmitted resulting in lesser change in pressure of the adjacent fluid element. As a result, change in temperature of the subsequent fluid elements will be inferior for a fluid with lower compressibility.

In order to quantify the above explanation, the ratio of maximum temperature change

Relative thickness of the TBL

Before moving further, the meaning of the relative thickness in the context of the present work is described as it forms the basis of understanding the see-saw motion. It is well known that the TBL is defined as a region where the temperature change is within 99% of the temperature difference between the bounding surface/wall temperature and the bulk. For example, consider a 1D case along the 𝑥 -direction where a fluid region is defined 0 ≤ 𝑥 ≤ 𝐿. Denoting the wall temperature (at 𝑥 = 0) by 𝑇 𝑊 and temperature of the bulk (𝑖. 𝑒. the region which is yet to experience the effect of thermal disturbance/temperature change at the boundary by virtue of thermal diffusion) by 𝑇 𝐵 , then the TBL represents a region (from 𝑥 = 0) for which the temperature 𝑇 obeys 𝑇 -𝑇 𝑊 ≤ 0.99(𝑇 𝐵 -𝑇 𝑊 ). It is thus clear that the definition of the TBL not only depends on the wall temperature (𝑇 𝑊 ) but also the bulk temperature (𝑇 𝐵 ). Therefore, for the same wall temperature but a different bulk temperature, it is intuitive to have different thickness of the TBL, say 𝛿 𝑇 𝐵1 and 𝛿 𝑇 𝐵2

where the subscript denotes different bulk temperatures (𝑇 𝐵1 and 𝑇 𝐵2 ) for two systems/regions 𝐵1 and 𝐵2, respectively. In the present work, the relative thickness is defined as the thickness of the TBL at a given location, say 𝐵2 𝑤. 𝑟. 𝑡 bulk temperature at some other location, say 𝐵1 and denoted 123 by 𝛿 𝑇 𝐵2 -𝑇 𝐵1 . Though, a more appropriate nomenclature will be 'relative bulk-temperature based thickness of the TBL', for the sake of brevity, it is referred to relative thickness. The need and justification for this description will become clearer in the following section wherein it will be shown that owing to the high compressibility and low thermal diffusivity of SCFs, noticeable and significant spatial variation in the bulk temperature may exist in the same closed cell due to the action of external forces. Thus, the use of relative thickness holds a special significance and forms the basis of see-saw motion thereby making the physical explanation not so obvious.

Mechanism of see-saw motion of the thermal boundary layer

The physical mechanism resulting in see-saw motion of the TBL is now explained and is further shown to be attributed to FPE in vibration. Consider a schematic of the process as shown in while the temperature profile along 𝑦-axis through these points will evolve as illustrated in Fig.

5-17(b).

Here, 𝑇 𝑤𝑎𝑙𝑙 denotes the wall temperature while the bulk temperature at that instant is represented by 𝑇 𝑏𝑢𝑙𝑘,𝑡 0 . The thickness of the TBL at all 𝑥 -positions will thus be defined by the temperature which satisfies 𝑇 𝑡𝑟𝑎𝑛𝑠𝑖𝑡𝑖𝑜𝑛 -𝑇 𝑤𝑎𝑙𝑙 = 0.99(𝑇 𝑏𝑢𝑙𝑘,𝑡 0 -𝑇 𝑤𝑎𝑙𝑙 ), where 𝑇 𝑡𝑟𝑎𝑛𝑠𝑖𝑡𝑖𝑜𝑛 denotes the temperature satisfying the above criteria. In order to analyze the evolution of the temperature in the presence of vibrational acceleration, a case when the cell moves in positive 𝑥-direction is analyzed. By virtue of FPE, the pressure along the line passing through 𝐶2 in 𝑦-direction will increase due to compression thereby resulting in a uniform increase in the temperature of fluid elements therein. The temperature will thus evolve as shown by dashed red line in Fig. 5-17(c). On the other side near the wall (𝐴𝐷), the temperature will decrease (due to expansion) leading to the temperature field along the 𝑦 axis as represented schematically in Fig. 5-17(c) by dashed blue lines. It can thus be observed that at any 𝑥, FPE in vibration causes a change in temperature in the bulk and the TBL. Thus, the TBL (say corresponding to 𝑥 -position at 𝐶1) will be represented by temperature along a particular location. In such scenario, the average bulk temperature would nearly be same as at the center owing to the symmetric effects (about the center) of FPE and thus it is more reasonable to ascertain it as the reference state. Therefore, relative to this reference state, the relative bulk temperature based TBL thickness or relative thickness of the TBL along the 𝑥 -direction is maximum (minimum) near wall 𝐴𝐷(𝐵𝐶) for positive direction of vibration as shown in Fig. 5 (b,d) that the transition from the TBL to bulk occurs at nearly same 𝑦 at various 𝑥 positions and thus the actual thickness of the TBL at respective 𝑥 remains unaffected. However, analyzing the system as a whole as depicted in Fig. 5-18(a,c), where the bulk average temperature would be nearly the same as the one in the middle (as explained above), it can be seen that the relative thickness (as defined previously) changes along the direction of vibration. It is to be mentioned that low value of amplitude has been purposely chosen so as to evade any effects of Rayleigh-vibrational instability.

As shown in Fig. 5-18(a), when the cell moves in the positive 𝑥 -direction, the relative thickness is higher at 𝐿𝑀 (𝛿 𝑇 𝐿𝑀 -𝑇 𝑀 ) while it is lower at 𝑅𝑀 (𝛿 𝑇 𝑅𝑀 -𝑇 𝑀 ) defined 𝑤. 𝑟. 𝑡 to temperature at the middle (or bulk average). The scenario reverses with change in the direction of vibration as shown in Fig. 5-18(c) and thus supports the above explanation.

It is also worth observing that the relative thickness decreases on moving from expansion to compression side (Fig. 5-18(a,c)). It was mentioned in §5.2.1 that the effect of FPE decreases on moving towards the center of the cell. As a result, the difference between the local bulk temperature and bulk temperature at the center and thus the relative thickness, which implicitly indicates this difference, decreases on moving towards the center as observed in Fig. 5-18(a,c) and represented schematically in Fig. 5-17(a) by 𝑀𝑁 𝑥+ and 𝑀𝑁 𝑥-. Thus, more fundamentally, the see-saw-motion can also be defined as change in the relative thickness of the TBL along the direction of vibration due to the temperature changes by virtue of FPE in vibration.

Factors affecting FPE in vibration and see-saw motion

The see-saw motion described above implicitly represents the temperature difference caused by virtue of FPE resulting in relative bulk temperature. Thus, analyzing the effect of various factors on temperature change caused by FPE in vibration can provide a first-hand means to characterize the see-saw motion. The effect of various parameters on ∆𝑇 𝑚𝑎𝑥 using a simple 1D analysis is thus presented here. An important remark not be missed here is that the see-saw motion will become more pronounced with advancement in time. This is ascribed to the fact that with growth in the thickness of the TBL at the center, the difference in the relative thickness on extreme ends will increase making this worth noticeable.

• Effect of proximity to the critical point and acceleration

• Effect of cell size

The observance of temperature changes by FPE in vibration (and thus see-saw motion) is physically attributed to the reaction force acting on the SCF when it impacts the wall. It was shown in section §5.2.1 that a higher momentum by virtue of higher acceleration resulted in higher ∆𝑇 𝑚𝑎𝑥 .

Another factor which may affect the momentum is the cell volume (and thus mass of SCFs). This aspect is also highlighted in Fig. 5-19 (a) wherein ∆𝑇 𝑚𝑎𝑥 is plotted for 2 different cell sizes, ℎ = 3.5 𝑚𝑚 and 7 𝑚𝑚. A higher value observed for a larger cell is ascribed to a higher mass and thus momentum which implies that for the same vibrational acceleration, the fluid impacts the wall with a higher force. Consequently, a higher reaction force causes greater change in pressure and thus temperature.

Chapter summary

In this chapter, two captivating phenomena when SCFs are subjected to simultaneous quench and vibration are explained which are observed along-with Rayleigh-vibrational instabilities ( §4.4).

The first one describes the drop in the temperature of the SCF below the imposed boundary temperature. These regions are defined as sink-zones. In order to explain these observations, several 1D cases with different boundary conditions are analyzed in the presence of constant acceleration and it is found that the temperature gradient must be parallel to the direction of acceleration for the formation of sink-zones. Unlike the 1D case where the direction of acceleration is along the temperature gradient, the appearance of sink-zones in case of 2D is attributed to oblique flows with direction of acceleration near the corners. The component of vibration acceleration along these flow lines leads to similar conditions as in 1D flow and thus the appearance of sink-zones.

The second phenomena pertain to see-saw motion of the TBL. It is shown that pressure changes due to external forces (vibrational acceleration in the present case) can cause significant changes in the temperature field due to diverging behavior of the thermo-physical properties in near-critical fluids and is termed as forced piston effect (FPE). The effect of FPE diminishes on moving from the walls towards the bulk which is ascribed to damping like action by virtue of high compressibility of SCFs. It is further shown that even though the actual TBL thickness may be same along the direction of vibration, the relative thickness (defined 𝑤. 𝑟. 𝑡 relative local bulk temperature) will change along the direction of vibration which over a period of vibration forms a see-saw motion of the TBL.

In the next chapter, the elementary test cases are presented using phase-field modelling for twophase fluid system near the critical point. An important aspect pertaining to this two-phase system is that, both the phases, liquid as well as gas, are highly compressible.

field parameter which can lead to diffusion. Even though it can be reduced by decreasing the mobility (𝑀 0 ) or refining the mesh, the complete removal remains a big challenge. In some cases, these can become as large as the characteristic velocities of the problem thereby altering the real physical behavior. the currents is higher. This is primarily attributed to a larger interface region and thus more diffusion to attain equilibrium.

Another aspect that affects spurious currents as described earlier is a higher surface tension and density ratio. In order to investigate this, the system is analyzed 1𝐾 below the critical point where surface tension and density difference between both the phases is relatively higher (see Table 6-1). Fig. 6-3 shows the evolution of maximum velocity for 𝑎 = 1 𝑚𝑚 and 𝐶𝑛 = 0.02 . It can be observed that at a higher value of mobility (𝑀 0 = 10 -7 S.I), the system initially attains a lower value of spurious currents only to overcome initial unsteadiness. Once that has reached, a higher mobility causes diffusion thereby resulting in increasing spurious currents. With the reduction in mobility, there is a slight increase in the spurious currents which however shows a diminishing behavior with time. It is to be mentioned here that the velocity plotted herein is the maximum velocity in the whole domain which will be prevalent only near the interface region while it is nearly zero in the rest of the domain. This is illustrated in Fig. 6-4 using contour plots of 𝑥 and 𝑦 velocity at different time instances for 𝑎 = 1𝑚𝑚, 𝐶𝑛 = 0.01, 𝑀 0 = 10 -13 S.I when 100𝑚𝐾 below the critical point.

Transformation of elliptical bubble to circular shape

The interfacial energy depends on the interface area which in case of 2D problem is given by the perimeter of the interface. In order to minimize the energy, the interface will finally try to attain a shape which minimizes its energy. Since for a constant area, the perimeter of a circle is lower, an initial elliptical vapor bubble will transform into a circular shape. This aspect is investigated for the schematic as shown in the Fig. 6-5, for 𝑎/𝑏 = 2 with 𝑎 = 2𝑚𝑚 ( 100𝑚𝐾 below the critical point).

Fig. 6-6 shows time evolution of the vapor bubble from elliptical to circular shape. One of the primary test in this regard is that once the bubble becomes stable, it should satisfy the Young-Laplace relation for the evaluation of surface tension given by, In order to verify the coalescence, the final radius of the bubble after coalescence is compared with the analytical value. In order to obtain the diameter (and thus radius) from simulation, the mean distance between 𝜙 = 0.5 along 𝑥 = 5𝑚𝑚 and 𝑦 = 5𝑚𝑚 was evaluated. The numerical simulation gave a value of 1.142 𝑚𝑚 as compared to analytical value of 1.131 𝑚𝑚 , an error of 1% which is well within the acceptable limits.

Separation of liquid and vapor phases

As a final study, which highlights one of major strength of phase-field modelling to model the appearance of interface, the spinodal decomposition of a liquid-vapor mixture system is investigated in weightlessness conditions. The initial conditions correspond to a square cavity with 𝜙 0 = 0.623 and 𝜙 𝑟𝑎𝑛𝑑 is a random distribution between -0.5 × 10 -2 and 0.5 × 10 -2 . This is similar to the study of Lamorgese and Mauri [START_REF] Lamorgese | Diffuse-interface modeling of liquid-vapor phase separation in a van der Waals fluid[END_REF], wherein Gaussian white noise was added for random perturbation. As discussed in §3.4.2, the mixture will separate into individual phases as it is an energetically favored condition. An important difference from the classical spinodal decomposition (as described in § 3.4.2) is that here, the momentum equation is also solved along with thereby including the effect of surface tension. It is thus possible to simulate the appearance of interface using the described phase-field model. The model can be further developed to study the continuous transition from single-phase (supercritical) to two-phase region (sub-critical) region. The challenges and possible solution methodologies is described in perspectives, §7.2.

Chapter Summary

The phase-field model based on mass fraction as phase-field parameter was investigated for various cases in sub-critical regime near the critical point. These primarily include spurious current in quiescent stagnant vapor bubble, transformation of an elliptical bubble to circular shape and coalescence of two droplets favored by minimizing of the energy. In addition, the spinodal decomposition of a liquid-vapor system at equilibrium with each other is also analyzed wherein appearance of interface is observed. This is an essential feature which is required to extend the model for continuous transition from supercritical (single-phase) to sub-critical regime (two-phase) wherein simulating the appearance of the interface forms a vital element.

CHAPTER 7. CONCLUSIONS AND PERSPECTIVES

Conclusions

In the present thesis, motivated by the long-term goal to model a continuous transition from supercritical to subcritical state in the realm of continuum mechanics, a compressible model is developed to analyze the flow behavior in supercritical fluids. Several intriguing phenomena, attributed to the behavior of thermo-physical properties near the critical point, were observed and physical explanations presented therein. A preliminary analysis is also presented for sub-critical state using phase-field methodology in isothermal conditions. The conclusions of the whole work can be summarized as follows,

• A compressible model, which calculates density directly from the mass conservation and without the need to use any pressure velocity coupling algorithm is described. This has been made possible by directly introducing the dependence of pressure on density and temperature into the momentum equation which precludes any unknown pressure present in the momentum equation. The model is subsequently verified by analyzing the propagation of acoustic wave in water and the results are found to match well with the analytical results.

In addition, the model is also validated by comparing numerical results with experimental data and a close match is observed. This highlights the capability of the model to analyze fluid system with very low compressibility (such as water) to highly compressible fluids such as supercritical fluids.

• The numerical investigations for the propagation of acoustic wave in water reveal interesting numerical phenomenon wherein the accuracy of the solution is observed to be affected by Courant number in unusual way, oscillatory solution for low values of Courant number while accurate numerical solutions at higher values. This has been explained by virtue of contribution to the error arising from phase-lag and difference between group velocity of numerical waves and physical waves, in addition to the conventional error growth rate. While the classical analysis renders implicit schemes as unconditionally stable for linear wave equation and stability decreases with increasing Courant number for nonlinear systems, the correct and complete error dynamics brings out the true behavior. It is thus shown that even though the scheme may be neutrally stable, the error growth rate being unity, it can still yield erroneous results due to phase-lag and dispersion errors.

• The mathematical model has been used to analyze the behavior of supercritical fluids when simultaneously subjected to thermal quench and mechanical vibration. This leads to primarily two types of instabilities, Rayleigh-vibrational instabilities and parametric instabilities, when the direction of temperature gradient and vibration are respectively normal and parallel to each other. The physical mechanism leading to these instabilities has been described and the effect of various parameters, such as quench rate, frequency, proximity to the critical point on the critical amplitude leading to onset of these instabilities is investigated. Further, the wavelength in Rayleigh-vibrational instabilities is found to decrease over a period of time and can increase at higher acceleration.

• An interesting phenomenon observed along-with Rayleigh-vibrational instabilities is the drop of fluid temperature below the boundary limit in supercritical. These regions have been termed as sink-zones. The problem is first explained in 1D configuration wherein the direction of acceleration is along the direction of temperature gradient. It is observed that by virtue of high compressibility and coefficient of thermal expansion, the temperature of the SCF increases by virtue of self-weight. The explanation is subsequently extended to 2D configuration, wherein it is observed that near the corners of the cell, the conditions are analogous to 1D configuration during one half time-period of the vibration.

• Another notable observation alongside Rayleigh-vibrational instabilities pertains to the seesaw motion of the TBL. It is shown that pressure changes due to external force (vibrational acceleration in the present case) can cause significant changes in the temperature field. We describe this phenomenon as forced piston effect (FPE) and its effect is observed to diminish on moving from the walls towards the bulk. It is further shown that even though the actual TBL thickness may be the same along the direction of vibration, the relative thickness (defined 𝑤. 𝑟. 𝑡 relative local bulk temperature) will change along the direction of vibration which, over a period of vibration, forms a see-saw motion of the TBL.

• A phase-field model is presented to analyze the flow behavior in sub-critical state for isothermal systems. The relevance and advantage of using phase-field method in the context of near-critical fluids is ascribed to the diverging behavior of the interface thickness on approaching the critical point which is in coherence with the assumption of finite thickness of the interface region in phase-field modelling. Another notable feature is that it permits the appearance of interface from a homogeneous mixture which is presented though the case of phase-separation. The model is further used to analyze various elementary test cases, such as the stability of a stagnant bubble, coalescence of two bubbles etc.

Perspectives

While the current work contributes to the study of the behavior of near critical fluids using mathematical and numerical method, it opens up new horizons of investigations which can be considered to be of significant importance, both from the fundamental as well application perspectives. Some of these can be described as,

Instabilities in supercritical fluid

• The influence of thermo-mechanical coupling has been shown to lead to interesting results.

However, in these studies, the external force is spatially uniform throughout the fluid cell by virtue of longitudinal accelerations. This aspect will be significantly different if the cell is subjected to angular rotation wherein due to centrifugal force a spatial varying force will exist in the system. This can lead to interesting flow features and needs to be investigated.

The relevance of this perspective can be considered from the physical fact that, even a small longitudinal vibration when acting eccentrically 𝑤. 𝑟. 𝑡 to a fixed point inherits rotational characteristics and can lead to additional changes in the flow behavior.

• In several applications, the fluid is not confined in cell but flows through channels, such as in the nozzle of rocket engines. While the current study was primary confined to closed systems, the effect of simultaneous vibration and thermal perturbation in open or flow system can be considered as the next step of investigation of instabilities in supercritical fluids.

Phase-field modelling in sub-critical regime with temperature change

The current work analyzed the behavior near the critical point in two-phase using phase-field modelling in isothermal conditions. With promising elementary results, especially pertaining to the separation of a mixture into the liquid and vapor phase, the current model can be further developed to model a continuous transition from supercritical to sub-critical regime and vice-versa. A brief discussion on the challenges and possible roadmap is presented now for future investigations:

• The current model considers the flow behavior in isothermal conditions, 𝑖. 𝑒. without any temperature or phase change. The first important perspective is to extend the model to nonisothermal systems. This is an important aspect especially attributed to highly compressible liquid and vapor phase which can thereby induce strong thermo-mechanical coupling.

• One of the primary investigation for non-isothermal system will be to model phase-change in the sub-critical regime, for example a test case of a saturated liquid drop suspended in superheated vapor which should subsequently be converted to vapor phase. However, an important point to be considered here is that, the mass-fraction of each phase will not be conserved. This calls for an alternative phase-field model, such as Allen-Cahn model.

• An important role played in the phase-field model is an appropriate choice of free energy function. Hence, it is essential to develop an appropriate function which can include all the dynamics of phase transition, such as it should represent minima corresponding to two phases at all temperatures below the critical point while above the critical point, it should have only single minima thereby favoring existence of a single phase. Further, owing to the possible phase-change within the sub-critical regime, the function should be able to include the effects of latent heat. It is to be noted that for non-isothermal conditions, it is necessary to have function as dependent on temperature conditions.

• It is worth noticing that parameters (𝛼) and (𝛽) are related to the interface thickness and surface tension. While surface tension tends to vanish on approaching the critical point, interface thickness tends to diverge. This additional constraint which arises from physical attributes needs to be considered when deriving the function. A simple means to achieve this could be to consider dependence 𝛽 on (𝑇 -𝑇 𝑐 ) with a higher power as compared to 𝛼.

It is worth mentioning here that while the experimental studies of thermo-vibrational studies in subcritical regime near the critical point are available in literature [START_REF] Gandikota | Frozen-wave instability in near-critical hydrogen subjected to horizontal vibration under various gravity fields[END_REF][START_REF] Gandikota | Faraday instability in a near-critical fluid under weightlessness[END_REF][START_REF] Gandikota | Dynamic equilibrium under vibrations of H2 liquid-vapor interface at various gravity levels[END_REF], the numerical investigations are still an open topic. Some future perspective in this regard can be described as,

• It is known that when approaching the critical point, the surface tension vanishes. Thus, it will be interesting to analyze thermo-capillary phenomena in near-critical fluids.

• An important aspect which governs the flow evolution in phase-field modelling is the choice of appropriate time step because of various time scales involved in the analysis, such as diffusion time scale due to the gradient of phase-field parameter, time scale of external force, viscous time scales. Further, when the system will be analyzed in non-isothermal conditions near the critical point, thermal diffusion and piston effect time scale will also play an important role. In order to effectively capture the exact physical mechanism, it is essential to have appropriate

C. Error propagation equation

In where k is the wavenumber with initial condition given by, 𝑢 𝑗 0 = ∫ 𝑈 0 (𝑘)𝑒 𝑖𝑘𝑥 𝑗 𝑑𝑘 [START_REF] Sengupta | Error dynamics: Beyond von Neumann analysis[END_REF][START_REF] Sengupta | High Accuracy Compact Schemes and Gibbs' Phenomenon[END_REF].

Defining the error growth rate as ) , (

) ) 1 ( , ( ) , ( t n x k U t n x k U t x k G j        
, the solution at 𝑛 𝑡ℎ time step as described by Eq. (C.2) can also be represented by [START_REF] Sengupta | Error dynamics: Beyond von Neumann analysis[END_REF][START_REF] Sengupta | High Accuracy Compact Schemes and Gibbs' Phenomenon[END_REF], It can be observed that, unlike conventional error analysis [START_REF] Charney | Numerical Integration of the Barotropic Vorticity Equation[END_REF], wherein the contribution to error arises only from error growth rate, the correct error dynamics has contribution from phase lag as well as the difference in energy propagation (or group velocity) of the numerical as physical wave..

𝑢

D. Numerical characteristics for second order central difference and implicit first order forward Euler in time numerical scheme

In order to draw out the characteristics of the numerical scheme used in the current work (as mentioned above), the finite difference discretization methodology is used here for the sake simplicity. The 1D wave equation (C.1) as described above, in discretized form using the aforementioned numerical scheme can be thus be written as,

E. Free energy of a non-homogenous system

The free energy functional for an isotropic system with non-homogeneous composition or density, as derived by Cahn and Hilliard [START_REF] Cahn | Free Energy of a Nonuniform System. I. Interfacial Free Energy[END_REF] is presented here. Even though their analysis was presented with composition as a parameter, it was reported that any intensive scalar property can be used as a parameter. Thus, let us consider a general parameter 𝜙, which can represent any intensive scalar property, such as composition or density etc., which varies spatially in the system and can be used to distinguish between two different phases or fluids.

Let 𝐹 𝑇 denote the total free energy of the system (here the free energy can be either Gibbs or Helmholtz depending on the case). As the free energy depends on the parameter 𝜙, we can thus write 𝐹 𝑇 in terms of 𝜙 and its gradients as, or 𝐹 𝑉 = 𝑓(𝜙) + 𝜅(𝜵𝜙) 2 . Here 𝑓 0 (𝜙) is contribution due to the bulk (𝑖. 𝑒. the bulk free energy) and 𝜅(𝜵𝜙) 2 is contribution due to the interface (interface energy). However, since we are interested mainly in interface dynamics, the free energy contribution from the bulk can be subtracted from the overall free energy and thus 𝑓 0 (𝜙) will then represent free energy 𝑤. 𝑟. 𝑡 base state which is the bulk energy.

𝐹

Note:

In general, the bulk free energy function, 𝑓 0 (𝜙), is in the form of a double well function of variable 𝜙 and thus dimensionless. A more generic expression of free energy is represented by 𝐹 𝑉 = 𝛽 ̅ 𝑓 0 (𝜙) + 𝛼 ̅ 2 |𝜵𝜙| 2 (E.9)

where the coefficients 𝛽 ̅ and 𝛼 ̅ are related to physical parameters, interface thickness (𝜉) and surface tension (𝜎) and have the necessary dimensions. Further, in the current work in §2.2.8, the free energy function, 𝐹, is defined per unit mass and thus we have 𝐹 𝑉 = 𝜌𝐹 while the coefficients are 𝛼 ̅ = 𝛼𝜌 , 𝛽 ̅ = 𝜌𝛽.

E.1 Interface Profile in 1D case

The free energy function as defined in previous section by Eq. (E.9) describes the dependence of free energy on parameter 𝜙 where 𝜙 = 𝜙(𝑥) and includes the effect of interface energy. The interface profile will thus depend on 𝜙 such that it minimizes the free energy, 𝐹 𝑉 . In the theory proposed by Cahn and Hilliard [START_REF] Cahn | Free Energy of a Nonuniform System. I. Interfacial Free Energy[END_REF], which was primarily described for composition as the parameter, the variational derivative of free energy was linked with chemical potential. On similar terms, the variational derivative of free energy with respect to 𝜙 can be defined as potential,

F. Coupling between phase-field (Cahn-Hilliard) and Navier Stokes

In order to couple the phase-field parameter (𝜙) with the Navier-Stokes equation, appropriate changes need to be incorporated in the governing flow equations so as to include the effect of interface dynamics. Among several means described in literature, the method adopted by Jacqmin is presented here [START_REF] Jacqmin | An energy approach to the continuum surface tension method[END_REF][START_REF] Jacqmin | Calculation of Two-Phase Navier-Stokes Flows Using Phase-Field Modeling[END_REF]. The main idea in this approach was based on the following two effects,

• the movement of interface can cause it to become thick/thin thereby causing a change in free energy.

• this effect of increase in energy is countered by a diffuse-interface force resulting in change of kinetic energy.

The change in energy due to convection will be attributed to change in parameter 𝜙 and thus can be represented by, Now, as suggested by Jacqmin [START_REF] Jacqmin | An energy approach to the continuum surface tension method[END_REF][START_REF] Jacqmin | Calculation of Two-Phase Navier-Stokes Flows Using Phase-Field Modeling[END_REF], the above Eq. (F.9) can be reduced as,

𝐸
𝑓

iiiTitre:

  Modélisation numérique des fluides fortement compressibles proches du point critique Résumé court Un fluide porté à une température et pression supérieures à celles du point critique est communément appelé fluide supercritique. Ce fluide possède des propriétés particulièrement intéressantes à cheval entre celles des gaz et celle des liquides. En effet, la masse volumique d'un fluide supercritique est proche de celle d'un liquide tandis que sa viscosité est proche de celle d'un gaz. Une des caractéristiques particulières de ces fluides quand ils s'approchent du point critique est que plusieurs des propriétés thermo-physiques montrent un comportement singulier (compressibilité divergente, diffusivité thermique évanescente etc). Dans ce travail, un modèle mathématique basé sur les équations de Navier-Stokes couplées à celle de l'énergie est proposé afin d'étudier les écoulements de ces fluides très proches de leur point critique. La validation du modèle a été effectuée sur un problème de propagation d'onde acoustique dans l'eau. Nous avons ainsi observé que des solutions précises avec des schémas implicites pour des systèmes non linéaires sont possibles avec des nombres de Courant élevés. L'étude des écoulements dans des fluides supercritiques, lorsqu'ils sont assujettis à une trempe thermique et à une vibration simultanées ont montré que de telles conditions pouvaient conduire à la formation d'instabilités thermo-vibrationnelles, en particulier les instabilités de Rayleigh-vibrationnelles et paramétriques.

Fig. 1 -

 1 Fig. 1-1: (a) 𝑃 -𝑉 -𝑇 surface plot for a pure substance which contracts on freezing (b) 𝑃 -𝑇 plot illustrating critical point and supercritical fluid region.

Fig. 1 - 2 :

 12 Fig. 1-2: Illustration of the Widom line, an extension to the coexistence line at supercritical pressure and as a marker of the crossover between supercritical liquid-like and gas-like states. It is a marker of the crossover between supercritical liquid-like and gas-like states. The dashed line at Z=PV/(RT)=0.95 denotes the transition to an ideal gas. The reduced density is 𝜌 𝑟 = 𝜌 𝜌 𝑐 ⁄ , 𝐶 𝑃,𝑟 = 𝐶 𝑃 𝐶 𝑃,𝐼𝐺 ⁄ , 𝐶 𝑃,𝐼𝐺 = γR/(γ -1), where

  development of mathematical and numerical models to analyze the flow behavior in near-critical fluids, primarily in supercritical state. The thesis is organized as follows, Chapter 2 describes the mathematical model developed to analyze highly compressible supercritical fluid. The current model directly incorporates the dependence of pressure on density and temperature into the momentum equation thereby circumventing the need of any pressure velocity coupling algorithm. In addition, the density is calculated directly from mass-conservation without the need of any explicit equation of state for the calculation of density. In order to analyze the twophase flow in sub-critical regime near the critical point, phase-field model is developed to analyze flow in sub-critical state under isothermal conditions. One of the primary advantage of phase-field modelling is ascribed to its capability to model the appearance of interface. Chapter 3 describes the numerical algorithm to solve the mathematical model developed in Chapter 2 for a single and two-phase flow which is succeeded by validation studies of these models. In addition, an in-depth error analysis is presented which highlights unusual behavior of Courant number on the accuracy of solution.

•)

  various thermo-physical properties, namely density, isothermal compressibility 𝜒 𝑇 = , 𝑘. A general mathematical model solving conservation laws[START_REF] Henk | An Introduction to Computational Fluid Dynamics: The Finite Volume Method[END_REF] can be written in non-conservative form (𝑑/𝑑𝑡 representing the material derivative) as, Conservation of mass (continuity equation) + 𝜵. (𝜇𝜵𝑽) + 𝜵((𝜆 + 𝜇) 𝜵. 𝑽) + 𝑭 (2.1.2) Here, 𝜇 and 𝜆 correspond to the classical shear viscosity and compression (bulk) viscosity or second coefficient of viscosity, respectively. Further F and 𝑃 denote volumetric force and the thermodynamic pressure, respectively. • Conservation of internal energy 𝜌 𝑑𝑒 𝑑𝑡 = -𝑃𝜵. 𝑽 + 𝑞̇+ 𝜵. (𝑘𝜵𝑇) + 𝜑 (2.1.3) or in terms of temperature, 𝑇 as, + 𝑞̇+ 𝜵. (𝑘𝜵𝑇) + 𝜑 (2.1.4)Here 𝑒 is the internal energy while φ = 𝜆(𝛁 • 𝑽)2 being the tensor of deformation rate (The Einstein summation convention on repeated indices is applied) while 𝑞̇ denotes the volumetric heat source.In addition, for a compressible fluid, the dependence of pressure on density and temperature, is usually governed by the equation of state, 𝑃 = 𝑃(𝜌, 𝑇) (2.1.5)Now, infinitesimal variations of temperature, pressure, and the density with time can be written as,

Fig. 2 - 1 :

 21 Fig. 2-1:Schematic illustration of sharp and smooth transition (interface).

Fig. 2 - 2 :

 22 Fig. 2-2: Schematic of double-well potential (bulk-free energy) considered in the present work.

2 𝑀.

 2 phase field model for the conserved parameter comprises of solving the Cahn-Hilliard form of equation. Here, we derive the same form of equation considering mass fraction as the parameter as described by Guo and Lin [55]. To begin with let us first define some terms required for the derivation. Consider a small control volume 𝑉 and the two fluids be represented by 𝑖 = 1,2 individually filling up a volume 𝑉 𝑖 . Let 𝐶 𝑖 denote the volume fraction of each fluid 𝑖. 𝑒. 𝐶 𝑖 = 𝑉 𝑖 𝑉 . Further, these fluids can mix with each other in the interfacial region (as per the theory of Gibbs [29]). Let 𝑀 denote the total mass of the mixture with 𝑀 𝑖 being the mass of 𝑖 th fluid. The mass fraction of each phase be denoted by 𝛾 𝑖 , 𝑖. 𝑒., 𝛾 1 The density of each fluid in each phase can be represented by 𝜌 𝑖 = 𝑀 𝑖 𝑉 𝑖 . In addition, we also define, an apparent or partial density (local volume-averaged mass density) taken over a 𝑉 for each fluid as 𝜌 ̃𝑖 = 𝑀 𝑖 𝑉

  solution of the mathematical model developed in the previous chapter is not feasible by analytical means owing to the coupled and non-linear terms involved therein. The general solution methodology in such situations pertains to discretizing the set of governing equations on discrete points in space and time and converting them into a set of linear algebraic equations. The discretized equations thereafter constrained by appropriate boundary conditions give solutions which represent the physical behavior of the system. This chapter presents the numerical methodology used to solve the mathematical model developed in the CHAPTER 2. A brief introduction to the Finite volume method is followed by the description of the home-made code Thetis, developed at the laboratory I2M, which is used in the current work. The developed models (mathematical and numerical) for the single-phase flow are validated by analyzing the propagation of acoustic wave in water. For the phase-field model, the Cahn-Hilliard equation is validated by analyzing the classical problem of spinodal decomposition.

Fig. 3 - 1 :

 31 Fig. 3-1: Schematic illustration of staggered grid in Finite Volume Method.

Fig. 3 - 2 :

 32 Fig. 3-2: Schematic arrangement of grid in Thetis.

  3.1, the Cahn-Hilliard equation.(2.2.29) is solved after the momentum equation. Thus, the additional source term to couple the phase-field equation and Navier-Stokes equations (as given by (2.2.30)) is discretized as, 𝐹 𝑆𝑇 (𝑛) = -𝛼 𝜵. (𝜌 (𝑛) 𝜵𝜙 (𝑛) ⊗ 𝜵𝜙 (𝑛) ) + 𝜵(𝛼𝜌 (𝑛) |𝜵𝜙 (𝑛) | 2 ) (3.3.7)

Fig. 3 - 3 :

 33 Fig. 3-3: 1D Schema of the test case for propagation of acoustic wave in water.

Fig. 3 - 4 :

 34 Fig. 3-4: 𝑃/𝑃 𝑎𝑡𝑚 at the inlet for 𝑁 𝑐 = 0.11,1 and 10 with corresponding grid sizes.

Fig. 3 -

 3 Fig. 3-4 shows the results ( 𝑃 𝑃 𝑎𝑡𝑚 ) obtained for different Courant numbers (𝑁 𝑐 = 𝑐 𝑠 ∆𝑡 ∆𝑥 where ∆𝑡 is the timestep and ∆𝑥 is the grid size ) given by the values, 0.11, 1 and 10, with corresponding grid sizes as marked in the figure ) at point 𝑀 marked in Fig. 3-3. In terms of time steps, these cases

  Error growth rate, 𝐺 𝑗 (𝑘∆𝑥, ∆𝑡) = 𝑈(𝑘∆𝑥,(𝑛+1)∆𝑡) 𝑈(𝑘∆𝑥,𝑛∆𝑡) at any node 𝑗, with |𝐺 𝑗 | = 1 representing neutrally stable scheme which is absolutely desirable for the removal of dissipation error (the last term on the right-hand side (RHS) of Eq. (3.4.3)); • Phase error proportional to (1 -[ 𝑐 𝑁 𝑐 ]) which quantifies the lagging or leading of numerical wave (𝑐 𝑁 is defined as the phase speed) with respect to the physical wave speed (𝑐)) (the first term on the RHS of Eq. (3.4.3)); • Ratio of group velocity (𝑉 𝑔𝑁 ) and physical wave speed (𝑐), 𝑉 𝑔𝑁 𝑐 , representing the energy propagation speed and hence measuring the dispersion error, as given by the second term on the RHS of Eq. (3.4.3).The significance of the three aforementioned parameters in evaluating the behavior of any numerical scheme comes from the expression of the error as given by Eq. (3.4.3). The first and last terms in the RHS of Eq. (3.4.3) explicitly represent effects of phase lag and error growth rate, respectively, whereas the second term implicitly represents, in energy propagation. Note that the error equation depends upon 𝐿𝑛|𝐺| as one of the contributor (third term on RHS of Eq. (3.4.3)). For all implicit methods, |𝐺| ≤ 1, and this term vanishes for specific combinations of numerical parameters when |𝐺| = 1. However, for such parameter combinations of the implicit method, numerical solutions are not always trouble-free, as the errors can be created due to the contribution of the first two terms in Eq.(3.4.3). In order to numerically capture the propagation of physical wave in an efficient way, it is desired to have a strictly neutrally stable scheme with vanishing phase error (1 -𝑐 𝑁 𝑐

Fig. 3 -

 3 Fig. 3-5 illustrates the contour plots of the characteristics of the aforementioned numerical scheme: modulus of error growth rate |𝐺 𝑗 |, energy propagation speed ( 𝑉 𝑔𝑁 𝐶 ) and phase error metric

  ) illustrates the operating points marked for different grid sizes for Courant number of 0.11 on the contour plots of |𝐺 𝑗 |, (

Fig. 3 - 6 :

 36 Fig. 3-6: Contour plots for 𝑁 𝑐 = 0.11 in (𝑁 𝑐 , 𝑘Δ𝑥)-plane (a) |𝐺 𝑗 | (b) 𝑉 𝑔𝑁 𝑐 (c) (1 -𝑐 𝑁 𝑐

4 . 1 . 1 ,

 411 𝑡 = 0.67 𝑚𝑠 corresponds to the analytical time at which sharppeak is observed. Thus, at all previous times, an accurate and precise numerical solution should lead to 𝑃 𝑃 𝑎𝑡𝑚 ≈ 15. However, as can be observed from Fig. 3-4 (for ∆𝑡 = 10 𝜇𝑠), 𝑃 𝑃 𝑎𝑡𝑚 deviates from the desired value for all Courant numbers considered in the study. In order to explain this behavior more clearly, we draw out the difference between the two extreme cases, smallest Courant number 55 with least spatial resolution (𝑁 𝑐 = 0.11 and 𝛥𝑥~125 𝑚𝑚) and largest Courant number with highest spatial resolution (𝑁 𝑐 = 10 and 𝛥𝑥~0.0149 𝑚𝑚).

Fig. 3 - 7 :

 37 Fig. 3-7: Contour plots for 𝑁 𝑐 = 10 in (𝑁 𝑐 , 𝑘Δ𝑥)-plane (a) |𝐺 𝑗 | (b) 𝑉 𝑔𝑁 𝑐 (c) (1 -𝑐 𝑁 𝑐 ). The marked points corresponding to fundamental wavenumber 𝑘 2𝜋 ~ 1 𝑚 -1 .

Fig. 3 -

 3 Fig. 3-8 (a-c): Contour plots of (a) 𝐺 𝑗 , (b) 𝑉 𝑔𝑁 𝑐 , (c) (1 -𝑐 𝑁 𝑐 ) in (𝑁 𝑐 , 𝑘Δ𝑥)-plane for 𝑁 𝑐 = 0.11 and ∆𝑥 ~125 𝑚𝑚. The marked point corresponds to second dominant wavenumber ( 𝑘 2𝜋 = 4 𝑚 -1 ) for solution at time 𝑡 = 0.6 𝑚𝑠.

Fig. 3 -

 3 Fig. 3-9(a-c): Contour plots of (a) |𝐺 𝑗 |, (b) 𝑉 𝑔𝑁 𝑐 , (c) (1 -𝑐 𝑁 𝑐 ) in (𝑁 𝑐 -∆𝑥) plane for 𝑁 𝑐 = 10 and ∆𝑥 ~0.0149 𝑚𝑚 for different wavenumbers for solution at time 𝑡 = 0.6 𝑚𝑠.

Fig. 3 -

 3 Fig. 3-10 showing the propagation of the pressure wave along the depth at various time intervals for 𝑁 𝑐 = 1 and 𝛿𝑡 = 0.1 𝜇𝑠. The distance travelled by the acoustic wave front after a period of 𝑡 1 = 0.2𝑚𝑠 is marked by point 𝑀 1 in Fig. 3-10 and matches well with the analytical values calculated using the speed of sound in water as mentioned in Table 3-1 (distance ℎ = 𝑐𝑡 1 ≈ 0.29𝑚).

Fig. 3 -

 3 Fig. 3-10:Pressure wave at various time instances along the depth of the container.

Fig. 3 - 11 :

 311 Fig. 3-11: Schematic illustration of stability regions.

Fig. 3 - 12 := 5 × 10 - 3

 3125103 Fig. 3-12: Time evolution of separation of phases (spinodal decomposition) for Cahn-Hilliard Eq. (2.2.29) without advection.

Fig. 4 - 1 :

 41 Fig. 4-1: Experimental observation of boundary layer fingering with supercritical hydrogen in weightlessness [81]. Here double arrow denotes the direction of vibration.

Fig. 4 - 2 :

 42 Fig. 4-2: Evolution of flow pattern with two miscible fluids subjected to vibration in horizontal direction in microgravity experiment 𝑓 = 8 𝐻𝑧, 𝐴 = 10.7 𝑚𝑚 at (a) t =1:4 s (b) 1.8 s (c) 2.9 s (d) 3.025 s. The Faraday waves (secondary instabilities) can be seen on the columns in (c) and (d). Reproduced from [88].

  1 . The system is analyzed for a wide range of vibration frequencies, varying from 5 𝐻𝑧 to 35 𝐻𝑧 while the proximity to the critical point and the quench rate, 𝑞 (= 𝛿𝑇 𝑇-𝑇 𝑐 × 100) are varied from 5 𝑚𝐾 to 2000

Fig. 4 -

 4 3 shows the schematic of the problem wherein a 2-D square cell (size ℎ = 7𝑚𝑚 ) with solid walls is filled with supercritical 𝐻 2 , implying an isochoric process. The dimensions of the cell are the same as considered in the work of Gandikota [89] .

Fig. 4 - 3 :

 43 Fig. 4-3: Schematic of problem for thermo-vibrational instabilities.

2 Fig. 4 - 4 :

 244 Fig. 4-4: Comparison between (a) numerical simulation and (b) experimental [89] observations for thermovibrational instabilities (𝑓 = 20 𝐻𝑧 and 𝐴 = 0.875 𝑚𝑚).

Fig. 4 - 5 :

 45 Fig. 4-5: Comparison of wavelengths, experimental [89] and numerical, for different proximities to the critical point for same 𝛿𝑇 = 400 𝑚𝐾, 𝐴 = 0.6125 𝑚𝑚 and 𝑓 = 30 𝐻𝑧.

Fig. 4 - 8 :

 48 Fig. 4-8: Critical amplitude (𝐴 𝑐𝑟 ) for the onset of Rayleigh-vibrational instabilities as a function of frequency and (a) quench rate for three proximities to the critical point 𝑇 𝑖 -𝑇 𝑐 = 2000𝑚𝐾, 500 𝑚𝐾 and 5 𝑚𝐾 (b) distance from the critical point for two different quench temperatures, 𝛿𝑇 = 5 𝑚𝐾 and 2.5 𝑚𝐾.

Fig. 4 -

 4 Fig. 4-8(a) shows the critical amplitude (𝐴 𝑐𝑟 ) as function of frequency (𝑓) and different quench rates (𝑞) for 𝑇 𝑖 -𝑇 𝑐 = 2000𝑚𝐾, 500 𝑚𝐾 and 5 𝑚𝐾. In order to ascertain the critical

Fig. 4 -

 4 8(b) shows the critical amplitude as a function of proximity to the critical point for the same quench temperature (𝛿𝑇). It is important to mention here that the comparison based on the same quench rate (𝑞) is not reasonable and may lead to erroneous interpretations that critical amplitude increases on approaching the critical point as in Fig.4-8(a). This is because, a different quench rate for different proximities to the critical point will induce a different value of 𝛿𝑇 in 𝑅𝑎 𝑣 and thus two independent parameters (𝛿𝑇 and 𝐴) in 𝑅𝑎 𝑣 . Thus, it is more judicious to compare the critical amplitude for the same temperature quench (𝛿𝑇) as shown in Fig.4-8(b) wherein a decrease in the value of 𝐴 𝑐𝑟 is observed on approaching the critical point. The diverging behavior of the expansivity results in a greater density gradient while diminishing diffusive forces (viscosity and thermal diffusivity) weakens the stabilizing effect. Consequently, the destabilization of the TBL occurs at a lower amplitude for the same vibration frequency as observed in Fig.4-8(b).

  of time and owing to the action of diffusive forces, adjacent fingers merge resulting in the change of the wavelength over a longer duration of time. This is illustrated in Fig.4-9(a-d) showing instabilities at various time instances for 𝑇 𝑖 -𝑇 𝑐 = 100 𝑚𝐾, 𝛿𝑇 = 10𝑚𝐾, 𝑓 = 20 𝐻𝑧 and 𝐴 = 4 𝑚𝑚.

Fig. 4 -

 4 Fig. 4-9(a-d): Change in wavelength in Rayleigh-vibrational instabilities with time for 𝑇 𝑖 -𝑇 𝑐 = 100 𝑚𝐾, 𝛿𝑇 = 10 𝑚𝐾, 𝑓 = 10 𝐻𝑧, 𝐴 = 4 𝑚𝑚.

Fig. 4 -

 4 Fig. 4-11(a-b) shows wavelengths as a function of acceleration for 𝑇 𝑖 -𝑇 𝑐 = 500 𝑚𝐾 and 𝑇 𝑖 -𝑇 𝑐 = 5 𝑚𝐾 for different quench rates. The dashed lines represent a power law fit to the measured wavelengths. It can be seen that the wavelength decreases (𝑖. 𝑒. the number of fingers increases) with the increase in quench rate which can be explained as follows.

Fig. 4 -

 4 Fig. 4-10 (a-d): Vector plot illustrating coalescence of fingers (marked in Fig. 4-9) with time ( 𝑇 𝑖 -𝑇 𝑐 = 100 𝑚𝐾, 𝛿𝑇 = 10 𝑚𝐾, 𝑓 = 10 𝐻𝑧, 𝐴 = 4 𝑚𝑚)

Fig. 4 -

 4 Fig. 4-11: Wavelength in Rayleigh-vibrational instabilities for different quench rates as function of acceleration for (a) 𝑇 𝑖 -𝑇 𝑐 = 500 𝑚𝐾 (b) 𝑇 𝑖 -𝑇 𝑐 = 5𝑚𝐾.

Fig. 4 -

 4 Fig. 4-12: Wavelength in Rayleigh-vibrational instabilities for 𝑇 𝑖 -𝑇 𝑐 = 500 𝑚𝐾 showing anomalous trend at higher accelerations (greater than ones considered in Fig. 4-11(a)) for 𝑞 = 10 % and 50 %.

Fig. 4 -

 4 Fig. 4-13: Wavelength in Rayleigh-vibrational instabilities for 𝑇 𝑖 -𝑇 𝑐 = 500 𝑚𝐾, 𝑞 = 10% for 2 different cell sizes.

Fig. 4 - 14 :

 414 Fig. 4-14: Wavelength in Rayleigh-vibrational instabilities with acceleration for different proximities to the critical point (a) 𝛿𝑇 = 100 𝑚𝐾 (b) 𝛿𝑇 = 5 𝑚𝐾.

Fig. 4 -

 4 Fig. 4-15: Critical Rayleigh-vibrational number with proximity to the critical point for 𝑓 = 10 𝐻𝑧 and 30 𝐻𝑧 on log-log plot.

Fig. 4 -

 4 Fig. 4-16: Parametric instabilities for one time-period of vibration for 𝑇 𝑖 -𝑇 𝑐 = 500 𝑚𝐾, 𝛿𝑇 = 50 𝑚𝐾, 𝑓 = 10 𝐻𝑧, 𝐴 = 10 𝑚𝑚.

Fig. 4 -

 4 Fig. 4-17: Schematic illustration of the mechanism leading to parametric instabilities (a) formation of thin TBL (b) equivalent configuration representing parametric (Faraday) instability like configuration (c) equivalent configuration representing stable arrangement 𝑤. 𝑟. 𝑡 Rayleigh-Taylor instability.

Fig. 4 -

 4 Fig. 4-18: 𝑥 -velocity in the bulk for 𝑇 𝑖 -𝑇 𝑐 = 500 𝑚𝐾, 𝛿𝑇 = 50 𝑚𝐾, 𝑓 = 10 𝐻𝑧, 𝐴 = 10 𝑚𝑚 for one time-period of vibration.

  -18 which shows velocity along 𝑥 -direction at 𝑦 = 3.5 𝑚𝑚 for one time-period of vibration. Here, 𝑥 -velocity profiles have been shown only in the bulk region (1 𝑚𝑚 ≤ 𝑥 ≤ 6mm) for the sake of clarity. The negative values are observed for 𝑡 = 2.65 𝑠 and 2.75 𝑠 which represent the movement of cell towards negative 𝑥 -direction at time instances when instabilities are observed near the left wall 𝐴𝐷 (Fig. 4-16(a) and (e)

Fig. 4 -

 4 Fig. 4-19: Parametric instabilities near left wall 𝐴𝐷 for 𝑇 𝑖 -𝑇 𝑐 = 500 𝑚𝐾, 𝛿𝑇 = 50 𝑚𝐾, 𝑓 = 20 𝐻𝑧, 𝐴 = 6 𝑚𝑚 with (a-c) symmetric and (d-f) wall conditions on the horizontal walls.

Fig. 4 -

 4 Fig. 4-20: Velocity vector near the corner (𝐴) with a wall condition on the top plate for 𝑇 𝑖 -𝑇 𝑐 = 500 𝑚𝐾, 𝛿𝑇 = 50 𝑚𝐾, 𝑓 = 20 𝐻𝑧, 𝐴 = 6 𝑚𝑚 illustrating (a) boundary layer velocity profile (b) formation of vortices (scale 1:8x10 5 cm/ms -1 ).

  -21(b) for 2 different quench temperatures 𝛿𝑇 = 100 𝑚𝐾 and 𝛿𝑇 = 5 𝑚𝐾 for initial temperature far (𝑇 𝑖 -𝑇 𝑐 = 2000 𝑚𝐾, 15000 𝑚𝐾 , 1000 𝑚𝐾 and 500 𝑚𝐾 ) and close (𝑇 𝑖 -𝑇 𝑐 = 100 𝑚𝐾, 50 𝑚𝐾 and 20 𝑚𝐾 ) from the critical point, respectively. As expected, the critical amplitude decreases on moving towards the critical point which follows the same argument as mentioned before.

Fig. 4 - 21 :

 421 Fig. 4-21: Critical amplitude (𝐴 𝑐𝑟 ) for the onset of parametric instabilities as a function of frequency and (a) quench rate for three proximities to the critical point 𝑇 𝑖 -𝑇 𝑐 = 2000𝑚𝐾, 500 𝑚𝐾 and 50 𝑚𝐾 (b) distance from the critical point for two different quench temperatures, 𝛿𝑇 = 100 𝑚𝐾 (for 𝑇 𝑖 -𝑇 𝑐 = 2000𝑚𝐾, 1500 𝑚𝐾, 1000 𝑚𝐾 and 500 𝑚𝐾) and 𝛿𝑇 = 2.5 𝑚𝐾 (for 𝑇 𝑖 -𝑇 𝑐 = 100𝑚𝐾, 50 𝑚𝐾and 20 𝑚𝐾).

Fig. 4 -

 4 Fig. 4-22: Critical amplitude for parametric instabilities for different cell sizes.

Fig. 4 -

 4 Fig. 4-23: Wavelength vs acceleration for 𝑇 𝑖 -𝑇 𝑐 = 500 𝑚𝐾 with different quench rates.

Fig. 4 -

 4 Fig. 4-24: Stability plot describing critical amplitude for onset of Rayleigh-vibrational and parametric instabilities as a function of quench rate and frequency.

4 .

 4 For the schematic of the problem as shown in Fig. 4-3, Fig. 5-1(a-b) shows these instabilities wherein temperature field has been shown for 𝑇 𝑖 -𝑇 𝑐 = 100 𝑚𝐾, 𝛿𝑇 = 10 𝑚𝐾, 𝑓 = 20 𝐻𝑧 and 𝐴 = 5 𝑚𝑚 at different time instances. A more interesting behavior can be observed in Fig. 5-2 which shows the profile of deviation of temperature from the wall (boundary) temperature at 𝑦 = 0.3 𝑚𝑚 for various time instances. It can be seen that this difference attains a negative value which subsequently grows with time.

Fig. 5 - 1 :

 51 Fig. 5-1: (𝑇 -𝑇 𝑤 ) contour plots for 𝑇 𝑖 -𝑇 𝑐 = 100 𝑚𝐾, 𝛿𝑇 = 10 𝑚𝐾, 𝑓 = 20 𝐻𝑧 and 𝐴 = 5 𝑚𝑚 at various time instances (a) 𝑡 = 7𝑠 (b) 𝑡 = 10𝑠. Insets highlight the growth of sink zones with time. The dashed lines in (a-b) represent the boundary where 𝑇 -𝑇 𝑤 = 0.

Fig. 5 - 2 :

 52 Fig. 5-2: Temperature profile at 𝑦 = 0.3 𝑚𝑚 at various time instances.

Fig. 5 -

 5 Fig. 5-3: : Temperature profile for 𝑇 𝑖 -𝑇 𝑐 = 100 𝑚𝐾, 𝛿𝑇 = 10𝑚𝐾, 𝑓 = 20 𝐻𝑧, 𝐴 = 5𝑚𝑚 at 10𝑠 (a) bottom wall adiabatic, sink-zones are formed (b) bottom wall isothermal (𝑖. 𝑒. at 𝑇 𝑖 ), no sink-zones are formed. The white dashed lines in (a) represent the boundary where 𝑇 -𝑇 𝑤 = 0.

Fig. 5 - 4 :

 54 Fig. 5-4: Temperature contour plots for 𝑇 𝑖 -𝑇 𝑐 = 100 𝑚𝐾, 𝛿𝑇 = 10𝑚𝐾, 𝑓 = 20 𝐻𝑧, 𝐴 = 5𝑚𝑚 at various instances of time along with streamlines depicting direction of motion of vibration. (a-b) for one-half period (c-d) growth of sink zone with time.

Fig. 5 -

 5 Fig. 5-5(a-b): Temperature contour plots illustrating sink zones (shaded regions) for various proximities to the critical point and other parameters as marked in the figure. The temperature contours are represented by lines.

Fig. 5 -

 5 Fig. 5-6: (a) Schematic of 1D case with both walls quenched (b) Evolution of the temperature field for 𝑇 𝑖 -𝑇 𝑐 = 20 𝑚𝐾, 𝛿𝑇 = 10𝑚𝐾, (solid lines): 100𝑔, (dotted lines): zero-g.

Fig. 5 -

 5 Fig. 5-6(b) illustrates the evolution of temperature field at various time instances for both the cases (𝑖)100𝑔 and (𝑖𝑖) zero-g as shown in Fig. 5-6(a). The results presented herein are coherent

Fig. 5 -

 5 Fig. 5-7: (a) Schematic of 1D case with bottom wall quenched and top adiabatic (b) Evolution of the temperature field for 𝑇 𝑖 -𝑇 𝑐 = 20 𝑚𝐾, 𝛿𝑇 = 10 𝑚𝐾, (solid lines): 100𝑔, (dotted lines): zero-g.

Fig. 5 -

 5 Fig. 5-7(b) and Fig. 5-8(b) illustrate the evolution of the temperature field for both cases shown in Fig. 5-7(a) and Fig. 5-8(a) respectively. It can be noticed that when the bottom wall is quenched (Fig. 5-7(b)) with top being adiabatic, the temperature at the top falls below the wall temperature. On the contrary, when the top wall is quenched with bottom wall being adiabatic (Fig. 5-8(b)), the temperature of the fluid at the bottom is higher than zero-g conditions and it does not

Fig. 5 - 8 :

 58 Fig. 5-8: (a) Schematic of 1D case with top wall quenched and bottom adiabatic (b): Evolution of the temperature field for 𝑇 𝑖 -𝑇 𝑐 = 20 𝑚𝐾, 𝛿𝑇 = 10 𝑚𝐾, (solid lines): 100 𝑔, (dotted lines): zero-g.

  Fig.5-7(b). The rest of explanation is similar to the case when both sides are quenched as mentioned earlier which leads to the drop of temperature near the top wall below the boundary temperature.

Fig. 5 - 9 :

 59 Fig. 5-9: Evolution of the temperature field (𝑇 -𝑇 𝑤 ) ((b) and (d)) for 𝑇 𝑖 -𝑇 𝑐 = 20 𝑚𝐾, 𝛿𝑇 = 10 𝑚𝐾, (solid lines) 100 𝑔 (dotted lines) zero-g for configurations shown in (a) and (c) respectively.

Fig. 5 -

 5 Fig. 5-12: Time of appearance of sink zones for the schematic as shown in Fig. 6(a) for a quench of 𝛿𝑇 = 10𝑚𝐾, (a) for different proximities to the critical point for an acceleration of 100 𝑔. (b) surface 3D plot for different accelerations and proximities to the critical point.

Fig. 5 -

 5 Fig. 5-13: Temperature profiles for 𝑇 𝑖 -𝑇 𝑐 = 20 𝑚𝐾 , 𝛿𝑇 = 2 𝑚𝐾, 𝑓 = 20 𝐻𝑧 and 𝐴 = 3 𝑚𝑚 for one-time period (50 𝑚𝑠). Here walls 𝐴𝐵 and 𝐶𝐷 are quenched. (a), (c), (e) temperature contour plots and (b), (d), (f) 1D plot of temperature field as a function of 𝑥 at 𝑦 = 0.3 𝑚𝑚 and 𝑦 = 1 𝑚𝑚. (a-b) 𝑡 = 14 𝑠 (c-d) 𝑡 = 14.025 𝑠 (e-f) 𝑡 = 14.05 𝑠.

Fig. 5 - 14 :

 514 Fig. 5-14: Temperature contours illustrating the decrease in the relative thickness of TBL (for example, marked approximately at two positions by 𝛿 𝑇𝐵𝐿,1 , 𝛿 𝑇𝐵𝐿,2 ) along the direction of vibration for 𝑇 𝑖 -𝑇 𝑐 = 20 𝑚𝐾 , 𝛿𝑇 = 2 𝑚𝐾, 𝑓 = 20 Hz and 𝐴 = 3 𝑚𝑚 (a) 14 s (b) 14.025 s.

Fig. 5 -

 5 Fig.5-15(a) shows the evolution of temperature difference at points 𝑃1 and 𝑃2 (as marked in the Figure), denoted by the variable 𝑇 𝑃1 -𝑇 𝑃2 for different amplitudes for 𝑇 𝑖 -𝑇 𝑐 = 20𝑚𝐾, 𝛿𝑇 = 2 𝑚𝐾, 𝑓 = 20 𝐻𝑧. In order to effectively capture the physics, a time step of 10 -6 𝑠 is used in this analysis (acoustic time is of ~ 30 𝜇𝑠 ). A nearly similar oscillatory behavior in the temperature field due to vibration has also been reported by Jounet et al.[START_REF] Jounet | Low-frequency vibrations in a near-critical fluid[END_REF]. While the difference in maximum temperature (∆𝑇 𝑝𝑒𝑎𝑘 ) is found to increase with amplitude, it is also observed that the maximum change in temperature (∆𝑇 𝑚𝑎𝑥 ) is not the same for fluid elements along the direction of vibration as shown for various 𝑥 positions at 𝑦 = 3.5 mm in Fig.5-15(b). The maximum change occurs near the wall while it decreases on moving towards the center. This peculiar observation has not been reported. These two observations can be elucidated as follows. The inertial force due to the vibrational acceleration cause the fluid to move in the direction of vibration. When the cell moves in positive 𝑥 -direction, the fluid element in the vicinity of the wall 𝐵𝐶 impacts the wall and is compressed by the resulting reaction force by virtue of high compressibility of the SCF. This compression causes an increase in the pressure in that region. However, the fluid element near the wall 𝐴𝐷 keeps moving towards the bulk fluid due to inertia and thus expands causing a drop in the pressure in that region. Consequently, a high value of isothermal compressibility (𝜒 𝑇 ) in conjunction with thermal expansion (𝛽 𝑃 ) causes a rise/drop of temperature in regions of compression/expansion as can be seen in Fig.5-15(a). The initial negative difference being due to

Fig. 5 -

 5 Fig. 5-15: (a) Difference in temperature at 𝑃1(0.5𝑚𝑚, 3.5𝑚𝑚) and 𝑃2 (6.5𝑚𝑚, 3.5𝑚𝑚) for different amplitudes of vibration for the schematic as shown on right. (b) Temperature at 𝑃2, 𝑃3( 5 𝑚𝑚, 3.5 𝑚𝑚). Plots are for 𝑇 𝑖 -𝑇 𝑐 = 20 𝑚𝐾, 𝑓 = 20 𝐻𝑧 and 𝐴 = 3 𝑚𝑚.

(

  ∆𝑇 𝑚𝑎𝑥 ) at two different 𝑥 positions (𝑃2 and 𝑃3, defined by 𝑟 = ∆𝑇 𝑚𝑎𝑥,𝑃2 ∆𝑇 𝑚𝑎𝑥,𝑃3,with same vibration parameters (𝑓 = 20 𝐻𝑧 and 𝐴 = 3 𝑚𝑚) is compared for two different proximities to the critical point (𝑇 𝑖 -𝑇 𝑐 = 500𝑚𝐾 and 20 𝑚𝐾). Since the compressibility decreases on moving away from the critical point thereby implying a lower damping-like action by the fluid (or a higher transmittance of force to the center), it is reasonable to ascertain that the difference between maximum change in temperature (∆𝑇 𝑚𝑎𝑥 ) near the wall and the center will be lower when far from the critical point. A higher value of 𝑟 for 𝑇 𝑖 -𝑇 𝑐 = 500𝑚𝐾 ( 𝑟 500 = 1.99997) as compared to 𝑟 20 = 1.99982 for 𝑇 𝑖 -𝑇 𝑐 = 20𝑚𝐾 thus substantiates the above explanation. In order to further illustrate the effect of compressibility, Fig.5-16 shows the evolution of difference in temperature at 𝑃1 and 𝑃2 for different proximities to the critical point for the same vibration parameters. It can be observed that on moving away from the critical point, the effect of FPE (maximum change in temperature due to vibration) reduces by virtue of decreasing compressibility. Thus, it can be remarked that change in temperature due to FPE in vibration therefore depends on the compressibility as well as external forces. A higher value of both these factors will entail a larger change in temperature.

Fig. 5 - 16 :

 516 Fig. 5-16: Difference in temperature at 𝑃1 and 𝑃2 for different proximities to the critical point for the schematic as shown on right for 𝑓 = 20 𝐻𝑧 and 𝐴 = 3 𝑚𝑚.

Fig. 5 -

 5 Fig. 5-17(a) where the line '𝑀𝑁' represents the thermal boundary layer at a time instant, say 𝑡 0 .The points marked by 𝐶1 and 𝐶2 represent the fluid elements near the left (𝐴𝐷) and right (𝐵𝐶) walls (see Fig.4-3), respectively, and are symmetrically located with respect to 𝑥 = 3.5 𝑚𝑚. The schematic of the temperature profiles along 𝑦 -direction is further illustrated in Fig.5-17(b-c). In the absence of vibration (𝑖. 𝑒. action of piston effect only) and only top and bottom wall being quenched, it can be well understood that the same temperature will persist at points 𝐶1 and 𝐶2

Fig. 5 - 17 :

 517 Fig. 5-17: Schematic illustration of the see-saw motion of thermal boundary layer (a) actual (𝑀𝑁) and relative thickness (𝑀𝑁 𝑥+ and 𝑀𝑁 𝑥-) of the TBL (b) (i) temperature profile along 𝑦 axis in the presence of piston effect only (ii) evolution of temperature profile in the presence of vibration for vibration in positive 𝑥-direction.

  -17(a) by 𝑀𝑁 𝑥+ (𝑀𝑁 𝑥-for negative 𝑥 -direction) This aspect can be further elucidated from Fig. 5-18 wherein temperature contour plots and temperature profiles at various 𝑥 are shown for 𝑇 𝑖 -𝑇 𝑐 = 20𝑚𝐾, with 𝛿𝑇 = 2 𝑚𝐾, 𝑓 = 20 𝐻𝑧 and 𝐴 = 0.5 𝑚𝑚. These have been depicted along-side each other for clear illustration of the relative thickness even though the actual thickness remains the same. It can be clearly seen from Fig. 5-18

Fig. 5 -

 5 Fig. 5-18: (a, c) Temperature contours illustrating change in relative thickness of TBL over half-period of vibration (b, d) Temperature profiles at various 𝑥 along 𝑦-direction illustrating nearly the same actual thickness of the TBL, for 𝑇 𝑖 -𝑇 𝑐 = 20 𝑚𝐾 , 𝛿𝑇 = 2 𝑚𝐾, 𝑓 = 20 𝐻𝑧 and 𝐴 = 0.5 𝑚𝑚.

Fig. 5 - 2 𝑔= 1 9

 5219 Fig. 5-19(a) shows a 3D plot for the maximum change in temperature (∆𝑇 𝑚𝑎𝑥 ) as function of different accelerations and proximities to the critical point when the SCF is only subjected to vibrations, as obtained by 1D analysis for the schematic shown in Fig. 5-19(b). For the sake of clarification, a 2D plot for 𝐴𝜔 2 𝑔= 1 is also shown in the inset to in order to describe decreasing behavior of ∆𝑇 𝑚𝑎𝑥 with the proximity to the critical point. As expected, a higher temperature rise is observed with increasing acceleration and decreasing proximity to the critical point. In addition, an increase in slope (obtained by best linear fit at a fixed acceleration) from -0.009 for 𝐴𝜔 2 𝑔 = 1 to

Fig. 5 -

 5 Fig. 5-19: (a) Maximum change in temperature, when SCF is subjected to vibration only, as a function of various proximities to the critical point and acceleration (𝐴𝜔 2 ) with inset illustrating the decrease of ∆𝑇 𝑚𝑎𝑥 for 𝐴𝜔 2 𝑔 = 1 (b) schematic of the problem used to evaluate maximum change in temperature.

Fig. 6 - 3 :

 63 Fig. 6-3: Maximum velocity as a function of time for 𝐶𝑛 = 0.02 on a 200 2 mesh for 𝑎 = 1𝑚𝑚 (1 𝐾 below the critical point).

Fig. 6 - 4 :

 64 Fig. 6-4: Velocity contours illustrating spurious velocity fields for a=1 mm (a-b) 𝑥 -velocity and (c-d) 𝑦 -velocity at various times (𝐶𝑛 = 0.01, 𝑀 0 = 10 -13 S.I).

Fig. 6 - 5 : 138 139Fig. 6 - 6 :

 6513866 Fig. 6-5: Schematic of test case for transformation of elliptical vapor bubble to circular bubble in a square cavity.

Fig. 6 - 7 :

 67 Fig. 6-7: Schematic of the test case for the coalescence of two liquid drops suspended in square cavity filled with vapor.

Fig. 6 - 8 :

 68 Fig. 6-8: Time evolution of coalescence of two liquid droplets just in contact with each other on 200 2 mesh (𝐶𝑛 = 0.02, 𝑀 0 = 10 -15 S.I).

( 2

 2 𝑚𝑚 × 2𝑚𝑚) filled at the critical density of CO2, 𝜌 𝑐 = 467.7 kg/m 3 at 𝑇 = 303.2 𝐾 (0.93 𝐾 below the critical point) with random perturbation added to the density. Thus, the initial conditions for 𝜙 thus correspond to, 𝜙((𝑥, 𝑦),0) = 𝜙 0 + 𝜙 𝑟𝑎𝑛𝑑(6.4.1)

Fig. 6 -

 6 9 schematically illustrates the thermodynamic representation of the initial state and densities of the liquid and vapor phase expected after the separation. The results from numerical simulation are shown in Fig. 6-10 over a period of time where contours are depicted for the density field highlighting the separation of liquid and vapor phase for 𝐶𝑛 = 0.025 and mobility as 5 × 10 -7 S.I. It can be observed that the mixture separates into liquid and gas phase with densities as sobered from contour plots in Fig. 6-10 (e-h).

Fig. 6 - 9 : 145 146Fig. 6 -

 691456 Fig. 6-9: Schematic illustration of initial conditions of the mixture and expected densities after separation of phases.

• 2 )

 2 time step in relation with the mobility in the Cahn-Hilliard equation which needs to be investigated. It will be interesting to analyze the behavior of fluid near the critical point in sub-critical regime when subjected to simultaneous quench and vibration which can help to gain insight into the several intriguing flow features ascribed to the anomalous behavior of thermo-physical properties near the critical point. B. Dispersion relation example Consider a 1-D wave equation given by, We can write solution 𝑢(𝑥, 𝑡) at any position 𝑥 and time 𝑡 as 𝑢 ̃𝑒𝑖(𝑘𝑥-𝜔𝑡) summed over all wavenumbers and frequencies. In order to derive the relation between 𝑘 and 𝜔 , we substitute the solution in Eq. (B.2) and we get ∂(𝑢 ̃𝑒𝑖(𝑘𝑥-𝜔𝑡) ̃𝑒𝑖(𝑘𝑥-𝜔𝑡) = -𝑐𝑖𝑘𝑢 ̃𝑒𝑖(𝑘𝑥-𝜔𝑡) ⇒ 𝜔 = 𝑐𝑘 (B.4) Eq. (B.4) thus represents relation between 𝑘 and 𝜔 for Eq. (1). We find that this is a constant linear relation. Similarly, dispersion relation can be derived for a given numerical scheme (as illustrated in (Appendix C)).

1 )

 1 order illustrate the nature of error propagation, consider a 1D wave equation, moving towards right with sound velocity 𝑐, mathematically represented by, Let the solution at any time 𝑡 𝑛 (= 𝑛∆𝑡) and position 𝑥 𝑗 be represented by, 𝑢 𝑗 𝑛 = ∫ 𝑈(𝑘∆𝑥, 𝑛∆𝑡)𝑒 𝑖𝑘𝑥 𝑗 𝑑𝑘 (C.2)

  𝑗 𝑛 = ∫ 𝑈 0 (𝑘) [|𝐺 𝑗 |] 𝑛 𝑒 𝑖(𝑘𝑥 𝑗 -𝑛𝛽 𝑗 ) 𝑑𝑘 = ∫ 𝑈 0 (𝑘) [|𝐺 𝑗 |] 𝑛 𝑒 𝑖𝑘(𝑥 𝑗 -𝑐 𝑁 𝑡 𝑛 ) 𝑑𝑘 (C.3) More formally, 𝐺 𝑗 is a complex number described as 𝐺 𝑗 = 𝐺| (𝑥=𝑥 𝑗 ) = |𝐺 𝑗 |𝑒 -𝑖𝛽 𝑗 with 𝑡𝑎𝑛 𝛽 𝑗 = -𝐺 𝐼𝑗 𝐺 𝑅𝑗 and 𝐺 𝐼𝑗 , 𝐺 𝑅𝑗 are imaginary and real parts of 𝐺 𝑗 , respectively. The modulus of error growth is |𝐺 𝑗 | = √(𝐺 𝐼𝑗 ) 2 + (𝐺 𝑅𝑗 ) 2 . Here, analogous to physical dispersion relation, 𝜔 = 𝑐𝑘, we have the numerical dispersion relation 𝜔 𝑁 = 𝑐 𝑁 𝑘, which gives 𝑐 𝑁 = 𝛽 𝑗 𝑘∆𝑡 . The error is defined as: 𝑒 = 𝑢(𝑥, 𝑡) -𝑢 𝑁 (𝑥, 𝑡), with 𝑢(𝑥, 𝑡) and 𝑢 𝑁 (𝑥, 𝑡) representing the exact and numerical solutions, respectively. Thus, we can write the spatial and temporal derivatives of 𝑢 𝑁 (𝑥, 𝑡) at any position 𝑥 and time instant 𝑡 using Eq. (C.3) as, (the subscripts 𝑗 and n have been dropped for clarity in the following equations) (𝑥-𝑐 𝑁 𝑡) 𝑑𝑘 (C.5) 157 Multiplying Eq. (C.4) by 𝑐 𝑁 and adding to Eq (C.5) we get, (𝑥-𝑐 𝑁 𝑡) 𝑑𝑘 + 𝑐 𝑁 ∫ 𝑖𝑘 𝑈 0 [|𝐺|] 𝑡 ∆𝑡 𝑒 𝑖𝑘(𝑥-𝑐 𝑁 𝑡) 𝑑𝑘 (C.6) Subtracting Eq. (C.6) from Eq. (C.1),we get, (𝑥-𝑐 𝑁 𝑡) 𝑑𝑘 -𝑐 𝑁 ∫ 𝑖𝑘 𝑈 0 [|𝐺|] 𝑡 ∆𝑡 𝑒 𝑖𝑘(𝑥-𝑐 𝑁 𝑡) 𝑑𝑘 (C.7) Now, using integration by parts on first term of Eq. (C.7), we get, 𝜕𝑒 𝜕𝑡 + (𝑐 -𝑐 𝑁 ) 𝜕𝑒 𝜕𝑥 = 𝑐 𝑁 ∫ 𝑖𝑘 𝑈 0 (𝑘) [|𝐺|] 𝑡 ∆𝑡 𝑒 𝑖𝑘(𝑥-𝑐 𝑁 𝑡) 𝑑𝑘 -′ (𝑥-𝑐 𝑁 𝑡) 𝑑𝑘 ′ ] 𝑑𝑘 -∫ 𝐿𝑛|𝐺| ∆𝑡 𝑈 0 (𝑘)[|𝐺|] 𝑡 ∆𝑡 𝑒 𝑖𝑘(𝑥-𝑐 𝑁 𝑡) 𝑑𝑘 -𝑐 𝑁 ∫ 𝑖𝑘 𝑈 0 [|𝐺|] 𝑡 ∆𝑡 𝑒 𝑖𝑘(𝑥-𝑐 𝑁 𝑡) 𝑑𝑘 (C.8)Cancelling the first and fourth term, we end up with error relation as, (𝑥-𝑐 𝑁 𝑡) 𝑑𝑘 (C.9)

2 (D. 11 ) 1 [Further)

 2111 solution any time 𝑡 𝑛 (= 𝑛∆𝑡) and position 𝑥 𝑗 , in the form of Eq. (C.2) can be written as 𝑢 𝑗 𝑛 = ∫ 𝑈(𝑘∆𝑥, 𝑛∆𝑡)𝑒 𝑖𝑘𝑥 𝑗 𝑑𝑘 . Further, we can define,𝜕𝑢 𝑗 𝑛 𝜕𝑥 = ∫ 𝑖𝑘𝑈(𝑘∆𝑥, 𝑛∆𝑡)𝑒 𝑖𝑘𝑥 𝑗 𝑑𝑘 (D.2)We can also write 𝜕𝑢 𝑗 𝑛 𝜕𝑥 at any 𝑥 𝑗 in terms of 𝑢 𝑗 𝑛 defined on other grid points as Here 𝐶 is the matrix obtained by writing the derivative in terms of 𝑢 𝑗 𝑛 defined on grid points. We can thus rewrite Eq. (D.1) as, 𝑢 𝑗 𝑛+1 -𝑢 𝑗 𝑛 = -𝑁 𝑐 [𝐶][𝑢 𝑗 𝑛+1 ] (D.3) Eq. (D.3) can now be further developed as follows, ∫ 𝑈(𝑘∆𝑥, (𝑛 + 1)∆𝑡)𝑒 𝑖𝑘𝑥 𝑗 𝑑𝑘 -∫ 𝑈(𝑘∆𝑥, 𝑛∆𝑡)𝑒 𝑖𝑘𝑥 𝑗 𝑑𝑘 = -𝑁 𝑐 ∫ ∑ 𝐶 𝑗𝑙 𝑀 𝑙=1 𝑈(𝑘∆𝑥, (𝑛 + 1)∆𝑡) 𝑒 𝑖𝑘𝑥 𝑙 𝑑𝑘 (D.4) which can be simplified to, ∫[𝑈 (𝑘∆𝑥, (𝑛 + 1)∆𝑡)𝑒 𝑖𝑘𝑥 𝑗 -𝑈(𝑘∆𝑥, 𝑛∆𝑡)𝑒 𝑖𝑘𝑥 𝑗 + 𝑁 𝑐 ∑ 𝐶 𝑗𝑙 𝑀 𝑙=1 𝑈(𝑘∆𝑥, (𝑛 + 1)∆𝑡) 𝑒 𝑖𝑘𝑥 𝑙 ]𝑑𝑘 = 0 (D.5) Since Eq. (D.5) is valid for all values of 𝑘, the integrand must be zero. Thus, we have, 𝑈(𝑘∆𝑥, (𝑛 + 1)∆𝑡)𝑒 𝑖𝑘𝑥 𝑗 -𝑈(𝑘∆𝑥, 𝑛∆𝑡)𝑒 𝑖𝑘𝑥 𝑗 = -𝑁 𝑐 ∑ 𝐶 𝑗𝑙 𝑀 𝑙=1 𝑈(𝑘∆𝑥, (𝑛 + 1)∆𝑡) 𝑒 𝑖𝑘𝑥 𝑙 (D.6) By dividing Eq. (D.6) by 𝑈(𝑘∆𝑥, 𝑛∆𝑡)𝑒 𝑖𝑘𝑥 𝑗 , it further gives, 𝐺 𝑗 -1 = -𝑁 𝑐 ∑ 𝐶 𝑗𝑙 𝑀 𝑙=1 𝐺 𝑗 𝑒 𝑖𝑘(𝑥 𝑙 -𝑥 𝑗 ) (D.7) Now, for central difference scheme, we can write, can write Eq. (D.7) for any 𝑗 as, 𝐺 𝑗 -1 = -𝑁 𝑐 𝐺 𝑗 [ 𝑒 𝑖𝑘∆𝑥 -𝑒 -𝑖𝑘∆𝑥 2∆𝑥 ] (D.10) After simplifications, Eq.(D.10) gives the expression of 𝐺 𝑗 , 𝐺 𝑗 = 1-𝑖𝑁 𝑐 sin (𝑘∆𝑥) 1+(𝑁 𝑐 sin (𝑘∆𝑥))Here 𝐺 𝑗 is complex with 𝐺 𝑅𝑗 = 1 1+(𝑁 𝑐 sin (𝑘∆𝑥)) 2 , 𝐺 𝐼𝑗 = -𝑖𝑁 𝑐 sin(𝑘∆𝑥) 1+(𝑁 𝑐 sin (𝑘∆𝑥)) 2 and |𝐺 𝑗 | = 1+(𝑁 𝑐 sin(𝑘∆𝑥)) 2 ] 1/2 . can now be evaluated using expressions obtained in Eq. (D.11). We know that 𝑉 𝑔𝑁 = 𝑑𝜔 𝑁 𝑑𝑘 , where 𝜔 𝑁 = 𝛽 𝑗 ∆𝑡 . Thus, using expression for 𝛽 𝑗 from Eqs. (D.11) and (D.12) we get, ∆𝑥 tan -1 (𝑁 𝑐 sin(𝑘∆𝑥)) (D.15)

5 ) 6 ).

 56 𝑇 = ∫ 𝐹 𝑉 𝑑𝑉 = ∫ 𝑓( 𝜙, 𝜵𝜙, 𝜵 2 𝜙 … . )𝑑𝑉 (E.1)Here 𝐹 𝑉 is free energy per unit volume, denoted by function 𝑓, which can be expanded using Taylor series expansion in terms of 𝜙 and 𝜵𝜙 as, 𝑓 (𝜙, 𝜵𝜙, 𝜵 2 𝜙. . ) = 𝑓 0 (𝜙) , it is required that all the odd rank tensors be identically equal to zero (𝑖. 𝑒. 𝛼 𝑖 = 0). Thus, we get, 𝐹 𝑇 = ∫(𝑓 0 (𝜙) + 𝛽 𝑖𝑗 𝜵 2 𝜙 + 𝛾 𝑖𝑗 (𝜵𝜙). (𝜵𝜙)+. . )𝑑𝑉 (E.4)Now, considering ∫ 𝛽 𝑖𝑗 𝜵 2 𝜙 𝑑𝑉, we have∫ 𝛽 𝑖𝑗 𝜵 2 𝜙 𝑑𝑉 = ∫ 𝛽 𝑖𝑗The first term in (E.5) can be reduced to ∫(𝛽 𝑖𝑗 𝜵𝜙) 𝑗 . 𝒏𝑑𝐴 which is equal to zero due to boundary condition that intensive property is conserved. Thus, we have,∫ 𝛽 𝑖𝑗 𝜵 2 𝜙 𝑑𝑉 = -The expression (E.4) thus reduces to, 𝐹 𝑇 = ∫ (𝑓 0 (𝜙) + ( 𝑓 0 (𝜙) + 𝜅 𝑖𝑗 (𝜵𝜙) 𝑖 (𝜵𝜙) 𝑗 )𝑑𝑉 (Owing to symmetry conditions, we have 𝜅 𝑖𝑗 = 𝜅 𝑗𝑖 and thus the above expression (E.7) can be reduced to 𝐹 𝑇 = ∫ (𝑓 0 (𝜙) + 𝜅(𝜵𝜙) 2 ) 𝑑𝑉 (E.8)

. 2 |𝜵𝜙| 2 ,

 22 The equilibrium will be attained when 𝜂 = 0. For 𝐹 𝑉 = 𝛽 ̅ 𝑓 0 (𝜙) + 𝛼 ̅ we get,𝜂 = 𝛽 ̅ 𝑓 𝑜 ′(𝜙) -𝜵. (𝛼 ̅ 𝜵𝜙) (E.10)which, for a constant 𝛼 ̅ can be reduced to,𝜂 = 𝛽 ̅ 𝑓 𝑜 ′(𝜙) -𝛼 ̅𝜵 2 𝜙 (E.11)Here, the variational derivative relation ( for a functional has been used. At equilibrium, we have 𝜂 = 0 which implies,𝛽 ̅ 𝑓 𝑜 ′(𝜙) = 𝛼 ̅ 𝜵 2 𝜙 (E.12)Considering a 1D case, we have, (writing 𝑓 𝑜 (𝜙) as 𝑓 𝑜 for simplicity),

3 ) 6 ) 2 𝜵(|𝜵𝜙| 2 )

 3622 written as (owing to conservation of 𝜙 over the entire domain, i.e.) 𝐸 𝑐𝑜𝑛𝑣 = -∫ 𝜂 𝜵. (𝑽𝜙)𝑑𝑉 = -∫ 𝜂 ∑ 𝜕(𝑉 𝑗 𝜙) 𝜕𝑥 𝑗 𝑗 𝑑𝑉 which can be further modified to,𝐸 𝑐𝑜𝑛𝑣 = -∫ ∑ 𝜕(𝜂𝑉 𝑗 𝜙)The first term first term in (F.2) is zero due to boundary condition as, (𝜂𝜙𝑽)𝑑𝑉 = ∫ 𝜂𝑽𝜙. 𝒏𝑑𝐴 = 0. Now, the change in energy (𝐸 𝑑 ) due to the diffusive force, (𝑓 𝑆𝑇 ) acting at the interface will be equal to its work done and can be represented by,𝐸 𝑑 = ∫ ∑ 𝑓 𝑆𝑇 𝑗 𝑉 𝑗 𝑗 𝑑𝑉 (F.4)As mentioned before, the change in energy due to convection in will be countered by work done by diffusive force, 𝑖. 𝑒. 𝐸 𝑐𝑜𝑛𝑣 = -𝐸 𝑑 . Thus comparing (F.3) and (F.4) we get,The above expression holds true for a compressible fluid and is in potential form. This can be converted to stress form as follows. We have,𝑓 𝑆𝑇 = -𝜙𝜵𝜂 = -𝜵 (𝜙𝜂) + 𝜂𝜵𝜙(substituting for 𝜂 = 𝛽 ̅ 𝑓 𝑜 ′(𝜙) -𝛼 ̅ 𝜵 𝟐 𝜙 ), we have𝑓 𝑆𝑇 = -𝜵 (𝜙𝜂) + 𝛽 ̅ 𝑓 𝑜 ′(𝜙)𝜵𝜙 -𝛼 ̅ 𝜵 𝟐 𝜙𝜵𝜙 (F.7)Now, we have 𝜵 (𝛽 ̅ 𝑓 𝑜 (𝜙)) = 𝛽 ̅ 𝑓 𝑜 ′(𝜙)𝜵𝜙 following which (F.7) reduces to,𝑓 𝑆𝑇 = -𝜵 (𝜙𝜂) + 𝜵 (𝛽 ̅ 𝑓 𝑜 (𝜙)) -𝛼 ̅ 𝜵 2 𝜙𝜵𝜙 (F.8)Using the vector identity, 𝜵. (𝜵𝜙 ⊗ 𝜵𝜙) = 𝜵𝜙𝜵 2 𝜙 + 1 for the last term in (F.8), we get, 𝑓 𝑆𝑇 = -𝜵 (𝜙𝜂) + 𝜵 (𝛽 ̅ 𝑓 𝑜 (𝜙)) -𝛼 ̅ 𝜵. (𝜵𝜙 ⊗ 𝜵𝜙)

  𝑆𝑇 = -𝜵 (𝜙𝜂) + 𝜵 (𝛽 ̅ 𝑓 𝑜 (𝜙)) -𝛼 ̅ 2 𝜵(|𝜵𝜙| 2 -𝛼 ̅ 𝜵. (𝜵𝜙 ⊗ 𝜵𝜙) + 𝛼 ̅𝜵(|𝜵𝜙| 2 ) (F.10) Since, 𝛼 and 𝛽 are constant, (F.10) can be modified to, 𝑓 𝑆𝑇 = -𝜵 (𝜙𝜂 -𝛽 ̅ 𝑓 𝑜 (𝜙) + 𝛼 ̅ 2 𝜵(|𝜵𝜙| 2 ) -𝛼 ̅ 𝜵. (𝜵𝜙 ⊗ 𝜵𝜙) + 𝜵(𝛼 ̅|𝜵𝜙| 2 ) (F.11)

Mots clés: ECOULEMENT COMPRESSIBLE, FLUIDES CRITIQUES, EFFET PISTON, MODÉLISATION NUMÉRIQUE, CHAMP DE PHASE xii NOMENCLATURE

		Symbols
	𝑎	Radius	mm
	𝐴	Amplitude	mm
		Area	m 2
	𝐴 𝑐𝑟	Critical amplitude	mm
	𝑐	Wave speed	m.s -1
	𝑐 𝑠	Sound speed	m.s -1
	𝑐 𝑁	Phase speed	m.s -1
		1

Table 1 -

 1 1: Universal power laws for thermodynamic and correlation properties[START_REF] Zappoli | Heat Transfers and Related Effects in Supercritical Fluids[END_REF]

	Property				Power law	Critical Exponent
	Isothermal Compressibility	𝜀 -𝛾	𝛾 = 1.239
	(𝜒 𝑇 =	1 𝜌	( 𝜕𝜌 𝜕𝑃	)	𝑇	)
	Coefficient of thermal expansion	𝜀 -𝛾
	(𝛽 𝑃 = -	1 𝜌	( 𝜕𝜌 𝜕𝑇	) 𝑃	)
	Heat Capacity at constant pressure (𝐶 𝑃 )	𝜀 -𝛾
	Heat Capacity at constant volume (𝐶 𝑉 )	𝜀 -𝛼	𝛼 = 0.11
	Sound speed (𝑐 𝑠 )	𝜀 𝛼/2

  2.2.16) Here, 𝑽 is the mean velocity of the fluid. Substituting (2.2.16) in (2.2.14) we get, Now, considering mass fraction of fluid 1 as the phase-field variable thereby designating 𝛾 1 by 𝜙 (which implies 𝛾 2 = 1 -𝜙 ), (2.2.22) can be reduced to,

	𝜕𝛾 1 𝜕𝑡	+ 𝑽. 𝜵𝛾 1 =	1 𝜌	𝜵. (𝐽 𝑑1 )	(2.2.22)
	𝜕𝜙 𝜕𝑡	+ 𝑽. 𝜵𝜙 =	1 𝜌	𝜵. (𝐽 𝑑1 )	(2.2.23)
	𝜕𝜌 ̃1
	𝜕(𝛾 1 𝜌) 𝜕𝑡	+ 𝛻. (𝛾 1 𝜌𝑽) = 𝜵. (𝐽 𝑑1 )	(2.2.21)
	which can be simplified to non-conservative form,

𝜕𝑡 + 𝜵. (𝜌 ̃1𝑽) = 𝜵. (𝐽 𝑑1 )

(2.2.17)

A similar relation can be written for fluid 2 in terms of its diffusive flux (𝐽 𝑑2 ) which yields, 𝜕𝜌 ̃2 𝜕𝑡 + 𝜵. (𝜌 ̃2𝑽) = 𝜵. (𝐽 𝑑2 ) (2.2.18) Adding (2.2.17) and (2.2.18), and using the relation from (2.2.12) we get, 𝜕𝜌 𝜕𝑡 + 𝜵. (𝜌𝑽) = 𝜵. (𝐽 𝑑1 + 𝐽 𝑑2 ) (2.2.19) Now, as in the interface region, the total mass is conserved which means, 𝐽 𝑑1 + 𝐽 𝑑2 = 0. Using this relation, we get the continuity equation, 𝜕𝜌 𝜕𝑡 + 𝜵. (𝜌𝑽) = 0 (2.2.20)

where the mean velocity V is the barycentric velocity, given by 𝜌𝑽 = 𝜌 ̃1𝑽 𝟏 + 𝜌 ̃2𝑽 𝟐 which is in coherence with the choice of the parameter (mass-fraction) considered to derive this relation. The above result shows that by considering mass fraction, the continuity equation holds throughout the domain. In order to derive the phase-field equation, substituting for 𝜌 ̃1 = 𝛾 1 𝜌 in (2.2.17), we get,

Table 3 -

 3

1: Physical properties and dimensions under consideration

  Dispersion RelationIt is the relation between the wavenumber 𝑘 and the frequency 𝜔. This means that not all possible values of 𝑘 can be taken for all values of 𝜔. There is a constraint imposed by the PDE which results in a relation between 𝑘 and 𝜔 as 𝜔 = 𝜔(𝑘) known as the dispersion relation. The dispersion relation also helps to evaluate the group velocity which is defined as 𝑉 𝑔𝑁 =

	𝑑𝜔 𝑑𝑘	. A simple
	example of analytical dispersion relation for 1D wave equation is given in the Appendix B. For
	accurate and precise solution, it is desired that numerical dispersion relation be completely coherent
	with analytical dispersion relation.

𝑐 𝑁 > 𝑉 𝑔𝑁 : normal dispersive medium ▪ 𝑐 𝑁 < 𝑉 𝑔𝑁 : anomalous dispersive medium •

Table 4 -

 4 1: Thermo-physical properties of 𝑛 -𝐻 2 with 𝜀 =

	𝑇-𝑇 𝑐 𝑇 𝑐	[10, 82]

Table 6 -

 6 Schematic of test case (square cavity filled with liquid CO2 with vapor bubble) for the investigation of spurious currents.

		1: Thermo-physical properties of CO2 near its critical point [28]
	𝑇 -	Equilibrium	Density		Viscosity		Isothermal		Surface
	𝑇 𝑐 (K)	Pressure	[ kg.m -3 ]		[ 𝜇Pa-s]		compressibility	Tension
		[MPa]					[1/MPa]		[mN.m -1 ]
			𝜌 𝐿	𝜌 𝑉	𝜇 𝐿	𝜇 𝑉	𝜒 𝑇,𝐿	𝜒 𝑇,𝑉	𝜎
	0.1	7.3606	525.17	411.68	37.579 29.164	27.344	40.359	0.003296
	0.93	7.2219	590.99	347.28	43.541 25.291	1.449	3.8128	0.055113
	1	7.2104	594.22	344.25	43.857 25.123	1.3114	3.5507	0.060375

The first term represents the pressure, while the remaining two terms contribute to the stress form. Thus, in the stress form, the coupling is achieved by adding 𝛼 ̅ 𝜵. (𝜵𝜙 ⊗ 𝜵𝜙) + 𝜵(𝛼 ̅|𝜵𝜙| 2 ) to the to the momentum equation. The first term is classically known as Korteweg stress tensor.

I would also like to thank them for providing me with enough opportunities to present my work at conferences, especially CEFIPRA (Centre Franco-Indien pour la Promotion de la Recherche Avancee) workshop-2017. thanks to Madame Valerie Thouard, for helping in all the administrative stuff especially

CHAPTER 4. THERMO-VIBRATIONAL INSTABILITIES IN SUPERCRITICAL FLUID

Among several applications and interesting phenomena ascribed to intriguing properties of SCFs, their behavior when subjected to simultaneous quench and vibration has piqued the interest of several researchers [START_REF] Jounet | Low-frequency vibrations in a near-critical fluid[END_REF][START_REF] Amiroudine | Thermovibrational instability in supercritical fluids under weightlessness[END_REF][START_REF] Gandikota | Rayleigh and parametric thermo-vibrational instabilities in supercritical fluids under weightlessness[END_REF][START_REF] Carlès | The unexpected response of near-critical fluids to low-frequency vibrations[END_REF], especially in weightlessness. This has been primarily motivated by g-jitters experienced by space vehicles and cryogenic reservoirs aboard the International Space Station (ISS). Further, the complex thermo-mechanical coupling leads to the evolution of several intriguing behaviors in the thermal boundary layer (TBL) resulting in various thermo-vibrational instabilities.

In this chapter, contemporary work pertinent to thermo-vibrational instabilities in supercritical fluids is reviewed followed by primary investigation of two types of instabilities, Rayleighvibrational and parametric instabilities. The analysis provides an insight into the physical mechanism causing these instabilities and the effect of various parameters (such as amplitude of vibration, temperature quench, distance from the critical point etc.) on the onset of these instabilities and the behavior of waves/finger like patterns. A stability plot is presented describing the type of instabilities to be expected for a given thermo-vibrational parameters.

Literature review

In highly compressible fluids such as SCFs, the action of vibration induces flow ascribed primarily to the high compressibility and thus captivates features due to a strong thermo-mechanical coupling. One of the initial pioneering work was carried by Carles and Zappoli [START_REF] Carlès | The unexpected response of near-critical fluids to low-frequency vibrations[END_REF] who analyzed the problem of a container filled with a near-critical fluid when subjected to mechanical vibrations by means of matched asymptotic expansions and presented various characteristic regimes of the fluid under vibration. While the acoustic regime, defined for vibration frequencies greater than the inverse of acoustic times (typically a few kHz) was found to be independent of the initial proximity to the critical point, a specific low-frequency regime was highlighted wherein the bulk part of the fluid almost behaves like a solid bouncing back and forth between two highly compressible thermal boundary layers. This response is specific to near-critical fluids. This phenomenon was further pertinent only to each type of instability. Thus, for RV instabilities, the top (𝐴𝐵) and bottom (𝐶𝐷) wall are quenched while for parametric only the vertical walls (𝐵𝐶 and 𝐴𝐷) are quenched. The remaining walls in both cases are maintained at adiabatic conditions. Further, the system is investigated under the conditions of weightlessness 𝑖. 𝑒. zero gravity. Further, for numerical analysis, an exponential grid with a total of 300 x 300 elements is used with smallest element being 1 𝜇𝑚 near the wall with time step as 10 -3 𝑠 , which corresponds to a compromise between accuracy and CPU time. Further, the force term in the momentum equation (2.1.17) corresponds to 𝑭 = 𝜌𝐴𝜔 2 sin (𝜔𝑡)𝒊 in the 𝑥 -direction. The thermo-physical properties are evaluated using relations obtained from Renormalization-Group (RG) Theory as described in Table 4-1.

Comparison with experimental observations

In order to ensure that the developed mathematical model well captures the dynamics of SCF, the results obtained from simulations are compared with experimental observation when SCF is simultaneously subjected to thermal perturbation (quench) and mechanical vibrations. For the experimental comparison, all the four walls of the cell are quenched while the vibration parameters corresponding to experimental conditions are 𝑓 = 20 𝐻𝑧 and 𝐴 = 0.875 𝑚𝑚 [89]. Further, as the temperature is not measured directly in the experimental study [89], the temperature boundary conditions are therefore prone to experimental errors. Thus, considering these values to be the exact conditions in numerical study is bound to yield different results. This is attributed to the considerable variation in the thermo-physical properties even with a slight change in the distance from the critical point. In lieu of these considerations, the temperature at the boundary is adjusted to attain a good agreement between the experimental and numerical observations. Subsequently, the error between experimental and numerical conditions is evaluated and the result shown here correspond to 22% error in the initial temperature (𝑇 𝑖 ) and 12% error in temperature imposed at the boundary, 𝑖. 𝑒. wall temperature (𝑇 𝑤 ), which is a good limit of acceptance. The vibration of the cell in the 𝑥-direction induces inertial forces in the TBL as well as the bulk as marked in Fig. 4-7(b). However, owing to the difference in density, a relative motion is induced due to the difference in inertial forces between these two layers. This difference results in a Bernoulli-like pressure difference due to the velocity difference between these layers. The

Analysis of sink-zones in two-dimensional case

In order to explain the formation of sink-zones in 2D wherein the direction of vibration (or acceleration) is normal to the direction of temperature, consider schematic of the process as presented in the Fig. 5-10(a). For the sake of simplicity, only top-left corner of the cell is shown and similar explanation follows for the other regions in the cell. shows that sink-zone grows (or is larger in size) when the direction of acceleration (component of acceleration in particular) is parallel to the temperature gradient near the wall. With change in direction of vibration, the direction of acceleration changes and the conditions pertaining to a stable configuration are attained. However, even though it may not help in the growth of sink-zone but refers to the configuration when temperature on adiabatic side is increased which thereby shrinks the sink-zone by a small amount.

It is to be mentioned that the illustration made using line 𝐵1 -𝐵2 is not the only line along which the acceleration will act. Depending on the flow field, several possible combinations throughout the whole cell are possible (as illustrated in Fig. 5-10(b)). In general, the lines mentioned here will be several streamlines along which these phenomena will occur and hence lead to the formation and the spread of sink-zones. It can be concluded from the above discussion that the time of appearance of sink-zone depends primarily on the following factors, magnitude of acceleration and the compressibility for the same quench conditions. Higher value of both these factors will entail a faster onset of the sinkzone. This can be explained by the fact that a higher value of both these factors will escalate the effect of self-weight. This explains why even in case of 𝑇 𝑖 -𝑇 𝑐 = 2000 𝑚𝐾, the formation of sinkzone (Fig. 5-5(a)) was observed though very late in time. This is further illustrated by analyzing the 1D case same as in Fig. 5-6(a) with perfect gas (H2) at ambient conditions ( 𝑇 𝑖 = 278 𝐾 and 𝑃 = 101325 𝑃𝑎). Fig. 5-11 shows the results for the perfect gas and it can be observed due to very low compressibility when compared to its supercritical state, even in the presence of very high acceleration, no phenomena of attaining sink-zones is observed. This is also attributed to the fact that due to very high diffusion, the system is homogenized very rapidly. Thus, it can also be remarked that the phenomena of the fall of fluid temperature below the imposed boundary value is peculiar to near critical fluids, favored by their anomalous thermo-physical properties.

Parameters affecting onset time of sink-zones

a region wherein 𝑇 𝑡𝑟𝑎𝑛𝑠𝑖𝑡𝑖𝑜𝑛,𝐶1 -𝑇 𝑤𝑎𝑙𝑙 = 0.99 (𝑇 𝑏𝑢𝑙𝑘,𝐶1 -𝑇 𝑤𝑎𝑙𝑙 ) where 𝑇 𝑡𝑟𝑎𝑛𝑠𝑖𝑡𝑖𝑜𝑛,𝐶1 and 𝑇 𝑏𝑢𝑙𝑘,𝐶1 represents the transition and bulk temperature at 𝑥 -corresponding to 𝐶1. It can be observed here that due to uniform effect of FPE 𝑇 𝑏𝑢𝑙𝑘,𝐶1 ≈ 𝑇 𝑏𝑢𝑙𝑘,𝑡 0 -∆𝑇 and 𝑇 𝑡𝑟𝑎𝑛𝑠𝑖𝑡𝑖𝑜𝑛,𝐶1 ≈ 𝑇 𝑡𝑟𝑎𝑛𝑠𝑖𝑡𝑖𝑜𝑛 -∆𝑇, ∆𝑇 being the temperature change due to FPE which depends on the position of 𝑥 from the cell wall. This implies that though the value of transition temperature marking a change from the TBL to bulk is different at each 𝑥 , the actual thickness of the TBL is independent of ∆𝑇.

A very important observation to be highlighted here is that the TBL at any 𝑥 is not defined 𝑤. 𝑟. 𝑡 the same bulk temperature. For example, along 𝑦 direction through 𝐶1 and 𝐶2, the bulk temperature is 𝑇 𝐶0 -∆𝑇 and 𝑇 𝐶0 + ∆𝑇 respectively, 𝑇 𝐶0 being the temperature at the middle. (As was described in previous section, the FPE has negligible effect at the center, thus 𝑇 𝑡𝑟𝑎𝑛𝑠𝑖𝑡𝑖𝑜𝑛 is nearly the same as 𝑇 𝐶0 though the latter is used for the sake of brevity). Thus, the temperature difference for which the TBL is defined in both the cases will be different as the temperature of the wall (𝐴𝐵) is constant. However, on defining the TBL along the 𝑥-direction 𝑤. 𝑟. 𝑡 the bulk temperature at some other 𝑥 position (and call this as reference temperature), the thickness of the TBL will no longer be same along 𝑥 and an interesting behavior can be drawn out. Thus, the usual definition of the TBL loses its significance and it becomes necessary to introduce an alternative nomenclature in the present context, described as the relative thickness in the previous section.

In order to elucidate this phenomenon more clearly, consider the reference position at 𝐶1

where the TBL is defined 𝑤. 𝑟. 𝑡 the bulk temperature therein at time 𝑡 as 𝑇 𝑏𝑢𝑙𝑘,𝐶1 = 𝑇 𝐶0 -∆𝑇. Now, for this bulk temperature (𝑇 𝑏𝑢𝑙𝑘,𝐶1 ) or temperature difference (𝑇 𝑏𝑢𝑙𝑘,𝐶1 -𝑇 𝑤𝑎𝑙𝑙 ), the thickness of the TBL along 𝐶0 and 𝐶2 will be represented by points 𝑅 and 𝑆, as marked on the temperature profile (and also projected on 𝑦 axis) in Fig. 5-17(c) and can be represented by 𝛿 𝑇 𝐶0 -𝑇 𝑐1 and 𝛿 𝑇 𝐶2 -𝑇 𝑐1 , based on the nomenclature defined in §5.2.2. It can therefore be observed that the thickness of the TBL when defined 𝑤. 𝑟. 𝑡 different bulk temperature at 𝐶1 changes from being maximum at 𝐶1 and minimum at 𝐶2. Alternatively, considering the temperature in the middle as the reference state (which is physically more relevant as explained below), then the relative thickness of the TBL increases at 𝐶1 and decreases at 𝐶2. With the change in the direction of vibration, 𝑖. 𝑒. towards negative 𝑥-direction, the temperature profiles at 𝐶1 and 𝐶2 will be swapped as the temperature now increases at 𝐶1. The term relative thickness of the TBL becomes significant when dealing with 2D case and it is essential to examine the cell as a whole instead of considering the evolution of

CHAPTER 6. PHASE-FIELD MODELLING IN SUB-CRITICAL FLUIDS

Hitherto, the analysis of the fluid behavior has been presented in supercritical regime using single-phase model. In order to study the behavior near the critical point in sub-critical regime, a phase-field model was presented based on mass fraction as the phase-field parameter in §2.2. In this chapter, some elementary test cases are presented in order to investigate the capability of the model to capture the flow behavior of highly compressible liquid and gas phases. This is a preliminary work which can be further extended to more intricate flow dynamics and thus attention is restricted only to isothermal systems.

The two-phase system (near the critical point) is analyzed primarily for the following test cases.

• Spurious currents in stagnant vapor bubble

• Transformation of elliptical bubble to circular shape in order to minimize the energy

• Coalescence of two drops in order to minimize the energy

• Phase separation into liquid and vapor phases at constant temperature In this chapter, near critical carbon-dioxide, CO2 (𝑇 𝑐 = 304.13 𝐾), has been used with properties taken from the NIST database [START_REF]Thermophysical properties of pure fluids database[END_REF] as described in Table 6-1. Furthermore, the simulation is performed with time step 𝛿𝑡 = 10 -4 𝑠 while a constant mesh (as per the problem) has been used.

The latter is an important condition so as to ensure that the mesh conditions for interface capturing, as described at the initial state, are prevalent throughout the computational time.

Spurious currents in a stagnant bubble

In a multiphase problem, an imbalance between the discretized forces in both the phases results in small amplitude of artificial velocity known as the spurious current. The intensity of these spurious currents increases with surface tension and density ratio of the two fluids. In the context of phase-field modelling, these spurious currents arise primarily due to the gradient of the phase-In order to analyze the effect of spurious currents in our model, a quiescent CO2 vapor bubble of radius 𝑎 is suspended in its saturated liquid 100 𝑚𝐾 below the critical point in zero gravity as shown in Fig. 6-1. Symmetric boundary conditions are imposed at all the walls. The initial conditions for the phase field parameter are set using a hyperbolic tangent profile (as described in Appendix E.1 for 1D case) and is defined as,

where (𝑥 0 , 𝑦 0 ) = (5𝑚𝑚, 5𝑚𝑚) denotes the center of the bubble. In order to evaluate the effect of spurious currents, the maximum of the velocity in either direction (𝑥 or 𝑦) in the domain was noted at each time, which in ideal situation should be zero. an error of 3.7 % which is well within the acceptable limit. This illustrates that the current model is well able to capture the dynamics of highly compressible liquid and vapor phases though accuracy can be improved upon by improving the numerical parameters.

Coalescence of two liquid drops

In this case, the coalescence of two identical quiescent liquid droplets which are at equilibrium with its vapor phase is analyzed as schematically represented in Fig. 6-7 for 𝑎 = 0.8𝑚𝑚 (100 𝑚𝐾 below the critical point). At the initial condition, the drop just touch each other such that 𝑑 = 2𝑎. As described by Lamorgese and Mauri [START_REF] Lamorgese | Diffuse-interface modeling of liquid-vapor phase separation in a van der Waals fluid[END_REF], the non-equilibrium developed due to attractive force by virtue of surface tension will result in coalescence of the two drops. This is further favored by minimizing of the total energy by the interface. Interface thickness is defined by the region where 𝜙 changes its value by 90% [START_REF] Jacqmin | Calculation of Two-Phase Navier-Stokes Flows Using Phase-Field Modeling[END_REF]. Thus, we have,

E.2 Surface Tension

Surface tension is defined as the surface energy per unit area. Thus, we have 𝜎 =