Finally, my PhD experience has arrived to its end. These have been years of great challenges and intense learning. It would not have been possible to do without the support and assistance that I received from many people.

This thesis work is the result of a collaboration between the LTeN (Laboratoire de Thermique et Energie de Nantes) of the University of Nantes and Segula Technologies within the framework of a CIFRE grant from the ANRT (Association Nationale de la Recherche et de la Technologie).

First of all, I am deeply grateful for the continuous support, insight and patience of my thesis director, Pr Lingai LUO, without her constant trust and, sometimes, gentle prodding, this thesis would not have been completed.

I would like to express my sincere gratitude to my thesis co-director Dr Yilin FAN, for his daily encouragement and enthusiasm. Through these years he has always been there in several disheartened moments, ready to discuss, giving valuable suggestions and guidance.

I am deeply grateful to Dr Mustapha BENAOUICHA for introducing me to the research world through my internship and later for giving me this opportunity to join the research team and work on such an interesting topic and for his insightful comments and suggestions.

I would like to extend my sincere thanks Pr Pierre SAGAUT and Dr Jean CASTAING-LASVIGNOTTES who reviewed this manuscript and participated as jury members during my thesis defense. Further, my sincere gratitude to Pr Roland SOLIMANDO who was the president of the jury, Pr. Vincent GOETZ, and Dr. Patrick BOT who were also the members of my thesis jury, and with whom the exchanges were very interesting during the defense.

I would like to offer my special thanks to Dr Thibault NEU for his many contributions, especially the experimental results. He has certainly helped me understand some key aspects of the Liquid Piston. Discussing work with him was always exciting as well as inductive.

I would like to thank the members of LTeN and my colleagues at Segula Bouguenais for their daily assistance.

Many thanks to all the former and current doctoral and post-doctoral students at LTeN. It was very enriching to discuss science, society or politics. It was also very nice to co-organize the PhD students' breakfast with Arthur and Rawad.

In addition, I could not have completed this dissertation without the support of my friends, Nabil, Hamza, Zakaria and Khalil who provided stimulating discussions as well as happy distractions to rest my mind outside of my research.

Finally, I would like to thank my parents, my brother and my sister for their wise counsel and sympathetic ear. You are always there for me.

Chapter 1 Introduction

All over the world, the electrical network is undergoing many changes linked to the awareness of the environmental impact of fossil energy sources and their rareness and to the exponential growth of the need for electrical energy. The prospect of massive penetration of renewable energies of all kinds in a stable and reliable electrical network is strongly compromised by the intermittency of production. Thus, storage solutions are proving to be a definite ally in supporting renewable energy and the stability of the electrical grid. It even appears to be a quasi-systematic and necessary condition for the integration of high-powered renewable energy production facilities. Moreover, with the development of renewable energies, decentralization and independence of the electrical network on the scale of an island, a region or even an entire country are becoming more and more frequent. The integration of storage systems is presented in this case as an essential means of managing the stability and flexibility of the grid. For the deployment of decarbonated production means, the ocean offers a strong potential to be exploited due to its many renewable energy resources available (wind, marine currents, waves, etc.).

In this context, Segula Technologies has been developing a solution for storing electrical energy in the form of compressed air in the ocean since 2013. This storage system, named REMORA (Réservoir d'Énergie en Milieu Océanique par Rétention d'Air), aims to absorb excess electricity generated in order to restore it during periods of low production or high demand (Figures 1.1, 1.2). One of they key components of the REMORA project is the Liquid Piston. It has proven its efficiency over solid piston especially for under-water applications. This technology, which is not a new concept lacks proper understanding of the physics related to the flow and heat transfer generated by the compression/expansion processes using it.

The flow and heat transfer associated to the compression/expansion using liquid piston are complex. It is a compressible, non-isothermal multiphase one.

This thesis is a collaboration between Segula Technologies and the LTEN (Labotratoire de Thermique et Energie de Nantes) of Nantes University and aims to model numerically the flow and heat transfer inside a liquid during the compression/expansion processes. The obtained Figure 1.1: Three main steps of energy storage stage in REMORA using a liquid piston [START_REF] Neu | Etude expérimentale et modélisation de la compression quasi isotherme d'air pour le stockage d'énergie en mer[END_REF] Figure 1.2: Three main steps of energy production stage in REMORA using a liquid piston [START_REF] Neu | Etude expérimentale et modélisation de la compression quasi isotherme d'air pour le stockage d'énergie en mer[END_REF] CHAPTER 1. INTRODUCTION results will therefore enable us to better describe the flow evolution and help understand the physics linked to it. The LBM method (Lattice-Boltzmann Method) was initially chosen to conduct this study because of its linear character which makes it efficient for large calculations. However, the method presented several scientific obstacles which proved to be very difficult to solve. It was therefore chosen to pursue this study with a more classical approach : the resolution of the Navier-Stokes equations with the Finite Volume Method coupled to Volume of Fluid Method for interface tracking.

Objectives of the thesis

The main result expected from this research work is to carry out a global study of the compression/expansion processes by liquid piston. This consists in carrying out numerical modeling in order to understand the evolution of the flow, velocity, temperature and pressure fields. Methods to measure the efficiency of a liquid piston over a complete cycle (compression-expansion) are expected in order to better choose the numerical model to use.

In addition, the physical exploitation of the results will allow a better understanding, a qualification and a quantification of the physical characteristics associated to the flow and heat transfer in the liquid piston. A parametric study on the effects of the different parameters of the liquid piston and the compression/expansion will enable to identify improvement points, as to prepare, in perspective of this work, for a topological optimization of the liquid piston air compression/expansion of the REMORA system.

Outline of the manuscript

The thesis manuscript is organized as follows. A literature review study that introduces the liquid piston, the physics linked to it, the different methods to improve its efficiency through Heat Transfer Enhancement and finally the different challenges towards an application for CAES systems are presented in Chapter 2.

The different attempts realized using the LBM to model the liquid piston compression flow are described in the chapter 3. It concerns the modifications and implementation of the Shan-Chen model, free energy and Volume of Fluid method to enable simulations of multiphase flows with high density ratios using the LBM.

The numerical model based on the Navier-Stokes equations which was implemented and applied to model compression and expansion processes is presented in chapter 4. The obtained results have been compared to experimental ones to validate the model. Further analysis of the flow and heat transfer has been carried out and key flow parameters are presented. We analyzed the flow structure settling, its evolution and the origin of its disruption. We also showed that the flow changes from structured to chaotic after the transition occurs. To enable a better understanding of each compression parameter, a parametric study for piston speed and ambient temperature has been conducted. This allowed to show these effects on the process efficiency.

As for the compression process, the expansion is studied through an analysis of the evolution of flow and heat transfer during the process in chapter 5. We showed, as for the compression, that the flow undergoes several stages and analyzed each one. Later, the compression-isobaric cooling-expansion cycle results are presented in Chapter 5. A parametric study for piston speed and ambient temperature has also been conducted and showed the effects of each parameter on expansion efficiency.

Two new numerical models to couple the heat transfer in the solid walls to the flow (Conjugate Heat Transfer) are developed and implemented in OpenFOAM, as described in Chapter 6. A comparison of each model has been performed to validate both. Effects of solid walls on the compression and expansion are highlighted through parametric study of three different solid materials.

Finally, main conclusions and perspectives are presented.

Chapter 2

Review on Liquid Piston technology for compressed air energy storage

Introduction

The energy consumption worldwide has increased by 21% from year 2009 to 2019 and is expected to grow with more than 50% by 2050 [START_REF]EIA projects nearly 50% increase in world energy usage by 2050, led by growth[END_REF]. To meet this demand, the world energy production reached 14 421 Mtoe (million tonnes of oil equivalent) in 2018, with more than 81% driven by fossil fuels (natural gas, coal and oil) [START_REF]Global Energy Review[END_REF]. In the meantime, awareness has been raised on the concomitant environmental impacts due to the greenhouse gas emission, evidenced by the signing of the Paris Agreement [165] in 2016 aiming at substantially reducing the global warming and the risks of climate change. In this context, renewable energies, with their increasing contribution year by year to the global energy market, are expected to play a major role to achieve the COP21 objectives for a more sustainable and decarbonized future [165].

Most renewable energy sources (e.g., solar, wind, tidal) show the intermittent character, bringing more difficulties for their exploitation at large scale. In this regard, the Energy Storage Systems (ESS) have become an essential element for the power generation plants driven by renewable sources so as to augment their reliability and dispatchability [START_REF] Zhou | A review of thermal energy storage in compressed air energy storage system[END_REF]. Moreover, the integration of ESS brings additional economic benefits by storing the excess power produced at times of low demand or low generation cost to be used during the peak periods with a higher price, accelerating thereby the pathway to a flexible, low-cost and electrified future [START_REF] Luo | Overview of current development in electrical energy storage technologies and the application potential in power system operation[END_REF][START_REF] Zhou | A review of thermal energy storage in compressed air energy storage system[END_REF]. Different types of ESS have been developed and utilized in the past decades, classified as follows based on their storage mechanisms.

• Electrical energy storage: electrostatic energy storage including capacitors and supercapacitors, magnetic/current energy storage including Super-conducting Magnetic Energy Storage (SMES), etc. [START_REF] Luo | Overview of current development in electrical energy storage technologies and the application potential in power system operation[END_REF];

CHAPTER 2. REVIEW ON LIQUID PISTON TECHNOLOGY FOR COMPRESSED AIR ENERGY STORAGE

• Thermal energy storage: sensible, latent and thermochemical heat storage [START_REF] Pelay | Thermal energy storage systems for concentrated solar power plants[END_REF][START_REF] Pelay | Integration of a thermochemical energy storage system in a Rankine cycle driven by concentrating solar power: Energy and exergy analyses[END_REF][START_REF] Pelay | Dynamic modeling and simulation of a concentrating solar power plant integrated with a thermochemical energy storage system[END_REF][START_REF] Edem | A review on long-term sorption solar energy storage[END_REF];

• Mechanical energy storage: kinetic energy storage (flywheels), potential energy storage Pumped Storage Hydrodynamics (PSH) and Compressed Air Energy Storage (CAES), etc. [START_REF] Guewouo | Identification of Optimal Parameters for a Small-Scale Compressed-Air Energy Storage System Using Real Coded Genetic Algorithm[END_REF];

• Chemical energy storage: hydrogen, methane, etc. [START_REF] Luo | Overview of current development in electrical energy storage technologies and the application potential in power system operation[END_REF];

• Electrochemical energy storage: batteries (lead-acid, lithium-ion, nickel-cadmium, etc.) and fuel cells [START_REF] Venkataramani | A review on compressed air energy storage -A pathway for smart grid and polygeneration[END_REF].

Among all these energy storage technologies, the PSH and CAES have been proven to be the most adapted one to store electricity in large scale [START_REF] Budt | A review on compressed air energy storage: Basic principles, past milestones and recent developments[END_REF]. PSH is a mature technology featured by its large power capacity (100-3000 MW), long storage period (1-24 h) and high RTE (71% to 85%) [START_REF] Chen | Progress in electrical energy storage system: A critical review[END_REF]. Therefore, it has dominated the energy storage market accounting for 95% of the global capacity (over 100 GW) [START_REF] Letcher | Storing electrical energy[END_REF]. Its main drawback lies in the rareness of available sites for two large reservoirs and one or two dams, a long construction time required (typically about 10 years) and high investment costs (typically hundreds to thousands of million US dollars). Moreover, removing trees and vegetation from the vast land prior to the flooding of reservoirs also appends additional environmental impacts. For these reasons, the planning of new PSH sites has been limited in the past years [START_REF] Chen | Progress in electrical energy storage system: A critical review[END_REF] despite that the overall new power installation is still dominating [START_REF]Will pumped storage hydropower expand more quickly than stationary battery storage? Analysis[END_REF]. In the meantime, the CAES technology has experienced a rapid growth in recent years, demonstrating the potential to become the next major ESS for power plants [START_REF] Budt | A review on compressed air energy storage: Basic principles, past milestones and recent developments[END_REF].

Compressed Air Energy Storage

Storing electrical energy using a compressed air system can be dated back to the early 20 th century [START_REF] Chen | Progress in electrical energy storage system: A critical review[END_REF], in a way that it decoupled a conventional gas turbine cycle into separated compression and expansion processes [START_REF] Haisheng | Compressed air energy storage[END_REF]. The energy was stored in the form of elastic potential energy of compressed air and could be restituted by expanding it through turbines to reproduce the electricity. The waste heat of the exhaust gas can be captured through a recuperator before being released to the atmosphere. In modern CAES systems, the compressed air can be stored either in man-made containers at the ground level or underground (salt caverns, hard rock caverns, saline aquifers) [START_REF] Li | Compressed air energy storage: Characteristics, basic principles, and geological considerations[END_REF][START_REF] Olabi | Compressed air energy storage systems: Components and operating parameters -A review[END_REF][START_REF] Matos | Overview of Large-Scale Underground Energy Storage Technologies for Integration of Renewable Energies and Criteria for Reservoir Identification[END_REF]. Offshore and underwater storage systems have also been tested and are under rapid development [START_REF] Budt | A review on compressed air energy storage: Basic principles, past milestones and recent developments[END_REF]. Capital costs for CAES facilities depend on the storage conditions, ranging typically between 400$ and 800$ per kW [START_REF] He | Optimal selection of air expansion machine in Compressed Air Energy Storage: A review[END_REF]. Above-ground systems usually request higher investments than underground ones [START_REF] Chen | Progress in electrical energy storage system: A critical review[END_REF][START_REF]Will pumped storage hydropower expand more quickly than stationary battery storage? Analysis[END_REF].

Since the gas compression and expansion processes involve a large amount of heat, the efficiency and storage capacity of the CAES systems depend strongly on the heat transfer CHAPTER 2. REVIEW ON LIQUID PISTON TECHNOLOGY FOR COMPRESSED AIR ENERGY STORAGE behaviors. Different variants of CAES exist, distinguishing themselves by the manner of thermal management.

• Diabatic-CAES (D-CAES): this technology pressurises and heats air by combusting a fuel (usually natural gas) and expands it though turbo turbines to generate electricity. It is the basic but the only industrialized version of the CAES at large scale [START_REF] Raju | Modeling and simulation of compressed air storage in caverns: A case study of the Huntorf plant[END_REF][START_REF] Luo | Overview of current development in electrical energy storage technologies and the application potential in power system operation[END_REF][START_REF] Budt | A review on compressed air energy storage: Basic principles, past milestones and recent developments[END_REF]. Nevertheless, gas emissions are inevitable due to the fuel combustion required to drive the gas compressor.

• Adiabatic-CAES (A-CAES): the main improvement for this second generation of CAES is to retain the heat of the compressed air rather than wasting it, to be used then in the expansion process. As a result, it renders higher efficiency than that of D-CAES systems.

• Isothermal-CAES (I-CAES) : this emerging technological variant tries to prevent the temperature increase in the compressors (during charging) and the temperature drop (during discharging) in the expansion devices, permitting isothermal cyclic operations.

The main advantage of I-CAES is that fossil-based backup is not needed to start the expansion, implying higher efficiencies and more environmental friendliness [START_REF] Budt | A review on compressed air energy storage: Basic principles, past milestones and recent developments[END_REF][START_REF] Luo | Overview of current development in electrical energy storage technologies and the application potential in power system operation[END_REF].

Liquid Piston technology

It is by far very difficult to achieve isothermal compression and expansion using conventional devices. The heat transfer of the system has to be very efficiently mastered to dissipate (almost) all the heat generated in the compression process, and vice versa, to provide necessary heat for the expansion. Therefore, all I-CAES concepts known so far are based on piston machinery capable of performing a comparably slow compression or expansion, leaving enough time for heat exchange to approach a near-isothermal condition. In this context, the interest over LP for CAES usage has been rapidly growing in recent years. In the LP compressor, a column of liquid (usually water) is pumped in to compress a fixed amount of gas (usually air). The main advantages of LPs compared to solid piston machinery is that the gas leakage can be avoided and the dissipations due to friction be largely reduced, leading to a higher efficiency. Many CAES projects with LP concept have been developed or are under construction, as listed in Table 2.1. For example, the GLIDES project (at ground level) by Oak Ridge National Laboratory using liquid spray injection to enhance heat transfer reported a 66% RTE for their LP prototype [START_REF] Odukomaiya | Thermal analysis of near-isothermal compressed gas energy storage system[END_REF]. More recently, other modifications to the system have been introduced and several tests have been carried out by replacing the noncondensable air with a condensable gas (CO 2 , synthetic refrigerants, hydrocarbon refrigerants, etc.) [START_REF] Odukomaiya | Near-isothermal-isobaric compressed gas energy storage[END_REF]. Near isothermal and near isobaric charging/discharging processes could be achieved. CHAPTER 2. REVIEW ON LIQUID PISTON TECHNOLOGY FOR COMPRESSED AIR ENERGY STORAGE Moreover, the LP concept suits especially the underwater/sub-sea CAES systems in that the heat exchange can be largely enhanced due to the cold water environment with high heat capacity. Park et al. [START_REF] Park | Analysis and Proof-of-Concept Experiment of Liquid-Piston Compression for Ocean Compressed Air Energy Storage (Ocaes) System[END_REF] have developed a simplified model (steady-state operation; ideal gas) and a proof-of-concept experiment for a compression process by LP. Their study was later completed by Patil & Ro [START_REF] Vikram | Modeling of liquid-piston based design for isothermal ocean compressed air energy storage system[END_REF][START_REF] Vikram | Design of Ocean Compressed Air Energy Storage System[END_REF]. It has been estimated that the RTE of their 2 MWh stored energy CAES system could reach 51% with a compressor/expander efficiency of 79%, and up to 73% with a higher compressor/expander efficiency of 95%. More recently, the REMORA project [START_REF] Neu | Experimental study of convective heat transfer during liquid piston compressions applied to near isothermal underwater compressed-air energy storage[END_REF][START_REF] Maisonnave | Optimal energy management of an underwater compressed air energy storage station using pumping systems[END_REF][START_REF]REMORA | Segula Technologies[END_REF] showed that near-isothermal condition could be achieved with a LP-based compression/expansion cycle due to the increased air storage pressure and the reduced heat loss. A 95% compression efficiency could be achieved by the LP, leading to a 70% RTE for the CAES system [START_REF] Neu | Experimental study of convective heat transfer during liquid piston compressions applied to near isothermal underwater compressed-air energy storage[END_REF][START_REF]REMORA | Segula Technologies[END_REF]. Table 2.1 lists different CAES projects, where it could be seen that D-CAES is the only implemented technology already in service. I-CAES, although still underdevelopment, could be promising because it can theoretically offer higher RTE. CHAPTER

Objectives of this chapter

In the open scientific literature, the LP concept has also received increasing research attention, indicated by the significantly increasing number of yearly documents since 2008 (Fig. 2.1). A great number of researches have been devoted to study the heat transfer behaviors through monitoring the evolution of gas temperature [START_REF] Van De Ven | Liquid piston gas compression[END_REF][START_REF] Vikram | Experimental investigation of heat transfer in liquid piston compressor[END_REF][START_REF] Neu | Experimental study of convective heat transfer during liquid piston compressions applied to near isothermal underwater compressed-air energy storage[END_REF] and by implementing various techniques to enhance the heat transfer [START_REF] Vikram | Design of Ocean Compressed Air Energy Storage System[END_REF][START_REF] Odukomaiya | Experimental and analytical evaluation of a hydropneumatic compressed-air Ground-Level Integrated Diverse Energy Storage (GLIDES) system[END_REF]. Few studies have focused on the flow patterns and their transition, showing that the fluid flow hydrodynamics are rather complex and difficult to predict especially when coupled with heat transfer [START_REF] Van De Ven | Liquid piston gas compression[END_REF][START_REF] Neu | Experimental investigation of internal air flow during slow piston compression into isothermal compressed air energy storage[END_REF]. The effects of some design and operational parameters have been studied but their relationships and coherence to the thermofluidic behaviors are still veiled. Moreover, the systematic state-ofart literature review to keep track with all the relevant studies and the new developments is lacking. This review chapter is thereby expected to be the first to gather and survey previous works on LPs for CAES, a research topic with rapid advances in recent years. In particular, a comprehensive and combined knowledge on the thermodynamics, the fluid flow and heat transfer mechanisms within the LP is provided. Special focus is given on the thermal management issue, and measures proposed and implemented by different researchers to enhance the heat transfer, aiming at achieving a highly-efficient LP for I-CAES. The contributions of this chapter are important because it provides a technology state-of-the-art and a progress overview of the LP technology covering from the fundamental transport phenomena to the R&D projects. It may enable interested researchers to timely grasp the latest advancements and existing challenges for the investigation and application of LP technology in future (underwater) CAES systems.

The rest of the chapter is organized as follows. The operating principles and the thermodynamic background of the LP are described in section 2. The fluid flow patterns and the heat transfer characteristics inside the LP during the compression are presented as well. Various techniques to enhance the heat transfer are summarized and compared in section 3. Finally, conclusions and perspectives are presented in section 4.

Liquid piston for energy storage

LP is in fact not a new concept but can be dated back to the Humphrey pumps in 1906 [START_REF] Humphrey | An Internal-Combustion Pump, and other Applications of a New Principle[END_REF], which is a large internal combustion gas-fuelled LP pump used for large-scale water supply projects. Later on, LPs were also used for Stirling engines and Stirling pumps. The main applications of these configurations were the usage of a LP as a cooler to serve superconducting computers or as a heat pump [START_REF] Gerstmann | Isothermalization of Stirling Heat-Actuated Heat Pumps Using Liquid Pistons[END_REF][START_REF] West | Liquid piston stirling engines[END_REF]. In 2009, Van de Ven & Li first proposed the LP concept for energy storage through an open accumulator [START_REF] Van De Ven | Liquid piston gas compression[END_REF].

Operating principle

The LP concept, compared to traditional solid piston, is featured by the use of a compression chamber where a liquid (often water) is used to increase or decrease the pressure of a gas (often air or hydrogen). The gas and liquid phases are naturally separated by their density difference, as long as the piston advance velocity has no particular effects on the gas-liquid interface. But when a portion of gas becomes trained in the liquid due to the high pressure, the decreased bulk modulus of the liquid could cause the cavitation problem in low pressure areas of the hydraulic system. Van de Ven & Li [START_REF] Van De Ven | Liquid piston gas compression[END_REF] proposed to use a blade diaphragm in the liquid column to separate the LP fluid from the working fluid of the hydraulic system, allowing the pressure transfer across the diaphragm with the minimal loss. Other technological alternatives also exist, i.e., by coupling a LP with a solid piston to handle the trained gas problem [START_REF] Qin | Liquid piston compression efficiency with droplet heat transfer[END_REF].

The general operation principle of LP, as has already been studied [START_REF] Van De Ven | Liquid piston gas compression[END_REF] and proved to be feasible for CAES projects [START_REF] Neu | Experimental study of convective heat transfer during liquid piston compressions applied to near isothermal underwater compressed-air energy storage[END_REF][START_REF] Vikram | End-to-end efficiency of liquid piston based ocean compressed air energy storage[END_REF][START_REF] Park | Analysis and Proof-of-Concept Experiment of Liquid-Piston Compression for Ocean Compressed Air Energy Storage (Ocaes) System[END_REF], is briefly described as follows:

• During the charging stage, the surplus electricity drives the hydraulic pump to fill the liquid into the closed compression chamber from the bottom (Fig. 2.2a). Since the AIR ENERGY STORAGE chamber is hermetically sealed, the trapped gas in the top region undergoes its pressure increase from P 0 up to the setting point P c until it is discharged to the gas reservoir by opening the outlet valve. Note that the ratio between P c and P 0 is called the CR.

• During the discharging stage, the compression chamber initially full of liquid is filled with compressed gas released from the gas reservoir. The expansion of the compressed gas in the chamber pushes the liquid to flow back through the turbine to generate electricity (Fig. 2.2b). Once a certain volume of gas is present in the chamber, the inlet valve is closed.

Valves are used accordingly to control and switch the process from compression mode to expansion mode by changing the inlet/outlet domains [START_REF] Neu | Experimental study of convective heat transfer during liquid piston compressions applied to near isothermal underwater compressed-air energy storage[END_REF][START_REF] Van De Ven | Liquid piston gas compression[END_REF]. Besides air or hydrogen, helium was also used as the working gas [START_REF] Schober | CFD-simulation and experimental validation of heat transfer in liquid piston compressors[END_REF]. Piya et al. [START_REF]Heat Transfer, Fluid Flows, and Thermal Systems, Parts A[END_REF] made a 2D finite difference modeling to test different combinations of liquid (DTE 25 hydraulic fluid (oil) or water) and gas (helium or air). Their results showed that helium/water combination had the best compression efficiency η c , going up to 61.9% compared to 54.31% for air/water and 23.58% for air/DTE25 based on the same configuration.

Thermodynamic background

Lemofouet & Rufer [START_REF] Lemofouet | Investigation and optimisation of hybrid electricity storage systems based on compressed air and supercapacitors[END_REF] showed that the usage of LP-based CAES can involve different cycles, including the isothermal cycle (A-D path as shown in Fig. 2.3, the same below), the Joule cycle (A-B-D-E-A: a combination of the adiabatic and isobaric processes), the Otto cycle (A-B-C-D-E-F-A: a combination of adiabatic and isochoric processes) and the polytropic cycle (A-B'-C'-D-F' and A-B'-D-E'-A). The gray area enclosed in each cycle's path represents the energy losses for that cycle. It can be seen that the isothermal cycle without energy loss AIR ENERGY STORAGE is the optimal one among these cycles. The closer the path to the isothermal one the less losses the system generates (Fig. 2.3). In both closed and open systems, the Joule cycle is more efficient than the Otto cycle for small values of CR [START_REF] Lemofouet | A hybrid energy storage system based on compressed air and supercapacitors with maximum efficiency point tracking (MEPT)[END_REF][START_REF] Lemofouet | Investigation and optimisation of hybrid electricity storage systems based on compressed air and supercapacitors[END_REF][START_REF] Dib | Thermodynamic investigation of quasi-isothermal air compression/expansion for energy storage[END_REF]. Dib et al. [START_REF] Lemofouet | A hybrid energy storage system based on compressed air and supercapacitors with maximum efficiency point tracking (MEPT)[END_REF][START_REF] Dib | Thermodynamic investigation of quasi-isothermal air compression/expansion for energy storage[END_REF] also studied the thermodynamics of different compression/expansion cycles with a LP and showed that a Joule cycle is more adapted for a LP gas compression for energy storage. In an open-accumulator CAES system air at initial state (T 0 , P 0 , V 0 ) is compressed by a LP through a P -V trajectory (ζ) to a second state (T c , P c , V c ) with a CR = Pc P 0 and ER = Pc P f . Then, the compressed air is ejected into the open accumulator and cooled isobarically to a smaller volume (V iso ), by adding liquid into the bottom of the accumulator until the air reaches the final stored state (T 0 , P c , V iso ). For an internally reversible steady flow device, the specific work could be calculated by:

W = Vc V 0 P dV (2.1)
The total amount of work input to store the air at the final state is then:

W c actual = Vc V 0 (P (t) -P 0 )dV Compression f low + 1 - 1 CR (P c V c -P 0 V 0 ) Isobaric cooling + (P c -P 0 )V iso P umping (2.2)
After the storage stage, the total potential energy E p that could be extracted from the stored compressed air is shown as the area shaded by the solid vertical lines (Fig. 2.4a).

E p = (P c -P 0 )V iso Air discharge to isobaric tank storage The compression efficiency is then defined as Eq. 2.4 [START_REF] Yan | Experimental study of heat transfer enhancement in a liquid piston compressor/expander using porous media inserts[END_REF]:

+ P 0 V 0 ln(CR) + 1 CR -1 Isothermal expansion (2.3) (a) (b)
η c = E p W c actual = (P c -P 0 )V iso + P 0 V 0 ln(CR) + 1 CR -1 Vc V 0 (P (t) -P 0 )dV + 1 -1 CR (P c V c -P 0 V 0 ) + (P c -P 0 )V iso (2.4)
When the actual compression trajectory approaches the isothermal one, the area shaded by vertical lines gets closer in magnitude to that by the horizontal lines, and η c increases. On the contrary, when the trajectory gets closer to the adiabatic one, the difference between those two areas becomes larger and η c decreases.

For the expansion stage, the stored compressed air (T 0 , P c , V 0 ) is firstly ejected from the open-accumulator isobarically and then expanded in the LP at (T e , P e , V e ). The process also follows a trajectory between the isothermal one and the adiabatic one (Fig. 2.4b). Similarly, the expansion efficiency η e can be calculated in Eq. 2.5 [START_REF] Yan | Experimental study of heat transfer enhancement in a liquid piston compressor/expander using porous media inserts[END_REF]:

η e = W e actual E p = Ve V 0 (P (t) -P f )dV + (P c -P f )V 0 P c V 0 (ln(ER)) -(P f (V iso -V 0 )) + (P c -P f )V 0 (2.5)
It can be seen from Eq. (2.4) that the compression efficiency η c depends on the compression time as well as the CR. It also depends on the modeled cycle. Since many cycles have been studied, different formulas to calculate the efficiencies exist. While efficiency can be calculated in comparison to an isentropic one. The formulas presented fo LP are generally isothermal efficiencies. Patil & Ro [START_REF] Vikram | Energy and Exergy Analysis of Ocean Compressed Air Energy Storage Concepts[END_REF]

η c = ln(CR) + 1 CR -1 CR n-1 n -1 n-1 + CR -1 n -1 + (CR -1) CR -1 n -1 CR (2.6
)

η e = 1 - ( 1 ER ) n-1 n n-1 -1 ER n-1 n + 1 ER ln(ER) + 1 ER -1 (2.7)
Mutlu & Kiliç [START_REF] Mutlu | Effects of piston speed, compression ratio and cylinder geometry on system performance of a liquid piston[END_REF] have numerically examined the effect of piston speed, the CR and the cylinder geometry on η c and required power of a LP. It was found that the η c decreased from 87.93% to 68.89% when CR increased from 5 to 80. The piston velocity had only a small impact (0.5%) on η c once a certain velocity (U pist = 0.05m/s) had been reached. Decelerating the piston velocity at higher pressures did not affect the total input work nor the efficiency of the process. Their study recommended to use long cylinders with small diameters to achieve a higher η c , e.g., up to 91.42% with the proper length/diameter ratio L/D [START_REF] Mutlu | Effects of piston speed, compression ratio and cylinder geometry on system performance of a liquid piston[END_REF]. For low pressure compression, an important flow rate (velocity) is usually imposed to compensate the power level and the system performance. Both the convective heat transfer coefficient and the frictional pressure drop will increase with the increasing piston advance velocity, leading to higher energy density at the cost of higher hydraulic loss. The same thing applies to the L/D ratio and a trade-off needs thereby to be made.

Flow patterns during air compression

The fluid flow regime inside a LP column during dynamic operation is rather complex, involving unsteady compressible multiphase flow (liquid/gas) with conjugated heat transfer. The Reynolds number Re is usually used to characterize the flow pattern within a cylindrical column (uniform cross-sectional area):

Re L = ρ L U D µ L = U D ν L (2.8)
Or :

Re G = ρ G U D µ G = U D ν G (2.9)
Where ρ is the density of the fluid, U is the advance velocity (of the interface) in the piston column, D is the (hydraulic) diameter of the column, µ and ν the dynamic and kinematic viscosity, respectively.

However, one should distinguish the Re for the liquid phase (Re L ) and that for the gas phase (Re G ). Since a constant liquid volume flow rate is often assumed by the liquid pump, the liquid velocity as well as the Re L are considered thereby as constant during the compression process. On the contrary, the gas flow is much more complex than a simple pipe flow because the gas does not have the same velocity as water all over the piston column, as CHAPTER 2. REVIEW ON LIQUID PISTON TECHNOLOGY FOR COMPRESSED AIR ENERGY STORAGE reported by Mutlu & Kiliç [START_REF] Mutlu | Effects of piston speed, compression ratio and cylinder geometry on system performance of a liquid piston[END_REF] and Schober et al. [START_REF] Schober | CFD-simulation and experimental validation of heat transfer in liquid piston compressors[END_REF]. It has been observed by Neu et al. [START_REF] Neu | Etude expérimentale et modélisation de la compression quasi isotherme d'air pour le stockage d'énergie en mer[END_REF][START_REF] Neu | Experimental study of convective heat transfer during liquid piston compressions applied to near isothermal underwater compressed-air energy storage[END_REF] that during compression many re-circulation zones could appear next to the walls (upper, bottom and lateral) of the piston column. The numerical results of Gouda et al. [START_REF] Mehdi | Méthode VOF pour la simulation numérique de l'écoulement de l'air comprimé et du transfert thermique associé dans un piston liquide[END_REF] also shown that the maximum local velocity of air could reach 8 times the velocity of water. The gas phase actually undergoes regime transition from laminar to turbulent along with the advancement of the piston (air and water are the working fluids), as illustrated in Fig. 2.5. At first, an ascending central air flow becomes established during the short start-up period of compression (Fig. 2.5 (a)). This central rising flow having a velocity magnitude greater than the water velocity will then disturb the stationary air at the upper part of the column, initiating the air movement. In a little while, the central rising air flow occupies the entire compression chamber while the falling flow occurs close to the walls of the column (Fig. 2.5 (b)). Recirculation zones are thereby present at both ends of the compression chamber, i.e. at the column head where the central rising air is deviated to the outer periphery and at the air-water interface where the descending air rejoins toward the center of the column. At the boundary between the ascending and descending air flows, small instabilities would cause small amplitude of movements between the two velocity regions. When the central air velocity is significantly higher than the water velocity, the Rayleigh-Taylor instabilities [START_REF] Langdon-Arms | Rayleigh-Taylor instability in oscillating liquid pistons[END_REF] become stronger and will trigger swirling flows and vortices (Fig. 2.5 (c)). This transition regime is rather short and in a brutal way, the global structure of the air flow disintegrates CHAPTER 2. REVIEW ON LIQUID PISTON TECHNOLOGY FOR COMPRESSED AIR ENERGY STORAGE into a multitude of small vortices as shown in Fig. 2.5 (d), indicating the turbulent feature of the flow. Such transition of flow regime has also been confirmed by Gouda et al. [START_REF] Mehdi | Méthode VOF pour la simulation numérique de l'écoulement de l'air comprimé et du transfert thermique associé dans un piston liquide[END_REF] through CFD (Computational Fluid Dynamics) simulation of the compression process with the VOF (Volume-Of-Fluid) model.

Many other formulations of Re number have therefore been proposed to better characterize the gas flow regimes and their transition. For example, Adair et al. [START_REF] Adair | Instantaneous Heat Transfer To the Cylinder Wall in Reciprocating Compressors[END_REF] and Liu & Zhou [START_REF] Liu | Heat Transfer Between Gas and Cylinder Wall of Refrigerating Reciprocating Compressor[END_REF] used the transient equivalent diameter and the swirl velocity (ω g ) instead of the averaged piston speed:

Re = ρ(D/2) 2 ω g µ (2.10)
Van de Ven & Li [START_REF] Van De Ven | Liquid piston gas compression[END_REF] used the fully-developed pipe flow model to calculate the Re (Eq. 2.11):

Re = U ave D ν (2.11)
where U ave is the average gas velocity, estimated as half of the piston (liquid) velocity and D is the characteristic length (diameter of the piston). Neu et al. [START_REF] Neu | Experimental study of convective heat transfer during liquid piston compressions applied to near isothermal underwater compressed-air energy storage[END_REF][START_REF] Neu | Experimental investigation of internal air flow during slow piston compression into isothermal compressed air energy storage[END_REF] observed that the transition from laminar to turbulent flow seems to happen at different stroke times depending on the operational conditions (D, L, U pist , etc.). Their experimental measurements of air temperature in the compression column also showed a switch from a slow evolution to a fast oscillation during the compression. An empirical correlation has been proposed to predict the position of the piston front at the moment of regime transition (Eq. 2.12):

L * tr = -0.0344L 0 + 109U pist D 2 + 0.0227 D P 0 P ref -0.645 * √ U pist (2.12)
Where L 0 is the initial piston length, L * tr is the non-dimensional piston position (water/air interface) when the transition occurs, U pist is the piston inlet velocity, P ref is the reference pressure, generally equal to the atmospheric pressure and P 0 the initial pressure. This correlation has been experimentally verified within the range of

P 0 = 101325 P a, 2 m ≤ L 0 ≤ 6 m, 0.08 m/s ≤ U pist ≤ 1.25 m/s and 3 × 10 -2 m ≤ D ≤ 10 -1 m.

Heat transfer analysis

The heat transfer characteristics inside a LP column is also complex, due to various heat transfer modes involved. The increasing gas (air) temperature during the compression will cause the heat transfer between the gas and liquid which are in direct contact. Moreover, the convection heat transfer, which is the most dominant, occurs between the piston wall and the gas or/and liquid. Depending on the advancement of the gas-liquid interface, the flow regimes, the thermal resistance of wall and the temperature of the external environment, the liquid may absorb heat from the wall, or vice-versa, release the heat towards the wall. AIR ENERGY STORAGE Nikanjam & Greif [START_REF] Nikanjam | Heat transfer during piston compression[END_REF] measured the wall temperature variation of a LP to estimate the heat flux. The liquid chosen was oil, and two gases (argon and air) were compressed with a thin film gauge placed at the end wall of the compression chamber. Only conduction heat transfer was supposed to happen at the wall. Kermani & Rokni [START_REF] Arjomand | Heat transfer analysis of liquid piston compressor for hydrogen applications[END_REF] developed a detailed heat transfer model considering both the heat transfer through the gas-liquid interface and through the wall. Some critical parameters have been identified to maximize the amount of heat that can be extracted from the compressed gas for a hydrogen compression technology. A detailed modeling of heat transfer inside the LP has been performed by Odukomaiya et al. [START_REF] Odukomaiya | Thermal analysis of near-isothermal compressed gas energy storage system[END_REF], as shown in Fig. 2.6. It has been supposed that the upper wall temperature in contact with gas has a uniform temperature (T W , G ) and the one in contact with liquid (T W , L ) is also uniform. These temperatures evolve in time as the piston front moves upward (compression) or downward (expansion). The internal energy variation of the gas or liquid phase can be calculated by Eq. 2.13 or Eq. 2.14, respectively. Note that Redlich-Kwong equation of state was used to model the gas inside the column.

CHAPTER 2. REVIEW ON LIQUID PISTON TECHNOLOGY FOR COMPRESSED AIR ENERGY STORAGE m G c v,G dT G dt change of the gas energy = -h G,L A G,L (T G -T L )
m L c p,L dT L dt change of the liquid energy = h G,L A G,L (T G -T L )
heat transfer with the gas -HA L (T L -T amb ) heat transfer with ambient through the wall

+ ṁL c p,L (T amb -T L ) mass flow (2.14)
Where m is the mass, c v and c p are the specific heat capacity at constant volume and constant pressure, respectively, h is the heat transfer coefficient, A is heat transfer surface area, H is the global heat transfer coefficient and ṁ = dm dt is the mass flow rate. The subscripts G, L and amb represent the gas phase, the liquid phase and the ambient, respectively.

The variation of wall temperature in contact with the gas phase or liquid phase can be calculated as:

m W,G c p,W dT W,G dt change of wall energy = h i,G A i,G (T G -T W,G )
heat transfer between gas and inner wall

-h o A o,G (T W,G -T amb )
heat transfer between outer wall and ambient

(2.15) m W,L c p,W dT W,L dt change of wall energy = h i,L A i,L (T L -T W,L )
heat transfer between liquid and inner wall

-h o A o,L (T W,L -T amb )
heat transfer between outer wall and ambient (2.16) Where the subscripts i, o and W represent the inner surface, the outer surface and the wall, respectively. The global heat transfer capacity HA G or HA L is calculated based on the sum of effective thermal resistances according to the convection on the inner and outer wall surfaces and the conduction through the wall.

HA G = 1 1 h i,G A i,G + th W λ W A ave,G + 1 hoA o,G
(2.17)

HA L = 1 1 h i,L A i,L + th W λ W A ave,L + 1 hoA o,L (2.18)
where th W is the thickness of the wall and the subscribe ave stands for average.

Kermani & Rokni [START_REF] Arjomand | Heat transfer analysis of liquid piston compressor for hydrogen applications[END_REF] developed a LP model based on the first law of thermodynamics to predict the evolution of temperature and pressure. The amount of heat transferred from the gas (hydrogen) towards the wall was estimated according to different heat transfer correlations. The results showed small changes in hydrogen temperature, between 0.2% and 0.4%, compared to the adiabatic case. Patil et al. [START_REF] Vikram | Experimental investigation of heat transfer in liquid piston compressor[END_REF] experimentally tested a LP compressor with different chamber materials and various stroke times of compression. It was observed that the heat transfer rate increased with the increasing piston velocity but at the expense CHAPTER 2. REVIEW ON LIQUID PISTON TECHNOLOGY FOR COMPRESSED AIR ENERGY STORAGE of additional compression work and a higher gas temperature. The convective heat transfer coefficient was observed to be high at the initial phase of the compression but then decreased to a stable value. In fact, the interior convective heat transfer coefficient h i,G for the gas phase is a key parameter that should be carefully estimated in the modeling of heat transfer behaviors of LP. For engineering applications, it is usually a convenient way to present and compare the results in terms of the dimensionless Nusselt number N u. Considering that h i,G varies along the piston (tube), it is usually the average value of h i,G used to calculate the N u:

N u D = 1 L 0 L 0 0 D λ h i,G (L)dL = h i,G D λ (2.19)
Where λ is the thermal conductivity of gas. Many correlations have been proposed to estimate the heat transfer performance of the working gas, mainly for a reciprocating piston application or for a compression fluid in tubes. In these correlations the N u number is usually written as a function of the Reynolds number Re and the Prandtl number P r, in a generic form shown in Eq. 2.20:

N u = ARe a P r b µ µ 0 c (2.20)
Where A, a, b and c are constants fitting different working conditions. µ and µ 0 are the fluid viscosity at fluid's average temperature and at the heat transfer boundary surface temperature, respectively. Some widely used empirical N u -Re correlations are listed in Table 2.3. Among them, Dittus & Boelter [START_REF] Dittus | Heat transfer in automobile radiators of the tubular type[END_REF] correlation was proposed for fully developed turbulent flow in circular tubes, with smooth surfaces and moderate temperature difference between the wall and fluid mean temperature. Sieder & Tate [START_REF] Sieder | Heat Transfer and Pressure Drop of Liquids in Tubes[END_REF] extended the Dittus & Boelter correlation by taking the wall effect into account for the condition of large temperature difference between the wall and the fluid. Hamilton [START_REF] Hamilton | Extensions of mathematical modeling of positive displacement type compressors[END_REF] proposed a new correlation more adapted to the compressor chamber, especially for L/D ≤ 60. Likewise, Annand [START_REF] Annand | Heat Transfer in the Cylinders of Reciprocating Internal Combustion Engines[END_REF], Adair et al. [START_REF] Adair | Instantaneous Heat Transfer To the Cylinder Wall in Reciprocating Compressors[END_REF], Liu & Zhou [START_REF] Liu | Heat Transfer Between Gas and Cylinder Wall of Refrigerating Reciprocating Compressor[END_REF] and Hsieh & Wu [START_REF] Hsieh | Experimental investigation of heat transfer in a highpressure reciprocating gas compressor[END_REF] proposed other correlations for the reciprocating compressors based on their proper experimental data.

More recently, Van de Ven & Li [START_REF] Van De Ven | Liquid piston gas compression[END_REF] numerically compared the evolution of Re, N u and gas temperature in a single stage LP and in a reciprocating solid piston. Their results showed that the LP could have significantly lower N u but much higher h i,G than that of the solid piston, mainly due to a smaller diameter of the compression chamber. Moreover, the LP can reach a near-isothermal compression process which can hardly be achieved by the solid piston. Two N u-Re correlations were proposed, one for laminar flow regime (Eq. 2.21) and the other for turbulent flow regime (Eq. 2.22) [START_REF] Cengel | Heat Transfer, A Practical Approach[END_REF].

N u = 0.664Re 1/2 P r 1/3 (2.21)

N u = 0.023Re 0.8 P r 0.3 (2.22)
Neu et al. [START_REF] Neu | Experimental study of convective heat transfer during liquid piston compressions applied to near isothermal underwater compressed-air energy storage[END_REF] experimentally studied a LP prototype for CAES. A set of configurations were tested with the column diameter D varying from 30 mm to 100 mm and the piston length L 0 from 2 m to 6 m. Based on the data of 73 compression tests, N u-Re correlation has been proposed for the laminar flow regime (Eq. 2.23) and for the turbulent flow regime (Eq. 2.24), respectively. Good agreement could be achieved between the correlation prediction and the experimental data [START_REF] Neu | Experimental study of convective heat transfer during liquid piston compressions applied to near isothermal underwater compressed-air energy storage[END_REF]. The transition from N u lam to N u turb could be determined by using the non-dimensional formula shown in the previous section (Eq. 2.12).

N u lam =6.

Heat transfer enhancement technologies for LP

The description in the above sections clearly indicates that how to handle the heat transfer to approach an isothermal gas compression and expansion is actually a key issue to increase the Liquid Pistons (LPs) efficiency and the RTE of Compressed Air Energy Storages (CAESs). Therefore, different techniques for HTE have been proposed and implemented within the LPs. This section introduces several representative types of HTE with special focus on the concept, the efficiency increase and the related N u-Re correlation.

Spray injection

Spray injection concept can be utilized to increase the heat transfer during the compression/expansion process where a high pressure spray of water droplets is injected into the top of LP column. The water droplets having a high heat capacity, scattered in the air, can absorb heat from the air during compression or inject heat into air during expansion, as explained in Fig. 2.7. The mean droplet diameter (MDD), the mass loading (ML defined as the ratio of water spray mass injected into a cycle and the mass of dry air drawn into the piston chamber) and the injected spray temperature are key operational parameters that determines the efficiency of this thermal management concept.

Qin & Loth [START_REF] Qin | Liquid piston compression efficiency with droplet heat transfer[END_REF] studied the efficiency of LP with droplet spray through a 1-D simulation validated with experimental data. Two ways of spray injection were considered: premixed injection upstream of the intake valve and direct injection within the compression chamber. For the direct injection, the nozzles were installed inside the piston to inject and uniformly distribute the liquid droplets. For the premixed injection, the nozzle was placed upstream the intake valve. Various conditions were tested including a CR ranging from 6 to 12, a MDD of 10, 20, 50 and 100 µm and a ML of 0.5, 1, 2 and 5. Their simulation results showed that the direct injection concept of spray cooling (CR=10; MDD=20 µm; ML=1) could reach a η c higher than 95% compared to about 71% for adiabatic compression. When ML increased to 5, the η c could be further augmented up to 98%, as shown in Fig. 2.8. A N u-Re correlation (Eq. 2.25) based on Ranz-Marshall form [START_REF] We Ranz | Evaporation from drops: Part 2[END_REF] was used to estimate the heat transfer characteristics of LP with liquid spray.

N u = 2 + 0.6Re 1/2 P r 1/3 (2.25)

With a similar configuration, Zhang et al. [START_REF] Zhang | Numerical Study of a Quasi-isothermal Expander by Spraying Water[END_REF] only obtained a 10% gas temperature abatement compared to the adiabatic case, mainly due to a small L/D ∼ 1. Guanwei et al. [START_REF] Guanwei | Micron-sized water spray-cooled quasi-isothermal compression for compressed air energy storage[END_REF] has built an experimental setup to test the spray cooling technique for LP. Experiments have been done in the summer and in the winter with different ambient temperatures to compare the effects of the nozzle diameter (0.1 -1 mm) and spray droplet diameters (10 -100 µm) on the heat transfer rate and thus on the η c . The best η c (92,8%) has been achieved with CR=2 and D = 0.3mm nozzle diameter. Further increasing the nozzle diameter required a higher work to generate the spray due to higher mass of water injected in the piston.

Odukumaiya et al. [START_REF] Odukomaiya | Thermal analysis of near-isothermal compressed gas energy storage system[END_REF] numerically and experimentally studied the effect of liquid spray for enhancing the heat transfer. Three configurations were tested as shown in Fig. Droplet Diameter (MDD)=20 µm). Adapted from [START_REF] Qin | Liquid piston compression efficiency with droplet heat transfer[END_REF] compression and a fall of 15 K during expansion. The introduction of spray droplets could attenuate the temperature variation, to 13 K (for configuration 2) and 3 K (for configuration 3), respectively. Nevertheless, this active heating/cooling for regulating the temperature of spray via an external heat exchanger needs extra energy input. Note that the same correlation of Eq. 2.25 has been used to model the heat transfer. The impacts of other parameters, i.e., the compression time, the injection pressure and the spray angle on the heat transfer and the efficiency have also been studied by Patil et al. [START_REF] Vikram | Efficiency improvement of liquid piston compressor using metal wire mesh for near-isothermal compressed air energy storage application[END_REF][START_REF] Vikram | Efficiency Improvement of a Liquid Piston Compressor Using Metal Wire Mesh[END_REF]. Their results showed the positive effect of injection pressure when it increased from 10 psi (69 kPa) to 30 psi (483 kPa), but the benefit became marginal beyond 30 psi. The nozzle angle at 60 • was shown to be the most efficient configuration whereas the compression time had no significant impact.

The spray injection concept has proven its efficiency in both multi-stages and one stage compressions especially for obtaining a trade-off between power density and efficiency. For one stage compression, the tested CR is relatively low (maximum 10), which are not adapted for an industrial application. Another issue is the pumping power for the spray injection at high pressures and high flow rates [START_REF] Odukomaiya | Thermal analysis of near-isothermal compressed gas energy storage system[END_REF]. The trade-off between all these factors should be determined carefully.

Porous media

Another largely used measure to enhance the heat transfer is the porous media inserts. The used porous media could have several forms (e.g., structured interrupted-plates or unstructured foams) and with the use of different materials and densities (Fig. 2.9). They can cover the whole compressed/expanded gas domain or just a part of it. It has been reported that an CHAPTER 2. REVIEW ON LIQUID PISTON TECHNOLOGY FOR COMPRESSED AIR ENERGY STORAGE insert close to the piston top is the most efficient [START_REF] Wieberdink | Effects of porous media insert on the efficiency and power density of a high pressure (210 bar) liquid piston air compressor/expander : An experimental study[END_REF]. On one hand, the porous inserts serve as a good meida for heat transfer between the gas and liquid. On the other hand, they can reduce the time needed to complete the compression/expansion process, increasing thereby the power density for a given thermodynamic efficiency, or augmenting the thermodynamic efficiency at a fixed power density. Note that the power density ( E tA ) is defined as the energy stored or recovered during a process (E), divided by the total time it took to complete the process (t) and the surface area of the piston column (A). It has been reported that in compression, the porous inserts could increase the power density by 39-fold at 95% efficiency, or augment the efficiency by 18% at 10 5 W/m 3 power density [START_REF] Zhang | Numerical analysis of heat exchangers used in a liquid piston compressor using a one-dimensional model with an embedded two-dimensional submodel[END_REF]. In expansion, the power density could be increased threefold at 89% efficiency or the efficiency be increased by 7% at 1.5 × 10 5 W/m 3 power density [START_REF] Zhang | Thermal analysis of a compressor for application to Compressed Air Energy Storage[END_REF]. Zhang et al. [START_REF] Zhang | Numerical investigation of metal-foam filled liquid piston compressor using a two-energy equation formulation based on experimentally validated models[END_REF] compared several heat transfer correlations to estimate the heat transfer in their studied geometry with porous media inserts. These correlations were proposed by many researchers, including Wakao & Kaguei [START_REF] Wakao | Heat and mass transfer in packed beds[END_REF], Kuwahara et al. [START_REF] Kuwahara | A numerical study of interfacial convective heat transfer coefficient in two-energy equation model for convection in porous media[END_REF], Zukauskas [START_REF] Zukauskas | Convective heat transfer in cross flow[END_REF], Nakayma et al. [START_REF] Nakayama | A study on interstitial heat transfer in consolidated and unconsolidated porous media[END_REF], Kamiuto & Yee [START_REF] Kamiuto | Heat transfer correlations for open-cellular porous materials[END_REF] and Fu et al. [START_REF] Fu | Measurement and correlation of volumetric heat transfer coefficients of cellular ceramics[END_REF] (cf. details in Table 2.3). Among them, the Kamiuto & Yee correlation [START_REF] Kamiuto | Heat transfer correlations for open-cellular porous materials[END_REF] fitted the best their experimental results with 10 PPI (Pores Per Inch) insert, and was then chosen for a further numerical study [START_REF] Zhang | Design of an interrupted-plate heat exchanger used in a liquidpiston compression chamber for compressed air energy storage[END_REF]. A new correlation has been proposed based on the obtained simulation results (Eq. 2.26):

N u = 9.700 + 0.0879Re 0.792 Dh P r 1/3 (2.26)

Later on, an improved correlation has also been proposed (Eq. 2.27) [START_REF] Zhang | Thermal analysis of a compressor for application to Compressed Air Energy Storage[END_REF]:

N u = 8.456 + 0.325Re 0.625 Dh P r 1/3 (2.27)

Where Re Dh is the Reynolds number based on the hydraulic diameter D h of the piston column.

Besides, the Stanton number (St) and a dimensionless heat flow rate term (Q * ) were used to characterize the effect of porous inserts.

St = h i,G ρc p U 2
(2.28)

A linear relation St versus Q * -1 has been found [START_REF] Zhang | Thermal analysis of a compressor for application to Compressed Air Energy Storage[END_REF], written as:

-St = a(Q * -1 ) + b (2.29)
Where a, b are coefficients depending on the compression velocity.

Yan et al. [START_REF] Yan | Experimental study of heat transfer enhancement in a liquid piston compressor/expander using porous media inserts[END_REF] conducted an experimental investigation to measure the effects of porous media inserts on the heat transfer, the power density and the compression/expansion efficiencies. Compression experiments were run using 5 different inserts and a baseline (Fig. 2.9). Power densities ranged from 3.8 × 10 3 W/m 3 to 1.82 × 10 5 W/m 3 while similar constant flow profiles were imposed. In compression, porous media could increase the power density by a factor of 39 at a constant efficiency, or increase the efficiency by as much as 18% to reach a maximum of 95% (Fig. 2.10). For the expansion, the results showed an increase in power density by a factor of 3 for a fixed efficiency, and an increased process efficiency by as much as 7% up to 90%. Fig. 2.10 (b) shows the expansion efficiency vs. the normalized power with or without porous inserts, indicating that an efficiency improvement of 5% could be achieved at a certain mean power density. More recently, a LP structure with a two layer porous medium was proposed by Ren et al. [START_REF] Ren | Experiments on air compression with an isothermal piston for energy storage[END_REF], allowing enhanced air cooling so as to approach the isothermal compression. It was reported that this configuration could also reduce the pressure and the compression work needed.

A thermodynamic and heat transfer analysis was conducted by Khaljani et al. [START_REF] Khaljani | Thermodynamic and heat transfer analysis of a Liquid Piston Gas Compressor (LPGC)[END_REF] to evaluate the use of porous media as an effective for thermal management. The porous medium has been modeled as metal parallel plates inserted in the piston column. The results showed that the inclusion of five and nine plates could reduce the final temperature of the air by 15 K and 20 K, respectively, compared to the base case without plate inserts. The porous insert also had a significant influence on the wall temperature and on the flow pattern evolution. With 9 plates insert, the η c could be increased by 2%, going from 86.91% for the base case to 88.91%. Nevertheless, the energy storage capacity would decrease, since a part of the piston volume was occupied by the porous media, leaving smaller space for the gas phase.

In brief, porous media insert could be considered as the most widely used HTE technology for LP, adapted to different CR for both the compression and the expansion. It can offer high power densities while maintaining the efficiency and can be easily combined with the other HTE concepts as will be introduced in following sub-sections. Further researches could be focused on finding the optimal geometries of the porous media by using shape or topology optimization methodologies [START_REF] Boichot | A genetic algorithm for topology optimization of area-to-point heat conduction problem[END_REF].

Optimal trajectories

Another proposed HTE concept consists in optimizing the manner of air compression from the view point of finite-time thermodynamics by choosing the best trajectory (ζ) (Fig. 2.4), with the purpose of significantly improve the energy storage/generation capacity. Moreover, it tries to find the trade-off between the HTE and the pressure drop increase by imposing a pressure-linked flow rate, i.e., a high flow rate at the beginning of the compression/expansion but gradually decreasing with time.

Saadat et al. [START_REF] Saadat | Optimal trajectories for a liquid piston compressor/expander in a Compressed Air Energy Storage system with consideration of heat transfer and friction[END_REF] proposed a general numerical approach to determine the optimal compression profile by using more general heat transfer models. The viscous friction and system constraints were also considered in the optimization process. A fully developed laminar pipe flow model was used to estimate the heat transfer coefficient inside cylinders with high L/D >> 1 (Eq. 2.30):

h = 3.66 λ D (2.30)
An average Nusselt number N u ave was also proposed based on Re air,ave for cylinders with small L/D ∼ 1 where the flow is supposed to still be laminar (Eq. 2.31):

N u ave = 0.664Re Where U air,ave (t) is the air average velocity and L(t) is the liquid/air interface position. Based on the Pareto optimization procedure, their results [START_REF] Saadat | Optimal trajectories for a liquid piston compressor/expander in a Compressed Air Energy Storage system with consideration of heat transfer and friction[END_REF] showed that the optimal compression profile increased the storage power up to 10% for a η c of 90% compared to the A-PI-A (Adiabatic-PseudoIsothermal-Adiabatic process) and more than 40% when compared to a conventional compression profile such as sinusoidal or constant speed compression trajectories (Fig. 2.11). Several other trajectories were tested and compared. Note that the study was limited to a relatively small CR=10, the effect of optimal trajectory could be expected to be more significant at higher CR and with more complex geometries.

Honeycomb geometry

The honeycomb concept consist in substituting the large piston column by many long and thin tubes, the liquid flowing through these tubes to compress air. The thin tubes could have a circular geometry (Fig. 2.12a), or as hexagonal cells stacked in layers (Fig. 2.12b).

The main advantages of this concept are the easy implementation of the inserts which gives a higher power density and an increased L/D in each cell which enhances the heat transfer. Zhang et al. [START_REF] Zhang | Heat transfer in a long, thin tube section of an air compressor: An empirical correlation from CFD and a thermodynamic modeling[END_REF][START_REF] Zhang | Numerical analysis of heat exchangers used in a liquid piston compressor using a one-dimensional model with an embedded two-dimensional submodel[END_REF] numerically studied a LP compressor inserted with honeycomb geometry by CFD approach. Each tube channel had a L/D ratio of 100. The liquid (moving boundary assumed as water) had a constant velocity and the walls had a constant temperature. In their study a transient Re Dh number ranging from 130 to 1300 has been calculated based on the diameter (of the largest cylinder), the mean flow velocity, and a transient den- honeycomb matrix [START_REF] Zhang | Numerical analysis of heat exchangers used in a liquid piston compressor using a one-dimensional model with an embedded two-dimensional submodel[END_REF] sity. Moreover, a N u-Re correlation was proposed for predicting the heat transfer based on the local Darcian velocity Re lc , where l c is the characteristic length [START_REF] Celata | Laminar Flow Friction and Heat Transfer in Non-Circular Ducts and Channels Part II-Thermal Problem[END_REF]:

N u = 2f r l c L d Re lc P r 1/2 (2.33)
Where f r is the friction factor defined as Eq. 2.34 :

f r = 0.564 1 + (0.664P r 1/6 ) 9/2 2/9
(2.34)

Zhang et al. [START_REF] Zhang | Numerical analysis of heat exchangers used in a liquid piston compressor using a one-dimensional model with an embedded two-dimensional submodel[END_REF] has compared two concepts, a honeycomb insert and a porous media interrupted-plate heat exchanger, on the efficiency improvement of LP. The impacts of different parameters (material, length scale of the plate elements, specific surface area and porosity) on the temperature evolution were studied. A honeycomb insert with higher power density showed better air temperature abatement, and thus better efficiency. Nevertheless, the difference of η c among all tested configurations was below 4% (between 60% and 64%).

In summary, the honeycomb concept with mini-tubes insert can take the maximum benefit from slow compression/expansion using a LP. They are easy to implement, but more investigations are still required to further investigate the effect of the materials and the CR on the heat transfer. Studies on the expansion will also help the application of this HTE concept.

Hollow spheres

The hollow sphere concept consists in adding a layer of floating spheres at the gas-liquid interface, which absorb the heat from the gas and transfer it to water. These hollow spheres move along with the interface during the compression or expansion. In this way, the overall heat transfer rate would be increased without adding much the compression work. ). An analytical model has been developed for the compression, by assuming the temperature-dependent physical properties of air. The viscosity has been calculated by using the Sutherland formula: (2.36)

µ(T ) = µ 0 T T 0 1.
Ra = gρ 2 c T λµ (T -T W )D 3 (2.37)
where Ra is the Rayleigh number, D the diameter of the spheres, and c the specific heat of air.

c(T ) = 1002.5 + 275 × 10 -5 (T -200) 2 (2.38)
The thermal conductivity of air λ is calculated by:

λ(T ) = 0.02624 T 300 0.8646 (2.39)
Many influencing parameters were tested including the spheres diameter, the stroke time, the CR and the sphere's material (for the experimental part). The results showed clearly the effect of the material of the spheres on the temperature abatement compared to the base case. HDPE (High-Density Polyethylene) and SiC (Silicon Carbide) showed better performance than PP (Polypropylene) (Fig. 2.14). Moreover, the smaller spheres diameter (D=6.5 mm) showed better η c than the larger one (D=9 mm) because of the higher total heat transfer surface area. A maximum temperature abatement of 32 K could be achieved by a fast compression with hollow spheres insert.

Although the study was limited to a small CR and a limited range of spheres' diameter/material, this concept has showed its effectiveness and easy implementation. Yet, more investigations are still needed to further test this concept under high pressures, high CR and high flow rates conditions.

Aqueous foam

This concept consists in generating foam (with water additives) from the bottom of the piston which rises to the liquid/gas interface and increases the heat transfer surface between the liquid column and the compressed gas. This novel technique has been proposed and patented by SustainX [START_REF] Bollinger | Demonstration of Isothermal Compressed Air Energy Storage to Support Renewable Energy Production[END_REF] and studied experimentally by Patil & Ro [START_REF] Vikram | Experimental study of heat transfer enhancement in liquid piston compressor using aqueous foam[END_REF] for temperature abatement in a LP compressor. Compared to the droplet spray, the mixture of air and water as a homogeneous foam had several advantages, including the higher contacting surface area which can be maintained during the entire compression or expansion stroke. With the use of aqueous foam, it has been observed that the rise of air temperature could be reduced especially in the end of the compression process with higher temperature gradients, i.e., about 7 to 20 K at a CR=2.5. This leaded to a η c increase of 4-8% to 86% of isothermal η c . Moreover, the use of several foam generators (up to 4) could further increase the η c up to 90 %. The operating principle and the influence of the aqueous foam on the gas flow dynamics are illustrated in Fig. 2.15.

Patil & Ro [START_REF] Vikram | Experimental study of heat transfer enhancement in liquid piston compressor using aqueous foam[END_REF] used Skelland correlation (Eq. 2.40) [START_REF] Skelland | Heat transfer in turbulent non-Newtonian flow[END_REF] and Attia correlation (Eq. 2.41) [START_REF] Attia | Convective heat transfer in foams under laminar flow in pipes and tube bundles[END_REF] to model the heat transfer with aqueous foam in the LP. Nevertheless, none of them could give appropriate results due to the complexity of the physical phenomenon and the geometry.

N u x = 1.41 3k+1 4k 1/3 D h ReP r x 1/3 and N u ∞ = 8(5k+1)(3k+1) 31k 2 +12k+1 (2.40)
where N u x and N u ∞ are the Nusselt number in the entry region and in thermally fully developed region, respectively. k is the flow behavior index and x is the length from the entrance.

N u f o,x = a C (χ) -1 6 (Ca f o ) 1 8
3k + 1 4k

1/3 D h Re f o P r f o x 1/3 (2.41)
where a is a constant to be determined from experimental data, C(χ) is an empirical function dependent on surfactant mass fraction ξ in weight percentage and Ca f o is the capillary number. The subscript f o denotes the foam. AIR ENERGY STORAGE

The presence of foams in the gas domain inside the LP shows relatively high heat transfer rate with a low mass flow rate. Therefore, the ratio of heat transfer to coolant density is higher for foams. However, the accumulation of residual foams after a few cycles could change the heat transfer characteristics and flow dynamics inside the LP and could lead to the corrosion in some parts of the system. Effects of cyclic operation and variability of foams geometries have to be further studied so as to obtain more experimental feedback. Better scalability of the foam-based HTE could be ensured by further studies focusing on the use of microfoams, controlled macrofoams generation and under higher CR operations.

Wire mesh

Metal wire meshes can also be installed in the piston column in contact with both the gas and the liquid to enhance the heat transfer (Fig. 2.16a). For cylindrical column, an Archimedean spiral form could be used, facilitating the heat transfer in both the axial and radial directions (Fig. 2.16b). This concept could in fact be considered as a special case of porous media inserts.

Patil et al. [START_REF] Vikram | Efficiency improvement of liquid piston compressor using metal wire mesh for near-isothermal compressed air energy storage application[END_REF] compared experimentally the effect of wire mesh on the pressure, the temperature and the efficiency of a LP compressor. Two types of wire mesh were tested, one in copper and another in aluminium (Fig. 2.16) (b). Their results showed a significant temperature abatement (about 25 K) by adding the wire mesh. Meshes with larger wire size (0.7112 mm or 3.5 wire/cm) were more efficient than that with smaller wire size (0.4572 mm or 6 wire/cm). The compression speed also had an important influence on the final gas temperature. At medium velocity condition (4 s compression time), about 28 K temperature abatement could be achieved (Fig. 2.17). In general, the isothermal η c could be increased from 82-84% to 88-92% by using the wire mesh [START_REF] Vikram | Efficiency improvement of liquid piston compressor using metal wire mesh for near-isothermal compressed air energy storage application[END_REF].

Only two materials (copper and aluminium) have been tested while more systematic investigations on many factors (e.g., wire diameters, mesh sizes, forms of the spiral and mesh position) are still needed, so as to achieve an optimal mesh design and to obtain the maximum efficiency. Also, the utility of metal wire meshes for high CR should be studied to ensure that the obtained efficiencies will also be maintained at high power levels.

Optimal geometry of the piston column

Another method to augment the compression/expansion efficiency is to optimize the shape of the piston column instead of the conventional cylinder one with uniform cross-section area, shown as an example in Fig. 2.18. Zhang et al. [START_REF] Zhang | Design analysis of a liquid-piston compression chamber with application to compressed air energy storage[END_REF] made a CFD-based design analysis for a LP compression chamber by varying five geometrical parameters, including the inlet radius R i , the maximum radius R m , the top cap radius R t , the location of the maximum radius X m with respect to the inlet and the column length L (Fig. 2.18). 16 designs were proposed and Generally speaking, parametric study or topology optimization (geometry, piston velocity, materials, CR) are the best way to better understand and optimize the processes to help achieve the highest efficiencies. While it is costly to do such studies experimentally, numerical modeling based on CFD simulations could be more adapted and play an important role in the future. Once the optimal geometry of the column is determined, another problem is the fabrication difficulty and cost. Therefore, trade-off needs to be achieved for this concept for large-scale industrial application considering the computational cost, the efficiency improvement, the production cost, the material resistance and the installation in a real site.

Summary of HTE concepts for LP

Patil et al. [START_REF] Vikram | Efficiency Improvement Techniques For Liquid Piston based Ocean Compressed Air Energy Storage[END_REF] have compared several HTE concepts used in LP for CAES application. The comparison was done on the RTE value through a typical configuration of a global CAES system, composed of 4 LPs having a L/D = 10 and a CR = 6. The total stored energy was 2 MWh [START_REF] Vikram | End-to-end efficiency of liquid piston based ocean compressed air energy storage[END_REF]. By improving nothing more than the compression/expansion efficiency through thermal management, the RTE of the CAES system could rise from 45% to up to 62% (Fig. 2.20). The RTE improvement is different by using different HTE concepts, i.e., 5% by the optimal trajectory concept, 7% by hollow spheres insert, and more than 17% by spray cooling or by porous media inserts. Further improvement of LP efficiency could be expected by combining two or several HTE methods, stepping towards near-isothermal compression/expansion. For example, Ahn et al. [START_REF] Barah Ahn | Temperature Abatement using Spray Injection and Metal Wire Mesh in Liquid Piston Compressor for Ocean Compressed Air Energy Storage Application[END_REF] combined wire mesh insert with the droplet spray but the effect of each HTE measure could not be identified clearly. Further efforts should thereby be devoted to the development and implementation of combined HTE methods. Table 2.2 recapitulates the studies on the HTE technologies in LP with their main results, conditions, advantages and shortcomings. Patil & Ro [START_REF] Vikram | Design of Ocean Compressed Air Energy Storage System[END_REF] 2020

N u f o,x = a C (χ) -1 6 (Ca f o ) 1 8 3k+1 4k 1/3 D h Re f P r f o x 1/3
Aqueous foam insert

Conclusion and perspectives

This review chapter is dedicated to gather and survey the advances of the LP technology for CAES. The basic principles, the thermodynamic background, the flow and heat transfer characteristics related to LP are presented in detail. Especially the technological state-ofthe-art on the HTE concepts implemented in LP has been surveyed, aiming at approaching the isothermal compression/expansion. The main findings may be summarized as follows.

• Compared to solid piston, LP can allow high L/D geometry and a flexible compression/expansion speed. Moreover, it is a promising way to achieve a near-isothermal process, which will have a bright foreground for I-CAES application.

• The fluid flow and heat transfer inside the LP during the compression or expansion operation are rather complex and difficult to characterize, especially for the transition and turbulent regimes.

• Many N u-Re correlations were employed to predict the heat transfer of the compression/expansion processes inside the LP, as summarized in Table 2.3. Almost all of them are based on the average value of N u number or convective heat transfer coefficient.

• Different concepts were proposed to enhance the heat transfer in the LP, including liquid spray, wire mesh, porous media, optimal trajectory, hollow spheres and optimal geometry of the piston column. By using these HTE techniques the compression/expansion efficiency of the LP could be significantly improved up to 10%, leading to the RTE increase of the CAES system. Each HTE concept has their own advantages and short comings while the combination of two or several concepts seems more promising. By implementing the proper HTE technique, the LP in a CAES system can achieve a 95% compression/expansion efficiency.

Some scientific and technological barriers remain to be overcome for the widespread industrial application of LP, which are also the key issues and challenges of the current research and development:

• Fine models that describe the evolution of the flow & heat transfer inside a LP during compression/expansion processes are still lacking. Therefore, systematic experimental and numerical studies are especially needed on the local coupled fluid flow and heat transfer behaviors, for better characterization, understanding and mastering the transient gas-liquid two phase flow transport phenomena.

• Critical geometric and operational parameters (flow rate, CR, compression time, geometry, energy density, etc.) and the suitable combinations of them remain to be identified. This requires the development and application of advanced modeling and optimizing algorithms and experimental verification.

• Most of the LP studies are still on the laboratory scale with small LP dimensions and small CR values. The field testing at pilot or industrial scale with realistic CR are still rare but fully necessary for the possible application of LP in CAES systems.

• Further investigations on the new or combined HTE concepts, their testing and validation at pilot-scale prototypes are still required.

• Finally, studies are still required on full cyclic compression/expansion cycles considering all the auxiliary components (pumps, valves, etc.). Detailed thermodynamic and techno-economic analyses can also help determine the optimal operation conditions and adapted control strategies especially at pilot or industrial level.

Chapter 3

Attempts to model the liquid piston flow using Lattice Boltzmann Method

Introduction

The LBM (Lattice Boltzmann Method), commonly abbreviated to LBM, is a relatively new approach for flows simulations and contrasts with the traditional (Navier-Stokes) approach to CFD by adopting a bottom-up approach to fluid modeling. Further details of the method can be seen in the Appendix 8.

Research barriers and main contributions

The liquid piston is a crucial component of REMORA system and it is important to study the flow behavior and optimization of this system. The LBM (Lattice Boltzmann Method) is first chosen to model the liquid piston, due to its high performance for parallel computation, topological and shape optimization [START_REF] Dugast | Topology optimization of thermal fluid flows with an adjoint Lattice Boltzmann Method[END_REF][START_REF] Wang | Heuristic optimality criterion algorithm for shape design of fluid flow[END_REF] and its innovative character [START_REF] Boichot | A genetic algorithm for topology optimization of area-to-point heat conduction problem[END_REF][START_REF] Boichot | A simple Cellular Automaton algorithm to optimise heat transfer in complex configurations[END_REF].

In this chapter, simulation of the liquid piston is performed to study the physics of it using LBM. There are multiple challenges in performing the simulation of the liquid piston.

Multi-component (multi-fluid) simulation is the first one that must be studied. Most of the approaches to compute the multi-component multi-phase (MCMP) are unstable for highdensity ratios. Different analogies and methodologies have to be included to improve the stability of these solutions. Steps have to be taken to increase the stability of the system for high viscosity ratio as well. Other complications are to study the compressibility, the temperature evolution and the turbulence in the flow. Different approaches to implement the multiphase simulation are studied and the suitable method is chosen. Firstly , Shan & Chen (S-C) method [START_REF] Shan | Lattice Boltzmann model for simulating flows with multiple phases and components[END_REF] [START_REF] Shan | Multicomponent lattice-Boltzmann model with interparticle interaction[END_REF] is studied and implemented to get the required results. The S-C method is stable for density ratio of around 10, but the simulation of the liquid piston has to include the interaction of water and air having the density ratio of around 1000. Different analogies like higher-order non-ideal equation of state has to be implemented. Proper forcing scheme should also be implemented to increase the stability of the simulation. These implementations are performed on a simple static drop suspension case and later applied on to the flow in a channel. Later, the free energy models (Lee-Lin [START_REF] Lee | Eliminating parasitic currents in the lattice Boltzmann equation method for non ideal gases[END_REF] and Inamuro [START_REF] Inamuro | Lattice Boltzmann methods for viscous fluid flows and for two-phase fluid flows[END_REF]) are implemented in, respectively, 2D and 3D in-house codes. Static droplet and two droplets collision are simulated to compare and validate the implemented models. Finally, a Volume Of Fluid like method coupled to LBM is improved based on free-surface models to enable a two-phase flow simulation. Test cases of incompressible and compressible flows are computed.

Different Methods to model Multi-component Multi-phase in LBM

Interaction between different fluids (multi-component) is a common phenomenon in nature, where one fluid comes in contact with another giving rise to different physical phenomena.

Interaction between different phases of the same fluid also has to be considered due to the co-existence of phases at the same time. To model these conditions of MCMP, there are different methods in LBM.

• The color-Method: This method was developed by Gustensen & Rothman [44][138] which was the first method introduced. In this method, the fluid components are represented by the colour function, where each of the components evolves using its own LBM equations and the interaction of the components using the local colour gradient.

• The Pseudo-potential method or Shan & Chen (S-C) method: This method was developed by Shan & Chen [START_REF] Shan | Lattice Boltzmann model for simulating flows with multiple phases and components[END_REF] which is based on the interaction forces between different fluid phases and the equation of state of the fluid.

• Free-energy method: Simulation of multi-fluid by modelling the interaction forces using Free-energy potential was developed by Swift [START_REF] Michael | Lattice Boltzmann Simulation of Nonideal Fluids[END_REF]. Further development was made to this method by Inamuro [START_REF] Inamuro | An improved lattice Boltzmann method for incompressible two-phase flows with large density differences[END_REF] by introducing higher-order parameters for the interface of the fluids and also adding the effects of Poisson equation into LBM.

• VOF method: VOF (Volume of Fluid)-LBM coupling is a novel method based on the LBM and VOF for interface tracking. It is used initially for free surface flows.

There are many other methods which are based or derived on the ones mentioned above. They mainly deal with the interface treatment or the interaction forces between the fluids.

Firstly, the S-C method is modified to adapt it to the studied flow [START_REF] Shan | Lattice Boltzmann model for simulating flows with multiple phases and components[END_REF] as it is simple to implement and the separation of the phases and components are the effects of the particle LATTICE BOLTZMANN METHOD interaction forces, avoiding expensive interface tracking systems. The models were implemented in Palabos. Then free-energy models (Lee-Lin [START_REF] Lee | A lattice Boltzmann algorithm for calculation of the laminar jet diffusion flame[END_REF], Inamuro [START_REF] Inamuro | Lattice Boltzmann methods for viscous fluid flows and for two-phase fluid flows[END_REF]) which are adapted for multiphase flows with high density ratios are implemented in a in-house program in Fortran. Finally, a 3D VOF-LBM coupling model was adapted from the free-surface model in Palabos [START_REF] Thürey | Optimization and stabilization of LBM free surface flow simulations using adaptive parameterization[END_REF].

Shan-Chen method

Original S-C method

S-C method is based on Pseudopotential method where the phase and component separation effects are incorporated into the LBM equation in the form of inter-particle interaction forces. The forces are developed based on non-ideal Equation of state (EOS) which represent the multi-phase effects.

The separation of the phase and the component is due to the interaction forces between the molecules. There are Two parts of interaction forces in multi-component simulation, one is the interaction force between the same component and the other interaction force is between the different components. The interaction inter-particle forces within the same component (i, i) can be expresses as :

F i,i (x) ∼ = -c 2 s ψ i (x)G ii ∇ψ i (x) (3.1)
The interaction of inter-particle forces between the different components (i, j) can be expressed as:

F i,j (x) ∼ = -c 2 s ψ i (x)G ij ∇ψ j (x) (3.2)
ψ j (x) is the effective mass which is a function of the local density. "G" is the strength coefficient (Greens function) of interparticle force, where G < 0 represents attractive force and G > 0 represents repulsive force between the particles. Effective mass ψ k (x) of the S-C method is expressed as local density or the exponent form of local density as shown below:

ψ k (x) = ρ o 1 -exp -ρ k (x) ρ o (3.3)
To include the interaction forces developed between the different phases and components in the S-C equation is included using the velocity shift method 8.28. where u i is the macroscopic bulk velocity of the fluid, τ i is the relaxation time, ρ i is the density of the fluid and F i is all the forces acting on the particles such as interaction forces between the same component, different components, surface tension forces and external forces.

u eq i = u i + τ i F i ρ i (x) (3.4) LATTICE BOLTZMANN METHOD
the Macroscopic bulk velocity u i which is required to compute u eq i is computed with the sum of all the components in the simulation at that lattice point.

u = i 1 τ i N α=0 f α,i e α i 1 τ i N α=0 f α,i (3.5)
In the original S-C method, the maximum density for the multi-component is in the range of ∼ 4 without any numerical error. As the density increases the numerical error increases resulting in numerical instabilities. The original S-C method is also unstable for high viscosity ratio. Different methods must be determined to increase the stability of the method for high-density ratio.

Implementation of new S-C approaches for high density ratio modeling

The density ratio of the S-C method can be increased by introducing the higher-order and more accurate non-ideal EOS, which can be incorporated into S-C method [START_REF] Yuan | Equations of state in a lattice Boltzmann model[END_REF]. Non-ideal EOS decreasing the error caused due to numerical instability and allows to model higher density ratios. Better forcing method has to be adapted to increase the high viscosity ratio and also the stability of the equations. there are other simulation problems such as high compression and proper utilization of the turbulence model to study the flow. Which will be studied once the proper results are obtained for high-density simulations.

The modified S-C method proposed in the literature [START_REF] Stiles | High density ratio lattice Boltzmann method simulations of multicomponent multiphase transport of H2O in air[END_REF] [10] [START_REF] Yuan | Equations of state in a lattice Boltzmann model[END_REF] to increase the density ratio of two-phases of the fluid will be implemented using Palabos code by modifying the original S-C method. The appropriate EOS is chosen to model the high-density ratio initially for a Single-Component Multi-phase (SCMP) using S-C model and velocity shift method for forcing will be used for a static case. Parameters are varied and the stable condition is studied for carrying out MCMP simulation. Once the stable solution of MCMP is obtained the simulation is performed in a channel where one fluid pushes another fluid.

Modified S-C Method for SCMP

The effective mass (ψ) has to be computed with higher-order non-ideal EOS to increases the stability of the numerical solution for high-density values. Here, we use the cubic non-ideal EOS where the effective mass depends on the critical pressure, temperature and density of the component.

Original S-C force is considered and the gradient of effective mass is computed using the cubic non-ideal EOS. The force F is computed between two different components i, j. G ij is the Greens function which controls the interaction force between the components and ψ is the effective mass of component i and component j.

F i,j (x) -c 2 s ψ i (x)G ij ∇ψ j (x) (3.6)
Interaction of inter-particle forces computed using the gradient of effective mass which is sixpoint numerical gradient dependent on nearest and next to nearest neighbor lattice. Where, ψ is the Effective mass, w 1 and w 2 are lattice weights.

∂ψ(i, j) ∂x = w 1 [ψ(i + 1, j) -ψ(i -1, j)] + w 2 [ψ(i + 1, j + 1) (3.7) -ψ(i -1, j + 1) + ψ(i + 1, j -1) -ψ(i -1, j -1)] (3.8) ∂ψ(i, j) ∂y = w 1 [ψ(i, j+1)-ψ(i, j-1)]+w 2 [ψ(i+1, j+1)-ψ(i+1, j-1)+ψ(i-1, j-1)-ψ(i-1, j-1)]
(3.9) The pressure can be computed using the non-ideal gas equation as shown below:

p(ρ) = ρc 2 s + 1 2 Gc 2 s ψ 2 (ρ) (3.10)
For implementation of effective mass of complex EOS the equation of original S-C method can be modified to include the non-ideal pressure equation as shown below:

ψ i (ρ) = 2 (p i -c 2 s ρ i ) c 2 s G ii (3.11)
The EOS gives the relation between the pressure, temperature and density. The pressure term p i can be computed with different EOS as shown by Yuan and Shaefer [START_REF] Yuan | Equations of state in a lattice Boltzmann model[END_REF].

Equation of State (EOS):

There are different EOS which can be used to reduce the instability and increase the density ratio. Few of the commonly used higher-order EOS are listed below:

Van der Waals (VdW) EOS

P = ρRT 1 -bρ -aρ 2 (3.12)

Carnahan-Starling(CS) EOS

P = ρRT 1 + bρ/4 + (bρ/4) 2 -(bρ/4) 3 (1 -bρ/4) 3 -aρ 2 (3.13)

Peng-Robinson(PR) EOS

P = ρRT 1 -bρ - aα(T )ρ 2 1 + 2bρ -b 2 ρ 2 (3.14) LATTICE BOLTZMANN METHOD
In VdW EOS the instability is not caused due to the spurious velocity but due to the EOS itself [START_REF] Yuan | Equations of state in a lattice Boltzmann model[END_REF]. Hence it is not suitable to model a liquid-vapor system. Whereas CS EOS will give better performance with a high-density ratio around 1000, wider temperature range and less spurious velocity compared to VdW [START_REF] Yuan | Equations of state in a lattice Boltzmann model[END_REF].

The Peng-Robinson EOS can also simulate density ratios higher than 1000 and reduces the spurious velocity around the interface similar to CS EOS. We choose PR EOS as it is stable for even higher density ratio than CS EOS. PR EOS uses the acentric factor allowing binodal curve adjustment for obtaining the better physical property of the particular fluid. [176] [24] [156] [START_REF] Bao | Lattice Boltzmann equation model for multi-component multi-phase flow with high density ratios[END_REF]. Where R is the gas constant, T is the temperature, ρ is the local density. Here, a and b are attractive and repulsive parameters.

The α(T ) of PR EOS can be computed using the equation given below. ω is the acentric factor and for water the value is ω = 0.344.

α(T ) = 1 + 0.37464 + 1.5422ω -0.26992ω 2 1 -T /T cr (3.15)
Here, a is the attraction parameter, b is the repulsion parameter which is computed based on the critical temperature T cr , critical pressure P cr and Gas constant R for the PR EOS.

a = 0.45724R 2 T 2 cr P cr (3.16) b = 0.0778RT cr P cr (3.17)

Modified S-C Method for MCMP

The simulation of the multi-component system such as water and air has to be implemented differently compared to the single-component multi-phase. The interaction force of one component is acted on the other component which must be modeled. In this case of simulation equation of water is considered as non-ideal gas, where PR EOS is used. The equation of air is considered as an ideal gas and original SC method is directly used to compute the effective mass.

The interaction force within water for different phases has to be considered, as the forces are grater, but for the interaction forces between the air can be ignored as the forces are week. The inter-component forces must be implemented to consider the interaction force between the air and water, that is, when computing the forces for air and also when computing for water components.

Here, let us consider the subscript 1 for water and subscript 2 for air. To compute the interaction force acting on the water by air and vice versa, the original S-C equation of effective mass ϕ is used as shown below:

ϕ 1 (ρ 1 (x)) = 1 -exp (-ρ 1 ) (3.18) LATTICE BOLTZMANN METHOD ϕ 2 (ρ 2 (x)) = 1 -exp (-ρ 2 ) (3.19)
The S-C force computed for water component is shown below:

F 1 (x) = -c 2 s G 11 ψ 1 (x)∇ψ 1 (x) -c 2 s G 12 ϕ 1 (x)∇ϕ 2 (x) (3.20)
Where, F 1 is the S-C force experienced by the water. G 11 is the interaction strength acting within the water component and ψ 1 is computed using PR EOS . The strength of the interaction force acting between the component can be altered using G 12 .

The S-C force computed for air component F 2 is shown below:

F 2 (x) = -c 2 s G 12 ϕ 2 (x)∇ϕ 1 (x) (3.21)
The value of G 12 is maintained similar for computing both components. If the value of G 12 is positive or zero the diffusivity is very small and it is difficult to reach the equilibrium and if the value of G 12 is negative then there is an attractive force between the components which makes the air to stay in the liquid phase of water [130][21]. The value of G 12 is maintained negative to cause the attraction force for the air increasing the diffusivity and the value of G 12 is maintained positive in the water to cause the repulsion to separate the 2 components [START_REF] Chen | A critical review of the pseudopotential multiphase lattice Boltzmann model: Methods and applications[END_REF].

Results of implemented S-C method for SCMP and MCMP flows 3.3.3.1 SCSP : Poiseuille Flow (2D)

Single component single-phase (SCSP) simulation of a classic case of Poiseuille flow in 2 dimensions will be carried out to check the accuracy, stability and characteristics of the LBM simulation. Flow in a channel is modeled for 2-dimensional incompressible flow between two stationary parallel walls and the flow is driven by the velocity imposed at the inlet. The domain considered for the simulation has a length of L x = 10m and the width of the channel is L y = 1m, where unit length is discretized into 100 grid points. The maximum velocity u max = 1m/s and the Reynolds number Re = 10. The boundary condition considered for the simulation is: The analytical solution for a laminar flow in a channel is defined as :

u x = u max • 4(y -y 2 ) u y = 0    at the inlet (3.
       u x (x, y) = 4 (y -y 2 ) u y = 0 P (x, y) = -8(x -2)/ Re (3.25)
Here, u is the analytical solution in velocity and p is the analytical solution in pressure inside the domain.

The analytical solution is validated with the simulation performed using the LBM and the velocity profile obtained in the simulation overlaps with the analytical solution. The simulation is performed using the different resolution and produces a similar result with very slight error and the error decreases as the resolution increases. The simulation is performed for the resolution of 20, 40, 60, 80 and 100 grid points for unit length. The weak compressibility effect of LBM can be illustrated using the case of channel flow. The simulation is launched with an initial condition of zero velocity U 0 = 0 m/s and the pressure is considered as constant. The simulation is performed with grid size of N = 100 for unit length and the Reynolds number of Re = 10. Due to this effect, the initial conditions which are imposed are very important to obtain accurate results in LBM during the unsteady simulation and requires more time to converge to a stable solution in case of steady flow.

Stationary SCMP simulation

The initial step in modeling multiple fluids is to start with a single-component with multiphase, as it is easy to implement the developed equations for higher density ratio. In this case, we consider a single non-ideal component and to separate the phases by changing the critical density values. The phase separation gives rise to other effects at high-density ratios which must be studied as well. Optimization of the simulation helps in increasing the stability and the accuracy of SCMP which can be applied on to the MCMP model.

The simulation of SCMP is performed on 2 dimensions for a stationary case with periodic boundary condition on all the walls. A high-density liquid phase drop is suspended and the low-density gaseous phase will be distributed surrounding the drop. The total size of the domain considered for the simulation is 200 lattice units. The initial condition is applied in such a way that the liquid phase density of ρ liq is applied at the middle in a circle of radius 30 lattice units and gas-phase density ρ gas in all the other location of the lattice. The external forces such as gravitational forces are ignored in the simulation. The simulation is carried out for 100000 time steps to reach the study state. Shift velocity method is used for computation of external and inter-node interaction force and the relaxation time τ is considered as 1. The lattice structure considered for the simulation is D2Q9 and BGK collision model to simulate LBM equation. All the units (MLT) are expressed in Lattice units. LATTICE BOLTZMANN METHOD The density ratio of the phases can be controlled by changing the critical parameters of EOS. Here, in this simulation the reduction temperature T r will be varied to obtain different density ratio, where the reduction temperature is a function of the local temperature T and critical temperature T cr of that particular fluid (In this case it is water) T r = T /T cr . The value of T r can be reduced to obtain high-density ratio between the phases.

The density distribution at end of the simulation for a static case for T r = 0.8 is shown in the figure 3.7, where the high density is the liquid phase and lower density is the gaseous phase. The intermediate values of density observed in the interface formed between the 2 phases. The radius of the drop depends on the initial density distribution value of the simulation, as the density value at the liquid phase and gaseous phase changes depending on the critical parameters of the EOS. The initial value of density distribution is crucial as the higher the difference between the study state density values and initial values, which causes numerical instability. The simulation is carried out for the different value of T r to check the stable initial conditions and also to determine the density ratio obtained. The different initial condition is applied for the particular value of T r and the simulation is carried out. Certain values lead to numerical instability and the solution diverge. If the values are closer to the steady state stable condition, the simulation converges to an accurate value.

The results are plotted from the middle of the domain in complete x-axis which is shown in the plot below. From the plot, we can observe that, as the value of T r reduces the density ratio increases between the liquid and the vapor phase due to the EOS. The critical density value changes as the value of T r is changed which contributes to changing the density ratio. To initialize with the appropriate initial value, the simulation is run for different values for gas phase and liquid phase density for a few time steps for each T r value to determine the stable value. These stable values are given as the initial condition to simulate for different T r values. The T r value is varied from 0.95 to 0.6 and the respective density ratio at the study state is obtained. From the table 3 The density ratio is around 5 for the T r value of 0.95 and the density ratio increase by lowering the value of T r . The density ratio in the range of 1000 and above can be noticed for the T r value 0.65 and below. The maximum density ratio is around 8700 for the T r value of 0.6. If the value is reduced lower than the value of 0.6 the solution leads to numerical instability. This instability is caused due to the development of high-density ratio between the interface. As the density ratio between the interface increases leading to the development of high forces, which in turn leads to an increase of velocity. When the velocity of the flow reaches the speed of sound of the lattice c 2 s the solution becomes numerically unstable [START_REF] Gong | Numerical investigation of droplet motion and coalescence by an improved lattice Boltzmann model for phase transitions and multiphase flows[END_REF][10]. The velocity which is developed due to these forces is commonly known as spurious Velocity which will be discussed in the next section.

Spurious Velocity

Spurious velocities are non-physical velocity that is caused due to the fluid-fluid interaction force occurring at the interface of the fluids. When the spurious velocity approaches the speed of sound in lattice c 2 s , the simulation will become numerically unstable. The instability of the simulation in most of the S-C method is due to this spurious velocity. The amplitude of the spurious velocity increase as the density ratio increases which limits to obtain a high-density ratio.

In the case of S-C EOS, it reaches high spurious velocity due to the absence of attractive and repulsive terms compared to the other EOS at low-density ratio [START_REF] Stiles | High density ratio lattice Boltzmann method simulations of multicomponent multiphase transport of H2O in air[END_REF].

The Spurious velocities can be observed from the simulation performed at T r = 0.8 as shown in the figure 3.9, which is caused at the interface of the two-phase and gradually reduces away from the interface. The magnitude of the velocity distribution is symmetric across the interface and can be observed higher in the vapor phase than the liquid phase. These might be due to low density at the vapor phase causing an increase in the force, in turn, the velocity. The velocity profile across x-axis is plotted in the figure 3.10 for the simulations performed at different T r values. From the plots, we can observe that the velocity varies from 1 × 10 -5 to 0.4 for T r vales of 0.95 to 0.6 respectively. These plots show the value of T r reduces, the spurious velocity amplifies due to the increase in density ratio. It can be noticed that if the value of T r is reduced lower than 0.6 which can lead the velocity amplitude greater than 0.4 which is far more higher than the speed of sound in the lattice, that is is 1/3 for the D2Q9 lattice structure.

There are other issues such as, the spurious velocity cannot be distinguished from the flow velocity causing difficulties in understanding the accuracy of the simulation. These spurious velocities can also influence the flow condition leading to improper results. As the interface thickness and the spurious velocity are interrelated and are important parameters to determine the stability of the system and thermodynamic instability of density LATTICE BOLTZMANN METHOD is also caused due to these parameters [START_REF] Gong | Numerical investigation of droplet motion and coalescence by an improved lattice Boltzmann model for phase transitions and multiphase flows[END_REF]. The interface thickness the spurious velocity can be altered by introducing parameter β in EOS which will be discussed in the next section.

Modifying pressure Term

The spurious velocity is highly distinctive in the vapor phase than the liquid phase. These can be noticed in the analytical solution of effective mass with respect to the change in the volume and there is 2 observation that is discussed in the article [START_REF] Stiles | High density ratio lattice Boltzmann method simulations of multicomponent multiphase transport of H2O in air[END_REF].

• The variation of ψ at liquid volume V liq is grater that around gas phase volume V gas .

• Discontinuous in the volume of phase, when the volume is less than the volume of liquid phase V < V liq .

From these observations is clear that gas-phase density is not well defined by inter-node interaction force model but the liquid phase density is properly defined [START_REF] Stiles | High density ratio lattice Boltzmann method simulations of multicomponent multiphase transport of H2O in air[END_REF] [START_REF] Gong | Numerical investigation of droplet motion and coalescence by an improved lattice Boltzmann model for phase transitions and multiphase flows[END_REF]. There are different methods to improve the fluid-fluid force implementation and incorporation of improved fluid-fluid interaction force into the LBM equations [START_REF] Kupershtokh | On equations of state in a lattice Boltzmann method[END_REF].

A simple method is to replace the pressure term p with p * = β × p in the equation 3.11, where β is a scalar quantity as defined in [START_REF] Kupershtokh | Criterion of numerical instability of liquid state in LBE simulations[END_REF], [START_REF] Stiles | High density ratio lattice Boltzmann method simulations of multicomponent multiphase transport of H2O in air[END_REF].When β < 1.0, it decrease the pressure amplitude across the interface and stabilises the internode interaction force at the vapor phase. The β value doesn't have a noticeable effect on the liquid phase density and smaller the value of β, lower temperature range can be obtained. The gas-phase density at low reduction temperature and lower β value are not thermodynamically correct.

The value of β is changed for the case with T r value of 0.6 to determine the effect and stability of the simulation. The value of β is varied from 1 which is the initial value to the value of 0.75 and 0.5. The value of β was reduced even further to the value of 0.25 and 0.125, but the density values of the vapor phase become unstable leading to an error. These errors are because the density values keep reducing as the time step increments which leads to a very low fraction of vapor phase density causing the numerical instability. Here, the results are discussed for the simulation with the β value of 1, 0.75 and 0.5.

The effect of β on the vapor phase helps in reducing the density of the vapor phase leading to higher density ratio but keeping the liquid phase density constant. Initially, the density value for T r was around the range of 10 4 and as the value of β decreases, the density ratio increases to 10 5 and 10 7 for the beta values of 0.75 and 0.5 respectively. This modification of the value of beta will help in obtaining a high-density ratio for the same T r . LATTICE BOLTZMANN METHOD The effect of β on pressure and velocity at the interface can be studied. The pressure variation plot is focused at the left interface of the simulation result, which is shown in figure 3.12a to study the variation of the pressure across the interface. The pressure value across the interface reduces from 0.47 for β value of 1 to 0.35 for the β value of 0.75 and 0.22 for the β value of 0.5.

The spurious velocity also decreases as the β value is decreased which helps in simulating higher density ratios. The velocity magnitude variation plot is focused at the left interface of the simulation in figure 3.12 to study the variation of the velocity across the interface. It can also be noticed that the thickness of the interface increases as the value of β decreases. The increases in the interface thickness will help in suppressing the spurious velocity as the density variation across the interface is gradual compared to the other case of β = 1 where the interface is thin which causes the large variation of density across the interface.

Stationary MCMP simulation

The multi-component simulation is performed using the data obtained from the SCMP cases, which were computed before. The simulation is performed similar to SCMP but in this case, we try to simulate for 2 components with the water drop at the middle and the air surrounding it. The different phases of the 2 components have to be distributed with respective position. Each lattice point should have both the components present at the given location, failing to distribute both the components causes the numerical error. The simulation is performed for the lowest T r value of 6 and the β = 0.5 to obtain the highest possible density ratio.

To simulate this case, the liquid phase of water is initialized with a drop of radius 30 lattice units and the vapor phase is distributed in the other location of the domain. The air is distributed in an inverse way compared to water, where the actual density of air is placed in the surrounding domain and a fraction of density of air is distributed inside the drop. The total domain considered for the simulation is 200 lattice units with periodic boundary condition on all the sides. The simulation is computed for 100000-time steps to reach a stable condition. Velocity shift method is used to compute the inter-node interaction forces. The viscosity of both the fluid is assumed similar to the relaxation time τ = 1 for both components. D2Q9 lattice structure and BGK collision model is considered to compute the LBM equations.

The appropriate initial density distribution must be provided to obtain stable and accurate results. Variation in the initial condition can lead to numerical instability. The density distribution of water inside the drop is given as 8.8 and the density of water vapor around the drop is given as 0.001. The density for air inside the drop is provided as 0.001 and the density outside the drop is provided as 0.01. As the simulation progresses the density values changes leading to the stable condition.

The density distribution of the initial condition has to be greater than that of the stable condition as the compression wave is created at the initial time step causing variation in the density distribution leading to numerical instability. Different initial values of density distribution have to be tried to check the stable initial condition. LATTICE BOLTZMANN METHOD The results at the end of the computation show, the density distribution of each phase and each component. The figure 3.14 show the distribution of air, where the actual density of air can be seen around the drop and the lower fraction of air density is seen inside the drop, the in-between range is the interface. There are 2 small circular drops which can be noticed near the top boundary. These drops are the condensed liquid water from the vapor phase which was present outside the drop. This is caused due to variation of density due to the compression wave and higher value of density present at the vapor phase. The velocity distribution of air and water component in the domain at the end of the simulation is shown in figure 3.17. From the figures, we can notice that the velocity magnitude of air is around 0.44 and the highest value can be noticed at the interface due to the interaction forces between the components. The velocity variation can also be noticed at the top of the simulation in the form of small bubbles representing the condensation of the water present in the domain outside of the drop. The velocity distribution of water reaches very high velocity around the interfaces of the 2 components. The velocity at certain regions goes higher than 1 which can cause numerical instability. The main reason for high velocity might be because LATTICE BOLTZMANN METHOD of the negative density of water and the difference between the density of 2 components is huge introducing large forces which in-turn producing high velocities. The compression wave is generated at the initial time step and continues to resonate. Due to this reason, the initial values of density distribution have to be higher than the stable density distribution. This resonance initiates the variation of density distribution causing the water vapor present in the air to reach the critical density of liquid and condensation of the vapors into liquid. These formations of small liquid droplets can be observed around the big drop in figure 3.18. The formation of small drops leads to the formation of the small interface around it leading it to develop spurious velocity. As the small drops are closer it leads to an increase in the velocity even more, which can lead to the numerical instability. Proper initial conditions can get rid of these condensation effects but due to the compression wave, the accurate initial conditions cannot be implemented. These problems coexist hence, different methods have to be studied to reduce the effect of compression wave or methods to increase the numerical stability which may allow introducing a proper initial condition.

The density of the water reaches negative value outside the drop and the reason for this phenomenon is unknown for this moment. The negative density distribution at the beginning and the end of the simulation is shown in figure 3.19. A higher concentration of the negative density can be observed at the interface of the two components and the negative density decreases. The negative density may represent the concentration of the other component is higher to maintain the interaction force between the 2 components. The Interaction forces computed for the air component is considered as attractive force towards the water component. This attraction force draws a larger concentration of air towards the interface of the two components. Due to this higher concentration of air and to maintain the equilibrium the density of the water is drawn to a negative value. At the other location, the concentration is negative as well, but the negative value is lower representing the concentration of air is not as high as at the interface. The reason might be due to the initial density of air is higher than the actual value and the additional density cannot be neglected. The density of the air is initialized higher than the actual value to stabilize the simulation during the compression wave. The compression wave causes variation in the density distribution which may lead to instability due to improper forcing between the components if the actual density of air is provided. The density of both components is added to represent the actual density of the system. In the total density distribution, the negative density cannot be noticed and the distribution is linear across the domain. the density varies across the domain giving the average density ratio of the domain is around 1000. The conditions analysed in the static drop simulation of MCMP will be applied to a channel flow problem where one component pushes the other component, which was the initial goal of this project.

MCMP in a channel flow

The Channel flow simulation of MCMP is carried out for a 2-dimensional case similar to the case of Poiseuille flow analyzed before. In this case, both the components are placed next to each other. The flow is simulated between two stationary parallel walls and the flow is driven by the velocity imposed at the inlet. The domain considered for the simulation is 100 lattice units in width and 500 lattice units in length. In the domain, the water component is considered at an initial 1/8 of the channel and the remaining domain is filled with the air component. The inlet velocity condition is considered at the water-side of the domain and the outlet is considered at the end of the air component as shown in the figure 3.21. The basic idea of the simulation is to make the water push air simulating the simplified liquid piston without compression. The maximum inlet velocity is 0.01 Lu/Lt (lattice units/lattice time) for the water component and 0 for the air component. The boundary condition applied for the simulation are shown below: The simulation is conducted for higher T r value as the simulation gets highly unstable and proper initial condition has to be determined, as the initial values of the previous case of the static drop cannot be applied to this case of channel flow. As the value of T r reduces, it is difficult to determine the accurate value of the initial density. The stable condition for which the simulation could be carried out is for T r = 0.75.

u x = u max • 4(y -y 2 ) u y = 0
To increase the stability of the simulation for the chosen value of T r , the initial condition of the walls are changed to symmetry boundary condition and the inlet and outlet is changed as the wall conditions. On the sidewall due to separations of two components, the interaction forces will be high and the simulation can become unstable if the wall exists. This instability can be reduced at the walls by proper treatment of the interaction forces. The computation of interaction forces are ignored at all the wall and are only computed in the interior.

The simulation is run with the side walls as symmetry condition and wall condition at the 2 ends of the domain in place of inlet and outlet until the stable condition is reached or for the compression wave to disappear. The computation of the inter-node interaction forces are ignored at the wall but computed at the symmetric condition. This condition is computed for 10000-time steps and the boundary condition is changed to the actual channel flow condition of walls on the side and inlet and outlet conditions after this time step. The inter-node interaction at the wall are ignored and the simulation is computed until the water component reaches the end of the channel.

The density distribution of the symmetric wall condition can be seen from figure 3.22a. The component of air contains straight fringes of density variation all across the air dominant LATTICE BOLTZMANN METHOD region. These fringes are caused due to the initial compression wave separating the water and air component and due to the symmetry condition, the forces are applied in one direction. The simulation is carried out until the variation of the fringes reduces. Once the simulation reaches the stable condition, the inlet, outlet and wall conditions are applied. After the proper BC are imposed the variation of the densities can be observed in the figure 3.22b. Due to the actual BC, the forces are applied on all the sides and small perturbation causes the water present in the air to form small droplets which were fringes initially. The velocity variation of air and water component can be seen as the inlet and outlet conditions are applied in figure 3.23a and 3.23b respectively. As discussed early in the static MCMP section the velocity variation can be seen across all the interface of the small drops and the large separation of water-air, the flow can be affected due to this local velocity created. In the velocity profile plot 3.25a, we can see the first case density distribution in the figure 3.22. The velocity varies a lot and does not follow the profile because the condition has just been applied hence the velocity profile has not developed. The peaks on the plot represent the velocity fluctuation caused due to the interface of the condensed drops present in the channel which causes the large deviation. The second case (middle) and the third case (end) are shown in figure 3.25b and the respective density variation is figures 3.24. From the velocity plot of the second case follows the trend of the analytical solution but has a peak due to the small droplet of water near the result location and the profile has changed because of the previous droplets of water present in the flow. In the third case, the velocity profile is adapting well to the profile of the analytical solutions in most of the places with a slight error at the wall. There are 2 reasons, the first is the presence of small water droplet next to the wall causing the deflection of the velocity on the right side of the profile and the second LATTICE BOLTZMANN METHOD reason is the inter-node interaction forces are ignored at the walls leading to the error values at the wall.

Results analysis of Shan-Chen method

The multi-fluid simulation must be carried out to simulate the liquid piston (water-air interaction). There are different methods to simulate the multi-fluid system in LBM. Here, S-C method is chosen to simulate the multi-fluid system (MCMP system). The original S-C method is stable for low-density ratio of around 10, hence modification of the EOS is carried out to increase the density ratio of around 1000 to simulate the water-air system. To simulate MCMP, the simulation must be stable for SCMP.

Preliminary analysis is carried out on SCMP to determine the relation of reduction temperature T r and density ratio obtained. The highest density ratio of 8700 is obtained for T r value of 0.6. The density ratio cannot be increased higher as the increase of density ratio generates spurious velocity at the interface of two phases. As the magnitude of the spurious velocity gets higher than the speed of sound c 2 s , the simulation becomes unstable. From the analytical solution of effective mass ψ, it can be noticed that the density variation is grater at the gas phase causing instability for a slight change in the density distribution of the gas phase. These instabilities can be decreased by introducing β parameter in the pressure term. Reducing the value of β, increases the stability allowing to obtain higher density ratio. The β value of 1, 0.75 and 0.5 is implemented allowing to obtain the highest density ratio of around 10 7 for the beta value of 0.5. Decreasing the β value below certain value causes numerical instability leading the vapor density to decrease for each time step reaching 0 eventually. Decrements of the value of β will also reduce the spurious velocities and increases the thickness of the interface allowing the gradual change in the density at the interface which increases the stability.

The stable conditions obtained from the SCMP is applied on to MCMP simulation. The providing initial condition to the MCMP are complicated compared to the SCMP, as four different initial densities of two components must be defined. The compression wave in the LBM stimulation causes the variation of the density which changes the interaction force and change in these forces may lead to numerical instability. The variation of the density will also lead to pressure variation causing the vapor present in the air to condense to form water droplets. The interfaces are created around these drops causing the spurious velocity. Accumulation of these velocities may also lead to numerical instability. The value of the water vapor reaches negative value during the simulation, which causes large uncertainty . The reason is unknown for this instance, the investigation has to be carried out to determine the reason for the negative value of the vapor phase.

The simulation is carried out for the MCMP in a channel flow where, one fluid pushes the other fluid. The simulation is highly unstable, hence it was unable to simulate of T r value LATTICE BOLTZMANN METHOD below 0.75. The velocity profile of the simulation is in understanding with the analytical solution at the core of the flow in the absence of small water droplets. There is a variation of results at the walls, as the interaction forces are ignored at the walls. The presence of small droplets of water in the air causes deviation in the results because of the spurious velocity, which is formed at the interface influencing the flow velocity, in-turn introducing the error in the simulation. Developing accurate forcing scheme which increases the numerical stability can allow simulating for higher density ratios.

Free-Energy Phase Field models

Swift et al. introduced a free energy model [START_REF] Michael | Lattice Boltzmann Simulation of Nonideal Fluids[END_REF]. The initial idea is to re-derive the equation for the equilibrium density function using a free energy function. Following this method, the equilibrium function can be made to include a variety of effects of multiphase flows. For multiphase and multi-component flows, the free energy functional usually consists of three terms:

Ψ = V [ψ b + ψ g ]DV + A ψ s dA (3.29)
Where b, g and s are functions of space and time. The first term, ψ b , describes the bulk free energy. This term, serves primarily to define an equation of state (isothermal ideal gas, Landau free-energy models, Peng-Robinson EOS, van der Waals EOS) that allows for the coexistence of several fluid phases and/or components. While ψ g serves to capture the interface. Finally ψ s defines the interaction between the fluid and its surrounding solid (wetting phenomena).

Inamuro Model

The Inamuro model [START_REF] Inamuro | A lattice Boltzmann method for incompressible two-phase flows with large density differences[END_REF] belongs to the class of free-energy models developed initially by Swift et al. [START_REF] Michael | Lattice Boltzmann simulation of nonideal fluids[END_REF]. It is characterised by the use of two distribution functions, the model aims to solve the N-S equation and interface tracking governed by the Cahn-Hilliard equation. The first model developed by Inamuro in 1995 [START_REF] Inamuro | A non-slip boundary condition for lattice Boltzmann simulations[END_REF] was later upgraded in 2006 [START_REF] Inamuro | Lattice Boltzmann methods for viscous fluid flows and for two-phase fluid flows[END_REF] and finally in 2016 [START_REF] Inamuro | An improved lattice Boltzmann method for incompressible two-phase flows with large density differences[END_REF]. The main advantage of this model is its capacity to model two phase flow with high density ratios and at relatively high Reynolds number. Using a Lattice Kinetic Scheme to improve the LBM scheme, the proposed method can simulate incompressible multiphase flows with up to 1000 density ratio, which is the case of water.air flow in a liquid piston. Thus the method is implemented in a in-house Fortran code. The method solves two particle distribution functions f i , g i . f i is used to calculate the order parameter φ (Eq. 3.41), which distinguishes the two phases. g i is used to predict the velocity field of the flow without a pressure gradient. A third particle DF h i is introduced later to solve the Poisson equation (velocity field correction using pressure field). The D3Q15 scheme LATTICE BOLTZMANN METHOD is used in the model. The first D.F is the equation 8.9. A source term S i is added to LBE to obtain the second D.F 3.30.

g i (x + e i ∆t, t + ∆t) = g i (x, t) - [g i (x, t) -g eq i (x, t)] τ g + S i (3.30)
Where the source term S i models the effect of the viscous stress tensor (3.31):

S i = 3E i e iα 1 ρ (∂ β [η(∂ α u β + ∂ β u α ]) ∆t (3.31)
The order parameter φ is obtained through equations (8.9 and 3.32) and the predicted velocity u * is calculated through equations (3.30 and 3.33).

φ = 15 i=0 f i (3.32) u * = 15 i=0 g i e i ( 3.33)
The equilibrium distribution functions f i and g i are defined as [START_REF] Inamuro | An improved lattice Boltzmann method for incompressible two-phase flows with large density differences[END_REF] :

f eq i = H i φ + F i p 0 -κ f φ∇ 2 φ - κ f 6 |∇φ| 2 + 3E i φe iα u α + E i k f G αβ (φ)e iα e iβ (3.34)
and (3.35)

g eq i = E i 1 + 3e iα u α + 9 2 e iα u α e iβ u β - 3 2 u α u α + 3 2 (τ g -0.5)∆x ∂u β ∂x β + ∂u α ∂x β e iα e iβ
+ E i κ g ρ G αβ (ρ)e iα e iβ - 2 3 F i κ g ρ |∇ρ| 2 (3.36)
Where :

G αβ (ρ) = 9 2 ∂φ ∂x β ∂φ ∂x α - 3 2 ∂φ ∂x γ ∂φ ∂x γ αβ (3.37)
The coefficients,E i , H i , F i are defined in the D3Q15 scheme as :

E 0 = 2/9, E 1 = ... = E 6 = 1 9 , E 7 = ... = E 14 = 1 14 H 0 = 1, H 1 = ... = H 14 = 0 F 0 = -7 3 , F i = 3E i (i = 1, ..., 14)
In the above equations, δ αβ is the Kronecker delta, κ f is a constant parameter determining the width of the interface, g is a constant parameter determining the strength of the surface LATTICE BOLTZMANN METHOD tension, and the other variables, ρ , ρ L , µ , and u * are defined below : p 0 is calculated by the Van der Waals EOS :

p 0 = φ ∂ψ ∂φ -φ = φT 1 -bφ -aφ 2 (3.38)
Where ψ(φ, t) = φT ln φ 1-bφ -aφ 2 , a and b are re free parameters determining the maximum and minimum values of φ. The density in the interface is obtained by using the cut-off values of the order parameter, φ * L and φ * G , for the liquid and gas phases with the following equation 3.39:

ρ =            ρ G φ > φ * G ∆ρ/2[sin φ-φ * ∆φ * π + 1] φ * G ≤ φ ≤ φ * L ρ L φ > φ * L (3.39)
where ρ G and ρ L are densities of gas and liquid phases, respectively. ∆ρ

= ρ L -ρ G , ∆φ * = φ * L -φ * G and φ * = (φ * L -φ * G )/2.
The dynamic viscosity η is calculated using the given parameters :

η = ρ -ρ G ρ L -ρ G (η L -η G ) + η G (3.40)

Cahn-Hilliard equation

The Cahn-Hilliard (CH) equation is usually used to track the interface explicitly [START_REF] Jacqmin | Calculation of Two-Phase Navier-Stokes Flows Using Phase-Field Modeling[END_REF]. In the used models, it is used to recover the order parameter φ instead of using an equation of state such as in S-C models described in previous section. The forth degree partial equation can be written as :

∂φ ∂t + ∇.(φu) = ∇.(M ∇µ) (3.41)
The diffusive term on the right-hand side accounts for the motion of the order parameter due to inhomogeneities in the chemical potential. In many cases, the mobility parameter M is taken to be constant, although in general it is a function of the fluid order parameter. u = (u α , u β , u g amma) In the above derivation the following formula are used:

i f eq i = φ, i f (1) i = 0 and i f eq i e iβ = φu β (3.42)
Because the incompressibility condition is satisfied in this method, the equation 3.41 can be written as : [START_REF] Inamuro | Lattice Boltzmann methods for viscous fluid flows and for two-phase fluid flows[END_REF] ∂φ ∂t

+ ∂φu α ∂x α = ∂φ ∂t + u α ∂φ ∂x α = λ ∂ 2 P αβ ∂x α x β (3.43)

Poisson Equation

Since u * is not divergence free, a correction of u * is required. The current velocity u which satisfies the continuity equation (∇.u = 0) can be obtained by using :

Sh u -u * ∆t = -∇P ρ (3.44) ∇. -∇P ρ = Sh ∇.u * ∆t (3.45)
where Sh = U/c is the Strouhal number and p is the pressure. The third distribution function h i is used to solve the Poisson equation :

h n+1 i (x + e i ∆t, t + ∆t) = h n i (x, t) - [E i p n (x, t) -h n i (x, t)] τ h - 1 3 E i ∂u α * ∂x α ∆x (3.46)
where n is the number of iterations and the relaxation time τ h is given by :

τ h = 1 ρ + 1 2 (3.47)
The pressure is obtained by Eq. 3.48:

φ = 15 i=0 h i (3.48)
Numerically, during the LBE solving, the iteration of Eq. 3.46 is repeated until |p n+1 -p n |/ρ < (e.g., = 10 -6 ) is satisfied in the whole domain.

Algorithm process :

The model described above is computed as follows 3.

Test case : Two droplets collision in 3D

In order to validate the implemented algorithm, a set of comparison tests have been ran and compared to results from the original paper [START_REF] Inamuro | Lattice Boltzmann methods for viscous fluid flows and for two-phase fluid flows[END_REF][START_REF] Inamuro | An improved lattice Boltzmann method for incompressible two-phase flows with large density differences[END_REF]. The method was applied to the simulations of three dimensions binary droplet collisions with the same diameter D for various Weber numbers 20 < W e < 80 where W e = ρ L DV 2 σ and σ is the surface tension. Two liquid droplets having the same diameter D are placed 2D apart in a gas phase, and they collide with the relative velocity U (Figure 3.27). Several density ratios have been tested, while the one retained for validation is ρ L ρ G = 50. The coalescence collision phenomena is simulated (Table 3.2). The model is validated through comparison with theoretical results and numerical LATTICE BOLTZMANN METHOD Table 3.2: Droplet collision : simulation parameters results obtained in Inamuro's papers [START_REF] Inamuro | Lattice Boltzmann methods for viscous fluid flows and for two-phase fluid flows[END_REF][START_REF] Inamuro | An improved lattice Boltzmann method for incompressible two-phase flows with large density differences[END_REF]. Figure 3.29 shows the velocity field at t * = tU D = 1, before the impact, and Figures 3.28a and 3.28b show the density field at t * = 1.5 at the impact and show the two droplets coalesce at t * = 5.5. Then, the two droplets take a toroidal shape at t * = 6.5. This shape is later (at t * = 8) stretched due the kinetic energy. This observation is consistent with Inamuro's results and with the experimental results of Ashgriz and Poo [START_REF] Ashgriz | Coalescence and separation in binary collisions of liquid drops[END_REF]. The model showed high instabilities at higher density ratios ( ρ L ρ G > 120). This could be due to the choice made for the mobility λ (Eq. 3.49) which Inamuro et al [START_REF] Inamuro | A lattice Boltzmann method for incompressible two-phase flows with large density differences[END_REF] defined as :

λ = (τ f -1/2 -1/3C)∆x (3.49)
There is a balance to find in order to maintain a stable LBE (Eq 8.9) and a small value for mobility which reduces the flow's damp caused by the diffusion [START_REF] Jacqmin | Calculation of Two-Phase Navier-Stokes Flows Using Phase-Field Modeling[END_REF]. For our model, The added term C to reduce the mobility was chosen to be null C = 0. Inamuro recommended it to be very small in some cases and can be null on other [START_REF] Inamuro | Lattice Boltzmann methods for viscous fluid flows and for two-phase fluid flows[END_REF].

Results analysis of Inamuro's model

Inamuro's method was implemented in an three-dimensional in-house code based on Inamuro's first approach [START_REF] Inamuro | Lattice Boltzmann methods for viscous fluid flows and for two-phase fluid flows[END_REF], and then adapted to recent Inamuro et al. developments [START_REF] Inamuro | An improved lattice Boltzmann method for incompressible two-phase flows with large density differences[END_REF].

The obtained results enabled a comparison and validation of the model and the implemented code. Compared to S-C method, the free-energy model showed no spurious current for multicomponent flows even at relatively high density ratios (100) and high Reynolds number (2000). While it is mentioned by the authors that the method enabled simulations of multiphase flows with density ratios that can go up to 1000 in some cases [START_REF] Inamuro | Lattice Boltzmann methods for viscous fluid flows and for two-phase fluid flows[END_REF], the implemented code and the tested cases could not achieve a density ratio higher than 120. While, it remains a high one, it is still not sufficient for a water/air multiphase flow. Added to this, the initial method can only simulate incompressible flows, few modifications have been introduced to the method to enable two phase compressible flows simulation, but could not give appropriate results. Further developments of the mathematical model behind the method goes beyond this thesis objectives. Globally, using a Poisson equation to correct the velocity made the computational LATTICE BOLTZMANN METHOD load of the simulations much heavier, and did not enable to take any benefits from using the LBM.

Lee-Lin model

Lee and Lin model is based on the phase-field theory, it aims to reduce the parasitic currents in the LBM for multiphases simulations [START_REF] Lee | Eliminating parasitic currents in the lattice Boltzmann equation method for non ideal gases[END_REF][START_REF] Amaya | Single bubble rising dynamics for moderate Reynolds number using Lattice Boltzmann Method[END_REF]. They introduced several changes to He-Chen-Zhang model to increase the model's stability for high density ratios flows [START_REF] He | A Lattice Boltzmann Scheme for Incompressible Multiphase Flow and Its Application in Simulation of Rayleigh-Taylor Instability[END_REF]. The proposed model comprises the low Mach number approximation, the use of stress and potential forms of the surface tension force, the incompressible transformation, and the consistent discretization of the intermolecular forcing terms. It solves two particle D.F, one for the order parameter and the other one models the pressure and momentum. Lee and Lin [START_REF] Lee | Eliminating parasitic currents in the lattice Boltzmann equation method for non ideal gases[END_REF] suggested two forms for the two Lattice Boltzmann equations f i and g i : the potential form and the stress form. The stress and surface tension are introduced as an external force in Eq. 8.9 defined as F α . Where :

F α = ∂(ρc 2 s -p) ∂x β δ αβ + κ ∂ x β ∂ρ ∂x γ ∂ρ ∂x γ δ αβ - ∂ρ ∂x α ∂ρ ∂x β (3.50)
Note that the pressure P is modified as :

P = P -κρ ∂ 2 ρ ∂x γ ∂x γ + κ 2 ∂ρ ∂x γ ∂ρ ∂x γ (3.51)
The second term on the right-hand side of Eq. 3.50 is the surface tension stress tensor and its principle axes are directed in and perpendicular to the tangent plane of the interface which means that the surface tension will not change the total momentum at the interface. Added to this, the modified pressure term leads to smoother interface. So, the stress form is used in the LBE for pressure and momentum. The potential form of the LBE is used for the order parameter to recover the Cahn-Hilliard equation macroscopically, since it is more suitable for phase separation. The potential form of the intermolecular force F α is written as :

F α = ∂(ρc 2 s -p) ∂x α δ αβ + κρ ∂ x α ∂ 2 ρ ∂x β ∂x β (3.52)
In their model, instead of using the non-ideal gas equation as in He et al. [START_REF] He | A Lattice Boltzmann Scheme for Incompressible Multiphase Flow and Its Application in Simulation of Rayleigh-Taylor Instability[END_REF], they used the free-energy scheme.

E f (ρ) = β(ρ -ρ G ) 2 β(ρ -ρ L ) 2 (3.53)
Where β is a constant, and ρ G and ρ L are the densities of gas and liquid phases. In a plane interface at equilibrium, the density profile across the interface at equilibrium is :

ρ(z) = ρ G + ρ L 2 + ρ L -ρ G 2 tanh 2z W (3.54) LATTICE BOLTZMANN METHOD Where the interface width W is W = 4 ρ L -ρ G κ 2β
, we can obtain, then, the surface tension force σ :

σ = (ρ L -ρ G ) 3 6 2κβ (3.55)
The discritized Boltzmann Equation (DBE) for mass and momentum (Eq. 8.9) is transformed into the DBE for hydrodynamic pressure and momentum in order to eliminate acoustic waves from two-phase DBE in the low frequency limit and to improve stability. The transformation is done through a definition of a new particle distribution function :

g α = f α + ( p c 2 s -ρ)Γ α (0) (3.56)
where :

Γ i (u) = w i 1 + e iα u α c 2 s + e iα u α e iβ u β 2c 4 s - u α u α 2c 2 s (3.57)
The model is implemented in three steps: the pre-streaming collision step, the streaming step, and the post-streaming collision step. For the complete discretization schemes of the model, see [START_REF] Lee | Eliminating parasitic currents in the lattice Boltzmann equation method for non ideal gases[END_REF][START_REF] Lee | A lattice Boltzmann algorithm for calculation of the laminar jet diffusion flame[END_REF].

The density, the velocity, and the hydrodynamic pressure are calculated below after the streaming step :

ρ = i f α (3.58) ρu = i e α f eq α + ∆t 2 κ ∂ x α ∂ρ ∂x γ ∂ρ ∂x γ - ∂ x β ∂ρ ∂x α ∂ρ ∂x β (3.59) p = i g i + ∆t 2 u α ∂ρc 2 s ∂x α (3.60)

Test case : Stationary liquid droplet in a gas domain in 2D

A liquid droplet inside a gas domain with periodic boundary conditions (cf. Fig. 3.7) is simulated in order to compare and validate the 2D model. The computational domain is 100 × 100. The pressure is initially set to zero everywhere. and the density ratio is 1000. The relaxation parameters for liquid and gas are τ L = 0.53 and τ G = 065, respectively. The simulation parameters as detailed in table 3.3. The relaxation time τ 1 = τ 2 = τ is a function of position and is calculated as [START_REF] Amaya | Single bubble rising dynamics for moderate Reynolds number using Lattice Boltzmann Method[END_REF]:

τ (x) = φ(x)τ L ) + (1 -φ(x))τ G (3.61)
The convergence of the model is defined by the stabilization of the surface tension σ.

E σ σ num (t) -σ num (t -∆t) σ num (t) < 10 -5 (3.62) LATTICE BOLTZMANN METHOD Domain R W ρ L ρ G Re Time(s) σ τ G τ L 100 × 100 20 4
1 0.001 ∼ 0 2 6 × 10 -6 0.55 0.65 Table 3.3: Simulation parameters of the droplet case An small initial velocity is imposed which leads to a faster stabilization of the flow [START_REF] Lee | Eliminating parasitic currents in the lattice Boltzmann equation method for nonideal gases[END_REF]. The 2D liquid droplet immersed in the gas phases shows that the flow tends to be organized into eight eddies with centers lying on the interface 3.30 which is consistent with Lee and Fischer's results [START_REF] Lee | Eliminating parasitic currents in the lattice Boltzmann equation method for nonideal gases[END_REF]. The pressure and velocities in the stable mode show a 10% to 25% margin error in comparison to reference case (Table 3.4). This margin can be due to a small mobility M in the Cahn-Hilliard equation 3.41. The implemented model enabled a simulation of a multiphase flow with a high density ratio (1000). While a significant margin error was noted, the flow behavior was similar to the reference case one. The model can also handle small compressible flows. Nonetheless, simulating dynamic cases with a relatively low number Re > 10 was not possible with this model.

Variable

Velocity magnitude Pressure Simulation 7.9 × 10 -7 3.9 × 10 -7 Reference [START_REF] Lee | Eliminating parasitic currents in the lattice Boltzmann equation method for nonideal gases[END_REF] 3.9 × 10 -7 2.9 × 10 -7

Table 3.4: Comparison of the simulation results and the reference case's results

Results analysis of Lee-Lin's model

In comparison to previous models, Lee-Lin model enabled a simulation of multiphase flows with high density ratios. The model presents the advantages of handling the spurious current, LATTICE BOLTZMANN METHOD which did not affect directly the obtained results. Further improvement can be implemented to the model to increase its stability at relatively low Re number can be achieved, but as mentioned in the previous section, the choice was made to try another approach instead of developing a mathematical model that can improve this scheme. In addition, the final objective of the study is to simulate the flow inside a liquid piston which is a more complex one.

Volume-Of-Fluid LBM coupling

Volume Of Fluid (VOF) model constitutes an interesting alternative to pseudo-potential [START_REF] Shan | Lattice Boltzmann model for simulating flows with multiple phases and components[END_REF] color-gradient [START_REF] Gunstensen | Lattice Boltzmann model of immiscible fluids[END_REF] or phase-field [START_REF] Jacqmin | Calculation of Two-Phase Navier-Stokes Flows Using Phase-Field Modeling[END_REF] approaches for multiphase flows with high density ratios or high viscosity ratios. It was first proposed to be coupled to Navier-Stokes solvers.

Unlike the phase-field theory based on diffusive interface, VOF method is a sharp interface one. It is a model to simulate a sharp, freely moving interface between a fluid and gas phase in a Cartesian lattice [START_REF] Körner | Lattice Boltzmann model for free surface flow for modeling foaming[END_REF].

In the presented model, based on Thürey's work, the description of the liquid-gas interface is very similar to that of volume of fluid methods [START_REF] Thurey | Physically based animation of free surface flows with the lattice Boltzmann method[END_REF]. An additional variable, the volume fraction of fluid α, defined as the portion of the area of the cell filled with fluid, is assigned to each interface cell. The initial model proposed by Korner et al. is intended for free-surface flows simulation.

The free-surface model considers two-phase fluid and is particularly adapted for high density ratios. It simulates the physics of the heavy fluid only, and replaces the effect of the light fluid onto the heavy one through a condition of zero parallel shear stress along two-phase interface [START_REF] Körner | Lattice Boltzmann model for free surface flow for modeling foaming[END_REF] (Figure 3.32. This model is implemented in Palabos [START_REF] Latt | Palabos: Parallel Lattice Boltzmann Solver[END_REF] and validated in other studies [START_REF] Mohd | Lattice boltzmann method for free surface impacting on vertical cylinder: A comparison with experimental data[END_REF][START_REF] Thorimbert | Coupling of lattice Boltzmann shallow water model with lattice Boltzmann free-surface model[END_REF]. After identifying each cell (fluid α = 1, interface0 < α < 1, no fluid(gas)α = 0), the algorithm consists of three steps :

• the computation of the interface movement A modification to this method is proposed, it computes, in addition to the flow in the fluid cells and interface, the gas phase also through a second particle DF. It is therefore implemented in Palabos [START_REF] Latt | Palabos: Parallel Lattice Boltzmann Solver[END_REF].

Computation algorithm

The VOF-LBM coupling algorithm computes several steps to track the interface, model the two phases and ensure mass conservation. It is implemented as follows : Interface tracking : The movement of the fluid interface is tracked by the calculation of the mass in each cell. For this, two additional values need to be stored for each cell, the mass m and the fluid fraction α. The fluid fraction can be computed with the mass by dividing it through the density of the cell α = m/ρ. The interface motion if tracked by computing the the fluxes between the cells. However, as the DFs correspond to a certain number of particles, the change of mass can be directly computed from the values that are streamed between two adjacent cells for each of the directions in the model. For the interface and fluid cells, the mass is given :

∆m i (x, t + ∆t) = f I (x, t + ∆te i , t) -f i (x, t) (3.63)

Mass update :

It follows from the definition of the volume fraction that the mass m of a cell has the same units as the populations f i and g i of the lattice.

ρ L = i f i (x, t) (3.64) ρ G = i g i (x, t) (3.65) m(x, t + ∆t) = m(x, t) + i                    f I (x, t + ∆te i , t) -f i (x, t) if x is fluid f I (x, t + ∆te i , t) -f i (x, t)        1 if x is fluid α(x)+α(x+∆tei,t) 2 if x is interface 0 if x is gas if x is interface g I (x, t + ∆te i , t) -g i (x, t) if x is gas (3.66)
Where e I = -e i In order to ensure fluid mass conservation, there is no mass exchange between the gas cells and the interface cells. The amount of fluid leaving one cell has to enter the other one, and vice versa. The mass change values for all directions are added to the current mass for interface cells, resulting in the mass for the next time step.

Linking the two phases:

The populations from the gas cells cannot be streamed to the interface cells, since their value is calculated apart. In order to link the two equations (f i and g i , the interface cells population is completed as :

g out I (x = g eq i (ρ G , u(x)) + f eq j (ρ G , u(x -e i ∆t)) -f out j (x) (3.67)
Where g out I (x is the missing population coming from the gas cell, the superscript out denote the post-collision populations (before streaming). Eq . 3.67 preserves the off-equilibrium populations, as After this completion step is applied to all the neighboring gas cells of the interface, a normal streaming step can be applied to the whole domain, without any flagdependant treatment of the cells. An important point during this step is to properly handle mass excess; indeed, when a cell switches to gas cell, its mass is actually lower than α, whereas it is larger than 1 when it becomes totally filled. The mass excess in the case of a cell becoming fluid, and m g = m -ρ in the case of a cell becoming empty) must be redistributed somewhere else in the domain, in order to preserve the total mass of the fluid. Flag update : Based on the interface advection, at the end of the streaming-collision processes, exclusively for interface points, two interface change flags may be set, either interface ⇔ fluid or interface ⇔ gas:

m L (x, t + ∆t) =      > (1 + η)ρ(x, t + ∆t) ⇒ set interface → fluid < (η -0)ρ(x, t + ∆t) ⇒ set interface → gas (3.68) LATTICE BOLTZMANN METHOD
Where η ∼ 10 -4 is a small number with the purpose to increase the inertia of flag update process and thus, prevent from fast switching between flags.

Surface tension calculation

Two modes of surface tension calculation have been tested, one with constant gas density calculated from the pressure [START_REF] Körner | Lattice Boltzmann model for free surface flow for modeling foaming[END_REF]:

p = c 2 s ρ G -2κσ (3.69)
Where ρ G is here a constant and κ is the local mean curvature. The other model, called "dynamic", calculated the pressure and thus, the surface tension based on the numerical results ρ G according to equation 3.65.

LES Smagorisky turbulence model

In order to capture turbulent structures in the flow, a large eddy model is used. In LES models, a spatial filter is applied to the velocity field. The filter should be fine enough that the large turbulent structures of the flow do not get filtered out. Only the effect of the small eddies on the large-scale flow structures needs to be modeled and is included in the simulation through an additional turbulent viscosity ν T . A dynamic Smagorinsky model in which ν T depends on the shear rate is used [START_REF] Janßen | A lattice Boltzmann approach for freesurface-flow simulations on non-uniform block-structured grids[END_REF]:

ν T = (C s ∆x) 2 ||S|| (3.70)
Where C s is the Smagorinsky constant and S is the stress tensor

Test cases : compressible and incompressible two-phase Poiseuille flow

The initial free-surface model implemented in Palabos is already validated in a previous study [START_REF] Mohd | Lattice boltzmann method for free surface impacting on vertical cylinder: A comparison with experimental data[END_REF]. The new implemented model, which can simulate two-phase flows with high density ratios, is tested through three cases, the first is a flow with two phases where the dense fluid (water) pushes and expels the light fluid (air) from the top of the geometry is an outlet surface, for the second and third ones, the top of the rectangular cuboid is closed i.e compression (Figure 3.33). The domain is initially filled with air (light phase) and water is injected from the bottom at constant speed. The proposed test cases required the development of a Zou-He inlet boundary condition [START_REF] Zou | On pressure and velocity boundary conditions for the lattice Boltzmann BGK model[END_REF], with specifying the phase that will be injected. To better handle the small relaxation time, a regularized LBM [START_REF] Latt | Hydrodynamic limit of lattice Boltzmann equations[END_REF] is used. The parameters of the tested cases are presented in the table 3. Table 3.5: Parameters of the three presented cases Figure 3.34a shows the profile of the Z component (U z ) of the velocity in the studied case of a flow without compression (Direction of flow from left to right). We can see that the parabolic shape of the flow is well established. While figure 3.34b shows the volume fraction α showing the interface Liquid/gas advancement.

While no similar case is available in the scientific literature, the water domain should be similar to a Poiseuille flow. Thus, in order to compare and validate, the velocity profile is presented 3.35. Both liquid and gas velocity profile show a parabolic profile. The results were validated according to the analytical solution presented in Eq 3. [START_REF] Suraj | Evaluating the performance of the two-phase flow solver interFoam[END_REF].

Other tests are carried out to simulate a gas compression using a liquid piston. Liquid (Water) is injected at constant velocity into a closed chamber 3.33 and compresses the gas Results show a loss oscillating between 1.5% and 5% of the total liquid mass. While

Conclusion

The main challenge through a liquid piston simulation using the LBM is the handling of the multiphase high density flow. The original methods that enable a multiphase flow simulation are adapted for a small density ratios [START_REF] Gunstensen | Lattice Boltzmann model of immiscible fluids[END_REF][START_REF] Shan | Lattice Boltzmann model for simulating flows with multiple phases and components[END_REF]. Several improved methods are available in the literature and three different approaches were implemented and tested (S-C with new EOS, free-energy models and finally a new challenging method which couples LBM to VOF). The methods were implemented in either Palabos open-source code (VOF, S-C) or in an LATTICE BOLTZMANN METHOD in-house code (Lee-Lin, Inamuro). S-C method's interface tracking was enhanced through the use of more complex EOS in the pseudo-potential forces, instead of the original ideal gas one. SCMP and MCMP simulations were realized. The method suffers from spurious currents especially for the MCMP flows resulting in unstable schemes. The Free-energy models implemented, enabled simulations of multiphase flows with high density ratios (up to 1000). Lee-Lin method was tested for a static droplet at 1000 density ratio but could not enable a model a dynamic case, while Inamuro's method, which has a very high computational cost due to the Poisson equation for pressure correction, enabled a simulation of two droplet collision with high Re. Finally, a novel VOF method for two-phase flow is tested. It allowed simulation of high density ratio/high viscosity ratios multiphase flows. The method suffers from lack of mass conservation. Table 3.6 sums the tested method, their advantages and main remaining challenges. All the tested method did not give us the expected results to continue towards a liquid piston flow modeling. A choice had to be taken between further the free-energy approach which can require a new mathematical model to improve its Galilean invariance, modeling the liquid/gas interface as a moving solid boundary or leaving aside the LBM and switching to the more advanced and popular Navier-Stokes approach. Further development of a new model goes beyond this thesis scopes and its context (REMORA project), in addition, the complex flow inside a liquid piston during compression/expansion process requires other challenging improvements (Turbulence, compression, non-isothermal simulation). Tests have been carried to model the solid piston using Immersed Boundary Method [START_REF] Peng | A comparative study of immersed boundary method and interpolated bounce-back scheme for no-slip boundary treatment in the lattice Boltzmann method: Part II, turbulent flows[END_REF], this approach suffered from leakage between the moving interface and the walls. Another approach based on a moving boundary has been investigated based on Bouzidi's and Lallemand's methods [START_REF] Bouzidi | Momentum transfer of a Boltzmann-lattice fluid with boundaries[END_REF][START_REF] Lallemand | Lattice Boltzmann method for moving boundaries[END_REF]. Interpolation problems have been faced during the interface advancement. Thus, due to the complexity of the flow inside a liquid piston and the context of the thesis, a Chapter 4

Modeling of the compression process using a liquid piston 4.1 Introduction

The energy storage stage in a CAES system using a liquid piston consists of several sequence processes. Figure 4.1 shows the successive stages of compression in the liquid piston. The chamber is initially empty. A hydraulic pump compresses the air to a pressure P c by injecting water at a certain flow rate. When the air reaches the storage pressure, the same pump then pushes the compressed air from the chamber to the storage tank. When the piston is full of water, a non-productive draining phase allows to reset the compression chamber to allow a new compression by opening the top valve. Note that several compression chambers can be used in parallel to enhance the performance of the system [START_REF] Maisonnave | Optimal energy management of an underwater compressed air energy storage station using pumping systems[END_REF].

Our flow and heat transfer modeling study, presented in this chapter, only concerns the first process which consists of compressing the air by injecting the water from the bottom of the piston.

Air compression

Air ejection to tank Water ejection Time

Flow and heat transfer in a liquid piston

Gouda et al. made a summary of several numerical and experimental researches that studied on the physics of the liquid piston (heat transfer, efficiency, mass transfer), while very few focused on the flow inside the liquid piston [START_REF] Mehdi | Review on Liquid Piston technology for compressed air energy storage[END_REF]. Numerical studies have shown the positive effect that have liquid piston on maximizing heat transfer [START_REF] Zhang | Design analysis of a liquid-piston compression chamber with application to compressed air energy storage[END_REF][START_REF] Mutlu | Effects of piston speed, compression ratio and cylinder geometry on system performance of a liquid piston[END_REF][START_REF] Schober | CFD-simulation and experimental validation of heat transfer in liquid piston compressors[END_REF][START_REF] Park | Analysis and Proof-of-Concept Experiment of Liquid-Piston Compression for Ocean Compressed Air Energy Storage (Ocaes) System[END_REF][START_REF] Farzad | Optimal Control Experimentation of Compression Trajectories for a Liquid Piston Air Compressor[END_REF][START_REF] Zhang | Numerical investigation of metal-foam filled liquid piston compressor using a two-energy equation formulation based on experimentally validated models[END_REF][START_REF] Zhang | Numerical analysis of heat exchangers used in a liquid piston compressor using a one-dimensional model with an embedded two-dimensional submodel[END_REF].

It is proven through a 2D axisymmetric Computational Fluid Dynamics (CFD) study that higher Reynolds number results in lower heat transfer during the compression process and thus lower efficiency [START_REF] Zhang | Heat transfer in a long, thin tube section of an air compressor: An empirical correlation from CFD and a thermodynamic modeling[END_REF]. Similarly, a further study showed that the compression efficiency decreases when the velocity of the piston increases, when the compression ratio is high, and especially when the L/D ratio of the chamber is large [START_REF] Mutlu | Effects of piston speed, compression ratio and cylinder geometry on system performance of a liquid piston[END_REF]. A geometry optimization study concluded that a gourd shaped geometry favors the heat transfer inside the liquid piston during compression process [START_REF] Zhang | Design analysis of a liquid-piston compression chamber with application to compressed air energy storage[END_REF]. The effect of the compression time was studied and showed experimentally and through a numerical model that a longer compression time results in a lower compressed gas temperature [START_REF] Park | Analysis and Proof-of-Concept Experiment of Liquid-Piston Compression for Ocean Compressed Air Energy Storage (Ocaes) System[END_REF]. One can refer to chapter 1 where all the numerical and experimental studies on liquid piston are listed.

Main objectives and contributions

The aim of this numerical simulation is to allow a better understanding and to give a physical background to the flow and heat transfer evolution during the compression process inside a Liquid Piston. To achieve this goal, the experimental results allowed us to validate the numerical model through a comparison analysis. Numerical model also provides the results PISTON that help us to understand the origin and evolution of the instabilities that shape the flow and also the timing and progress of the transition from laminar to turbulent. It also offers better insights into the influence of flow dynamics on heat transfer, as well as the temperature field distribution of the air domain and finally the heat exchanges between air and water, and air and piston walls.

A 3D numerical study of air flow during compression within a liquid piston has not yet been documented in the literature. It is indeed a particular and complex flow since it is a compressible multiphase one with high density ratios generating a heat transfer. A transition from laminar to turbulent low regime can also be observed as well as strong local velocity gradients.

Numerical methods for multiphase flows

Two-phase flows are of great interest in natural and industrial applications such as rain formation [START_REF] Shaw | Particle-Turbulence Interactions in Atmospheric Clouds[END_REF], breaking waves [START_REF] Pizzo | Current generation by deep-water breaking waves[END_REF], spray-atomization [START_REF] Sirignano | Fluid Dynamics and Transport of Droplets and Sprays[END_REF] and bubbly flows [START_REF] Loisy | Buoyancy-driven bubbly flows: ordered and free rise at small and intermediate volume fraction[END_REF]. The purpose of this chapter is to model numerically the flow and heat transfer inside a liquid piston. As described in Chapter 1, this flow is a compressible and multiphase one. There are generally two different topologies for multiphase flows : dispersed flows and separate flows. The first means flows where the size of the particles of one phase (dispersed phase) is very small compared to the length scale characteristic of the flow of the carrier fluid, such as for example air bubbles in water. Separate flows are those which consist of two separate phases (stratified or pocketed) by a deformable interface. The main challenge in performing multiphase flow simulations is the ability of the numerical method to handle change in fluid properties across the interface. Many methods and approaches have been developed to solve the numerical phenomena. Mirjalili made a review of the most used methods, where they are summarized in the Figure 4.2 [START_REF] Mirjalili | Interface-capturing methods for two-phase flows : An overview and recent developments[END_REF].

VOF methods, which are the most used and robust ones, come in two variants:

• Geometric VOF methods; It uses geometric operations to reconstruct the fluid interface inside a cell and to approximate the water fluxes across faces.

• Algebraic VOF methods; It relies on the limiter concept to blend first and higher order schemes in order to retain sharpness and boundedness of the time advanced VOF field.

Geometric VOF schemes are typically much more accurate, but also computationally more expensive, complex to implement, and restricted to certain types of computational meshes, such as hexahedral meshes. Algebraic VOF schemes, on the other hand, are less accurate, but often faster, easier to implement, and developed for general mesh types. A combination of these two schemes is present in OpenFOAM [START_REF] Lin | Communication Optimization for Multiphase Flow Solver in the Library of OpenFOAM[END_REF], which is an open source library written in object-oriented C++ for Computational Fluid Dynamics (CFD) and is extensively applied in a wide range of engineering fields, such as aerospace, shipbuilding, Methods that use a sharp-interface approach have a white background while diffuse-interface approaches have a gray background [START_REF] Mirjalili | Interface-capturing methods for two-phase flows : An overview and recent developments[END_REF] oil exploration [START_REF]OpenFOAM | Free CFD Software | The OpenFOAM Foundation[END_REF]. In OpenFOAM, the interFoam solver family is dedicated to flows, where the pressure-velocity solver is coupled with the advection of a sharp fluid interface. These solvers apply the OpenFOAM specific algebraic Volume-of-Fluid (VOF) scheme called MULES (Multidimensional Universal Limiter with Explicit Solution) for the task of advecting the sharp interface. The MULES scheme has a number of desirable properties [START_REF] Suraj | Evaluating the performance of the two-phase flow solver interFoam[END_REF] :

• It is a mass conserving scheme, which preserves the volume of fluid, i.e. it does not artificially create or destroy fluid.

• It keeps the volume fraction field in the physically meaningful range between 0 and 1

Numerical model for multiphase flow and heat transfer in a liquid piston

In the liquid piston, the water/air interface is flat throughout the compression, and no bubbles or droplets appear. The water/air interface that acts as a piston is the only interface in the system.

For our numerical study, Finite Volume Method is used to solve the Navier-Stokes equations and is coupled to VOF Method for interface tracking, which is the most suitable one for this type of flows. Firstly, we define the volume fraction of the water component, which is given by Eq. (4.1) :

α = V water V cell (4.1)
Where V water is the volume of water in a cell of the computational grid V cell . The mixture density ρ and the viscosity µ is given by Eq. 4.2 and Eq. 4.3, respectively :

ρ = αρ water + (1 -α)ρ air (4.2) µ = αµ water + (1 -α)µ air (4.
3)

The compressible Navier-Stokes equations are written as :

• The mass conservation equation is given by:

∂ρ ∂t + ∇ • (ρU) = 0 (4.4)
where U represents the velocity vector.

• The momentum conservation equation is given by :

∂ρU ∂t + ∇ • (ρUU) -∇ • (µ∇U) = σκ∇α -g • x∇ρ -∇P d (4.5)
where σ is the surface tension coefficient; the curvature of the interface is calculated as κ = ∇ • ∇α |∇α| ; P d = P -ρg • x is the dynamic pressure; g is the gravitational acceleration, withg = gz and x is the position vector.

• The energy equation is given by: is the specific kinetic energy; c v,water and c v,air are the specific heat capacity at constant volume for the water and air phase, respectively. A transport equation for the water volume fraction used to capture the free surface is given by: ∂α ∂t

∂ρT ∂t +∇•(ρUT )-∇•(µ∇T ) = - α c v,water + 1 -α c v,air ∂ρe k ∂t + ∇ • (ρUe k ) -∇ • (UP )
+ U • ∇α + ∇ • U r α(1 -α) = 0 (4.7)
where ∇ • U r α(1 -α) is an anti-diffusion term utilized to sharpen the interface and U r = U water -U air is the relative velocity's vector, designated as the 'compression velocity'.

The time step is adjusted so the Courant Number does not exceed Co = 1 which is a necessary for the stability of the pressure-velocity algorithm PISO (Pressure-Implicit with Splitting of Operators) [START_REF] Issakhov | Numerical simulation of the movement of water surface of dam break flow by VOF methods for various obstacles[END_REF]. The numerical model is implemented through the open source code OpenFOAM using the CompressibleInterFoam solver [START_REF]OpenFOAM | Free CFD Software | The OpenFOAM Foundation[END_REF].

Thermodynamic Model

Despite the near-isothermal objective of the transformation, non-negligible local variations in air temperature appear. As this can affect its thermodynamic properties, Suhterland equation 4.9 is used to model the viscosity's evolution during the compression process. The transport equations of mass, momentum and energy are closed by the equations of state for density and enthalpy. Since temperatures don't exceed 400K and pressure are 5 × P atm , the air can still be considered as an ideal gas. Van de Ven & Li [START_REF] Van De Ven | Liquid piston gas compression[END_REF], Zhang et al. [START_REF] Zhang | Numerical investigation of metal-foam filled liquid piston compressor using a two-energy equation formulation based on experimentally validated models[END_REF], Mutlu & Kiliç, [START_REF] Mutlu | Effects of piston speed, compression ratio and cylinder geometry on system performance of a liquid piston[END_REF], Gouda et al. [START_REF] Mehdi | Méthode VOF pour la simulation numérique de l'écoulement de l'air comprimé et du transfert thermique associé dans un piston liquide[END_REF], used the ideal gas law to model the air. This model is chosen for this work. It is expressed by the following equation of state, which relates density to temperature and pressure:

ρ = 1 RT P (4.8)
Viscosity is a temperature related variable following the Sutherland equation :

µ(T ) = µ 0 T 0 + C T + C T T 0 3 2 (4.9)
Where (Table 4 The values of the initial and final average viscosities of air according to the equation 4.9 are listed in Table 4.2. An increase of 25% can not be considered as negligible since the term of viscosity intervenes in the resolution both in the momentum equation and in the energy equation. The air's heat capacity has a negligible relative difference, so this property can be considered constant for air during the simulation. Mass transfer between water and air domain is neglected.

As observed in previous work [START_REF] Mehdi | Méthode VOF pour la simulation numérique de l'écoulement de l'air comprimé et du transfert thermique associé dans un piston liquide[END_REF][START_REF] Mehdi | 2D versus 3D numerical modeling of flow and heat transfer in a liquid piston air compressor[END_REF], water maintains a constant temperature throughout compression in the liquid piston, and its pressure increases similarly to the air's one. The initial and final states are therefore respectively (T 0 = 293K, P 0 = 10 5 P a) and (T f = 293K, P f ∼ 5 × 10 5 P a). The relative difference of water density ρ water between the initial and final state shows that the flow inside the water column can be considered as incompressible since this difference is negligible. With the same reasoning, water's heat capacity c v,water , kinematic viscosity µ water , and Prandtl number P r = µwatercv,water λwater are modeled as constant throughout compression. OpenFOAM code (version 7) with the CompressibleInterFoam solver [START_REF]OpenFOAM | Free CFD Software | The OpenFOAM Foundation[END_REF] was used to solve the governing equations. A fixed value of water velocity at 0.033 m/s was imposed at the inlet boundary while all the other boundaries were treated as isothermal no-slip walls (Fig. 4.4). Detailed initial and boundary conditions for the simulation are listed in Table 4.3.

The selected discretization schemes and other numerical parameters can be found in Appendix 8.0.0.7. PISTON

Turbulence model for liquid piston compression processes

For the calculation of the Reynolds Number of the flow inside the liquid piston, one should distinguish the Re air based on air properties and Re water based on water column properties (Table 4.5).It is the velocity of the liquid piston is used as the characteristic velocity [START_REF] Van De Ven | Liquid piston gas compression[END_REF]. While Reynolds numbers are relatively small, Gouda et al. [START_REF] Mehdi | Méthode VOF pour la simulation numérique de l'écoulement de l'air comprimé et du transfert thermique associé dans un piston liquide[END_REF] showed the importance of choosing a turbulence model to simulate the liquid piston compression flow. PISTON The LES model showed that it is the most suitable one for catching and modeling the flow through the different regimes it passes through [START_REF] Mehdi | Méthode VOF pour la simulation numérique de l'écoulement de l'air comprimé et du transfert thermique associé dans un piston liquide[END_REF].

WALE and Smagorinsky sub-grid scale for LES turbulence model

The Large Eddy Simulation (LES), can offer a solution to balance the high computational cost of DNS (Direct Numerical Simulation) and the low accuracy of RANS (Reynolds-averaged Navier-Stokes) [START_REF] Germano | A dynamic subgrid-scale eddy viscosity model[END_REF]. LES spatially filters the Navier Stokes Equations and directly resolves the governing equations for large eddies (larger than filter size (generally mesh scale) and their turbulent energy, while it uses a Sub-grid Scale (SGS) model to simulate the effect of small eddies (smaller than filter size) on the mean flow and large eddies [START_REF] Lodato | A compressible wall-adapting similarity mixed model for large-eddy simulation of the impinging round jet[END_REF]. Computing resource requirement of LES is greater than RANS but is much less than DNS [START_REF] Lodato | A compressible wall-adapting similarity mixed model for large-eddy simulation of the impinging round jet[END_REF].

In this study, two SGS models in LES are discussed and employed in the LES analysis, the Smagorinsky model,and the Wall-adapting local viscosity model (WALE). Both of the two SGS models use the eddy-viscosity assumption, which try to simulate the diffusive transport and dissipation of kinetic and magnetic energy. The description of the general theory of the LES model and its filtering can be found in [START_REF] Leonard | Energy Cascade in Large-Eddy Simulations of Turbulent Fluid Flows[END_REF]. The two studied SGS models have a common definition of the subgrid-scale stress τ ij , which is defined as :

τ ij = 1 3 δ ij τ kk -2ν t S ij (4.14)
where δ ij if i = j and zero otherwise, ν t is the sub-grid scale eddy viscosity, S ij = 1 2 ∂u i ∂x j + ∂u j ∂x i is the deformation tensor of the resolved field and over-bar denotes the resolved quantities.

Smagorinsky model

Smagorinsky eddy-viscosity SGS model is the first employed in LES, and is still one of the most used SGS [START_REF] Smagorinsky | General Circulation Experiments With the Primitive Equations: I. the Basic Experiment[END_REF]. An assumption on τ ij can introduce the SGS kinetic energy k t = 2τ kk to Eq. 4.14, which can now be written :

τ ij = 2 3 k t δ ij τ kk -2ν t S ij (4.15)
The SGS eddy-viscosity in this model is assumed to be proportional to the sub-grid characteristic length :

ν t = (C s ∆) 2 |S| (4.16)
Where C s is a constant called the Smagorinsky coefficient, |S| = 2S ij S ij , ∆ is the grid size that defines the sub-grid length scale (grid spacing).

Wall-adapting local viscosity model

In this model, the SGS viscosity is evaluated using the square of the velocity gradient tensor; it is also negligible in laminar flow regions and in the viscous sub-layer and it has proper near wall behavior [START_REF] Nicoud | Subgrid-Scale Stress Modelling Based on the Square of the Velocity Gradient Tensor[END_REF]. It should be noted that the WALE model has the distinct advantage that it automatically recovers proper scaling in the near-wall region. ν t is defined as :

ν t = L 2 s (S d ij S d ij ) 3/2 (S ij S ij ) 5/2 + (S d ij S d ij ) 5/4
(4.17)

Where S d is the trace-less symmetric part of the square of the velocity gradient tensor and L s = C w ∆ a characteristic length depending on the wall-treatment model the near-wall mesh size.

Both the LES SGS models require a fine mesh near the wall (dimensionless wall distance Y + ∼ 1). This criteria has been respected for all the used meshes.

Through the comparison of average air temperature, the difference between the different models is not clear. Other comparison criteria were tested in the following sections. A separate study was conducted to choose the better LES sub-grid model between the Smagorinksy and WALE (Wall Adapting Local Eddy-viscosity). The WALE model is the most appropriate SGS model for this study, since it handles well the transition and also is more accurate in treating the near-wall region. (See appendix 8.0.0.7). A logarithmic law of the wall is chosen to model the boundary layer.

Results of numerical modeling of flow and heat transfer during the compression process

2D modeling was carried out in order to study the influence of certain simulation parameters, the advantage of 2D being the significantly shorter computational time than that of 3D. In addition, results of these simulations show the differences between the flow modeled in 2D and the flow modeled in 3D. The 3D modeling results were used for comparison with experimental results and for further flow analysis.

Mesh convergence study

A grid independence study was firstly conducted to determine the appropriate mesh size and structure for the CFD simulations for both the 2D and 3D models. Mesh comparison was based on the air's average air evolution and for the temperature at the end of compression when the temperature gradient was high. The meshes and the relevant spacing have been chosen respecting three requirements :

• Maintain a low-pass filter cutoff length in order to perform correctly resolved LES. PISTON

• Ensure proper resolution of the wall layer while keeping grid stretching to a minimum (Figure 4.6).

• Retain a reasonable computational cost due to the long computational times.

Mesh convergence study for the 2D model

The number of elements used for each study is presented in the table 4.6. All the 2D meshes are structured meshes throughout the domain, with refinement near the walls for the 400k mesh. The average air temperature curves resulting from the tests carried out during the mesh study are visible in Figure 4.5. Between 0 and 1s, the three resulting temperature curves are almost identical to each other. After 1s and until the end, the gap grows between the 60k curve and the other two. The final temperatures of the 240k mesh and the 350k mesh are 410K and 408K respectively, against 420K for the 60k one. The final temperature decreases as the number of meshes increases. The mesh sensitivity study, based on the average air temperature, therefore reveals that the mesh converges from 240k since the temperature no longer changes significantly.

Type 2D Rectangular

Number of elements ×10 3 50 180 400

Table 4.6: 2D meshes used in the mesh convergence study

Mesh convergence study for the 3D model

The conducted mesh convergence study included two types of mesh (Hexahedral and Prisms), different number of elements and different distribution of the meshing 4.7. To compare all the used meshes, two criteria were considered : the average temperature of the air and local velocities.

Grid independence study is conducted first to determine the correct mesh size and the mesh distribution in the numerical simulations. Two types of meshes are used (Prisms, Hexahedral). Different parameters are used to determine the mesh convergence : average air temperature, local velocity at defined points. For the hexahedral mesh, a structured O-Grid mesh is used, with mesh refinement at the wall. The hexahedral mesh shows globally a faster computational time for approximately the same number of element, it also allows a better local mesh refinement near the walls which enhances the accuracy of the solved model. At the end of the compression, where the temperature gradient are high and at mid-compression time where the highest velocities are PISTON Table 4.7: Types of 3D meshes used in the mesh convergence study observed, the hexahedral meshes shows. They also allowed a faster computing time. The Hexahedral mesh with 9 × 10 5 element was chosen for further analysis and for the rest of the study.

2D vs 3D comparison

The interest of 3D modeling is highlighted by comparing the velocity and temperature fields in 2D and 3D. The flow does indeed show significant differences between the two cases. The results of a simulation performed in 2D with a 400k mesh square structured and 3D 900k hexahedral meshes are used to explain these differences.

During the first phase, the 2D flow is identical to the 3D, except that the velocity extremes and the average temperature do not have the same values(Figure 4.8). After 5s, it is the area of air near the piston that deforms first, unlike in 3D (Figure 4.9). It is just after the moment of the regime transition that the differences between 2D and 3D models are the greatest. In 2D, vortices form and position themselves regularly on either side of the central current The benefit of switching to 3D modeling is strongly visible when comparing the air's average temperature curves (Figure 4.11). A comparison between numerical and experimental results [START_REF] Neu | Experimental investigation of internal air flow during slow piston compression into isothermal compressed air energy storage[END_REF] shows that the final temperature of the 2D simulation is 364K against 327K for the 3D. The difference with the final experimental temperature therefore goes from 41K to less than 2K. In fact, the 2D model assumes that the flow and the associated heat exchanges do not depend on the direction normal to the simulation plane and that the problem is equivalent to a compression performed between two infinite planes. Indeed, the amount of heat in the air exchanged with the outside depends on the surface area of the walls. However, for the same volume of compressed air, the wall surface area available for heat exchanges in Cartesian modeling is different from the experimental case. If we name C 2D and C 3D the compactness ratios of 2D and 3D cases, we have : The compactness ratio of the 2D model is thus twice less than that of the experimental case. Therefore, less heat will be exchanged with the external media (walls) [START_REF] Ramesh | Classification of Heat Exchangers[END_REF].

C 2D = V 2D S 2D = DL 2L = D 2 (4.18) C 3D = V 3D S 3D = π D 2 2 L 2π D 2 L = D 4 (4.19)
Another alternative is the 2D axisymmetric model, it has been proven lately (next sections) that the flow is not axisymmetric after the transition occurs resulting in totally 3D flow. This shows the need for a 3D simulations instead of 2D plan or 2D axisymmetric ones, despite the difference in computational cost that is related to 3D modeling (more complex meshing and computational time). Further analyses of the flow and the heat transfer are presented in the following sections based on the 3D simulations. PISTON

Experimental validation of the numerical model

The CFD results have been compared with the experimental data of Neu and Subrenat [START_REF] Neu | Experimental investigation of internal air flow during slow piston compression into isothermal compressed air energy storage[END_REF] under the same air compression condition (U pist = 0.033m/s, CR = 4.8, T walls = 293K), so as to validate the used numerical model. Both the evolution of average air temperature and pressure and the local air velocity fields were compared, as reported below. Figure 4.11 shows a comparison on the evolution of air average temperature (T ) between the numerical and experimental results. The two curves show a similar evolution along with the compression process, i.e., all the simulated values of air temperature are within the uncertainty range of the experimental data σ = ±1.6K. The maximum gap observed at the end of the compression (t = 21s) is below 3K. It should be noted that the average water temperature does not show any noticeable evolution over the compression time. Perfect agreement on the evolution of average air pressure can also be observed between the numerical results and experimental data, as shown in Fig. 4.12. At t = 21s of compression time, the piston travels 76% of its length L 0 , the final pressure reaching P f = 4.8 × P 0 .

The air velocity fields obtained by CFD method have also been compared with the experimental data of Neu and Subrenat [START_REF] Neu | Experimental investigation of internal air flow during slow piston compression into isothermal compressed air energy storage[END_REF] using Particle Image Velocimetry (PIV) technique. Fig. 4.13 shows a case of established flow structure at t = 3.15s. Very similar axial velocity profile can be observed between the CFD and experimental results. In contrast, a transi-PISTON

Physical Analysis

Flow and temperature fields description

The evolution of the temperature is a key element to understand the evolution of the compression process. The temperature of the air goes from T 0 = 293K at the beginning of the compression to T f = 325.5K at the end of it at t f = 21s. The process can be qualified as a near-isothermal compression compared to an adiabatic one T f = 525K (See Chapter 1). Different stages are observed in the temperature graph in Figure 4.11. The different stages that the air's average temperature go through can also be seen in the evolution of the flow. These stages are :

• First stage (A), from 0s to 5s, where the temperature increases by about 2K every second, and the curve is concave.

• Second stage (B), from 5s to 8.4s, where the temperature is almost constant and decreases slowly.

• Third stage (C), from 8.4s to 10.4s, where the slope decreases and the average temperature is almost stabilized at 305K at the inflection point.

• Fourth and last stage (D), from 10.2s to the end of compression, where the temperature increases more and more rapidly until the end of the compression where it reaches 325.5K, the curve this time being convex.

These different stages correspond to the flow structures which follow one another during the compression, and are visible on the velocity field. Indeed, in the first stage (A) (between 0s and 5s), the flow presents the first structure, illustrated by various figures at 2s, where the liquid piston is at L pist = 0.066m from its initial position L = 0m. During this first part, the flow is totally axisymmetric.

The air domain's average temperature at this time is 298K, which is 5K higher than the initial temperature of the chamber and the temperature imposed on the walls. A temperature gradient is present at all the walls that surround the air, where the temperature gradually goes from 293K (embedded to the wall) to 300K (about R = 0.010m from the wall). The temperature of the rest of the air domain is uniform at 300K, while water remains uniform at 293K. The first flow structure is as follows: a central current is surrounded by an annular stream moving in the opposite direction (dark blue area in Figures 4.15,4.16). The central stream flow is the fastest: the maximum velocity is 0.096 m/s, or three times the water column velocity.

The first flow structure is as follows: a central current is surrounded by an annular stream moving in the opposite direction (dark blue area in Figures 4. [START_REF] Bouzidi | Momentum transfer of a Boltzmann-lattice fluid with boundaries[END_REF] At the end of this first stage (A), after t = 5s, the water/air interface is L pist = 0.165m from its initial position. The average pressure in the air is P = 1.3 × P f . At this time, the central stream is faster: the maximum positive U z velocity is 0.350 m/s, more than 10 times the velocity of the piston (U z = 10 × U pist ). The higher velocities area of the flow is no PISTON longer towards the piston but has moved towards the cylinder head. The annular stream is also faster: the maximum negative velocity is -0.17 m/s, or 5 times the inlet velocity. At the cylinder head, the central stream is slightly widened. The average air temperature at this time is 302K, which is 9K higher than the initial temperature. The thermal boundary layer present on the walls is thinner, the maximum temperature reached at the center of the chamber now being 310K. The hottest area is closer to the cylinder head than to the interface, it overlaps the faster velocity area. At the air zone near the air/water interface (up to about L z = 0.080m from the interface), the air along the central axis of the chamber is slightly cooler than the air around it (closer to the walls).

The temperature gradients in the flow can be explained by the following phenomenon: the particles that are in the central flow move to the end of the cylinder head, and during this movement they are surrounded by hotter particles since they are the furthest from the walls which are colder (Figure 4.20b). Therefore, they arrive at higher temperature at the top of the chamber and then returns into the annular stream to move in the opposite direction, close to the walls. They cool down during this movement because of their proximity to the walls. These particles return back into the central stream at the water/air interface of the liquid piston which is also colder. These movements create then an radial annular stream flow.

CHAPTER 4. MODELING OF THE COMPRESSION PROCESS USING A LIQUID PISTON

This structure then gradually deforms between t = 5s and t = 9s (stage (B)). In the central flow, ripples starts to be observed throughout the chamber. After t = 9s, the piston has traveled L pist = 0.297m. The average air pressure is 1.5 times greater than the initial pressure. The maximum positive velocity measured during compression is reached : U z = 0.4m/s, and the reverse stream, faster in the area close to the cylinder head, reaches 0.25m/s (Figure 4.18) The average temperature is 305K, which is 12K higher than the initial temperature. The hottest area is always closer to the cylinder head, and the flow in the center of the cylinder, near the interface, is always a little cooler than the air at the piston center (R = 0) (Figure 4.20b). The thermal boundary layer on the walls is still present but is thinner (about 2 × 10 -3 m) in the area close to the cylinder head, and the maximum temperature is 320K. The ripples and deformation of the flow visible on the velocity field are also visible on the temperature field, it is no longer homogeneous (Figures 4.18 The average temperature of the air at the end of compression is T air,ave = 325.5K. The thermal boundary layer is thin, it is hardly visible. The area with the highest temperatures remains the area close to the cylinder head, where local temperatures can go up to T air = 360K (Figure 4.20d). The disruption of the flow can be seen as well through the velocity field as well in the temperature field which becomes non homogeneous.

Velocity profile

The velocity profile evolves during the compression process (Figure 4.21). When the flow structure is established (from 1s to 10s), the flow is separated into two parts, closer to the wall where the velocity is a negative, and closer to the piston center where the velocity is positive. In this part the velocity profile is parabolic. This results in a particular profile, as shown in the figure 4.21. Also, as described in the flow analysis part, the air velocity evolves during the compression process, where it can go from the piston inlet velocity to more than 10 times that velocity (U > 10 × U pist ). During the compression, the air density evolves in concordance with the pressure evolution, resulting in a ρ f ∼ 5ρ 0 . The updated Reynolds number can result into a one that is more than 50 time the initial calculated one in Table 4.4. The results confirm the model where the Re is calculated based on the piston velocity (U pist ) for the air flow modeling is not suitable.

Temperature profile

The heat transfer evolves similarly to the flow. Once the structure is established (after t = 2s), the temperature profile adopts a similar structure, where the temperature at the wall is the one imposed through the boundary conditions and then increases towards the piston center where it is maximal. It should be noted that before the structure disruption, where the air velocities are at their highest values, the temperature profile at the center is slightly concave. The highest radial value is observed at the 0 velocity zone are (R = 0.016m and R = -0.016m) (Figure 4.21,4.15). Figure 4.22 shows the evolution of the temperature profile radially from the beginning of the compression until the structure disruption.

Local velocities at L

z = 0.216m
The temporal evolution of the gas velocity in the z direction (The compression direction) is presented in Figure 4.24 for the axial position L z = 0.216m during compression. During the first moments, the air velocity increases to catch the piston (water) velocity. 1s after, the first negative velocities next to the wall can be seen, which marks the flow structure formation. The higher velocities can be seen in the piston center, and where a certain parabolic form can The figure 4.28 shows the evolution of the U z component of the velocity for several radial positions at the axial position L z = 0.817m. In this position same dynamics, zones and stages are observed as in the position L z = 0.650m, although higher velocities are observed especially closer to the cylinder center. A comparison with experimental results is also achieved. It shows that from t = 0s to t ∼ 6.5s the two curves are coincident. A similar PISTON Figure 4.28: Local velocity through time at L z = 0.817m maximum velocity and profile can be observed. Although, in the experimental results, the velocity starts decreasing once the maximum is achieved (Figure 4.29). Even though, the transition is observed at similar times (between 10 and 11s). This phenomena was observed in all the positions where the comparison was done. This difference can be explained by several factors which are the influence of the experimental setup on the flow, the roughness of the materials used experimentally which can favor the structure disruption. Moreover, a solid piston is used in the experimental study. Numerically, the heat transfer boundary condition at the walls is an isothermal wall, this choice can have an impact lead to a more stable flow and an increased heat exchange retarding the transition inside the flow. PISTON

The origin of structure disruption

The transition that occurs in the flow and the heat transfer which moves it from a structured flow to a chaotic one is a key element in understanding the air compression using a liquid piston. It is shown in figure 4.11 that the evolution of heat transfer is maximal moments before this transition occurs. Neu and Subrenat [START_REF] Neu | Experimental investigation of internal air flow during slow piston compression into isothermal compressed air energy storage[END_REF] made some hypothesis to explain the origin the structure disruption. The numerical study allows a better tracking of the flow. Thus, all the velocity's components (U x , U y , U z ) of the flow were studied to locate the position and the time where the first change in the structure is observed.

The main velocity component is U z , which is the compression direction shows that once the maximum velocity is reached around t = 8s, it starts slowly decreasing before the first non-axisymmetry of the flow is observed at t = 9s. Rapidly after that, the flow's structure is destroyed and the flow becomes a chaotic one (Figure4.28). The flow is totally axisymmetric from t = 1s to t = 6s all over the air domain (except the turnaround zone close the water-air interface and the piston head) (Figure 4.31a). The first perturbations are observed on the U * z velocity component at : t = 8.2s at L z = 0.7m (Figures 4.30, 4.31b). The velocity profile loses its axisymmetric property. The perturbation is then spread in both the directions towards the water/air interface and towards the cylinder head at the same time as shown in the figure 4.31b. It can be seen in Figure 4.30 that higher velocities are observed at higher L z . It is shown in the physical analysis section, both the maximum velocities and the higher local temperatures are observed in the zone between L z = 0.7m and L z = 0.8m. The origin of the flow structure dislocation could be a result of high velocity gradients and the shear that results at the section between positive velocities and negative ones (light blue zones in Figure 4.30). Neu and Subrenat [START_REF] Neu | Experimental investigation of internal air flow during slow piston compression into isothermal compressed air energy storage[END_REF] compared the phenomena to a Kelvin-Helmotz instability.

It remains uncertain whether this chaos is turbulence. Nonetheless, it is a very reproducible phenomena where negative velocities are still dominant close to the wall.

Evolution of air's flow rate

Neu and Subrenat [START_REF] Neu | Experimental investigation of internal air flow during slow piston compression into isothermal compressed air energy storage[END_REF] made an assumption for the evolution of the air flow rate. It is based on considering an isothermal compression where the air pressure increases linearly with the volume of air domain (Eq. 4.20).

Vz (t) = Vpist L 0 -L z L 0 -L pist (t) (4.20)
With Vz (t) is the total flow rate in a section at z height of the chamber, Vpist ∼ 7 × 10 -5 m 3 /s is the flow rate imposed by the movement of the water and L z is the position of the section of interest at the z = L z height.

The assumption is verified while the temperature rise is very slow at the beginning of the compression in all L z positions.

PV diagram

P-V diagram of the air compression by LP shows clearly that it is a near-isothermal one. At the beginning of the compression (stage A 4.11), the actual trajectory is very close to an isothermal one. This can also show that the heat exchanges during before the transition are maximal and enable us to reduce the temperature evolution and thus get closer to isothermal compression (Figure 4.33). We can also see that at the end of the process, this gap starts to PISTON increase rapidly. We then propose to use a variable and adaptive piston speed especially at the end of the compression when the temperature gradients are high. This will enable us to increase the heat exchanges with the external medium. The polytropic index is a key parameter in determining the character of the process [START_REF] Vikram | Energy and Exergy Analysis of Ocean Compressed Air Energy Storage Concepts[END_REF]. It varies from n = 1 for an isothermal compression to n = 1.4 for an adiabatic one.

n = log

P f P 0 log V 0 V f (4.21)
For our studied case, we found that n = 1.09 for the compression. Further parametric study can enable us to enhance the polytropic index as shows the next sub-section.

Parametric Study

In order to better understand the flow and heat transfer evolution and the physics related to compression process using a liquid piston, a parametric study is conducted, firstly by varying the piston advancement speed and then by changing the temperature at the walls which are still considered as isothermal.

Effects of Piston speed

Effects of compression speed have been studied by Mutlu and Kilic through a 2D axisymmetric CFD study [START_REF] Mutlu | Effects of piston speed, compression ratio and cylinder geometry on system performance of a liquid piston[END_REF]. They showed an efficiency decrease while increasing the velocity speed, PISTON this was due to higher power needed for the faster compression. In order to understand the exact effect of piston speed, a comparison between two different U pist is performed, the first is the one presented in the previous sections, while the second one is 3 times faster 4.8. The piston travels the same L pist = 0.69m at the end of the compression. In order to enable comparison, the compression is dimensionless as t * = t t f and both cases have the same final air's volume.

Configuration

slow fast Piston speed U pist (m/s) 0.033 0.1 Compression time t f (s) 21 6.9

Table 4.8: The parameters of the two studied piston speeds Figure 4.35 represents the air's average temperature evolution of these two studied cases. It can be seen that the temperature increases rapidly in the fast configuration, and the final T air,ave gap at the end of the compression is 33K. As seen in the figure 4.35, pressure data were very close, with a small gap at mid-compression time due to volume differences and a similar P f = 5 × P 0 . The results show that higher piston speed does not take advantage of the liquid piston advantages and thus reduce the heat exchanges with the external media resulting in higher air's temperature. Effects of the piston speed on the efficiency are presented in the next sub-section 4.7.1.

Effects of wall temperature

Effects of the external media temperature on the flow and heat transfer have not yet been studied in the literature. Although, as shown in Chapter 1, the convection at the walls is significant and represent one of the biggest advantages of the liquid piston, since it allows a higher L/D ratio and thus a higher heat exchange surface. In order to study the effect of the external media temperature modeled here as isothermal walls, three different T walls = 288K; 293K; 298K are studied (Table 4.9). Same numerical parameters as in table 4. Table 4.9: The parameters of the three studied walls temperature average temperature evolves rapidly to establish a 5K difference between each curve (Figure 4.36). All the three cases show a similar temperature evolution with a transition around t = 10s. Similar stages as in Figure 4.11 can be seen on the two other cases. At the end of the compression (t = 21s), the case at T walls = 288K called "cold" shows a T air,ave = 320.5K which is ∼ 5K lower than the "reference" case and ∼ 10K lower than the "hot" case. This results in a temperature increase of 32K whatever the walls temperature is. This shows that the wall temperature is a very important parameter that influences the compression 137 CHAPTER 4. MODELING OF THE COMPRESSION PROCESS USING A LIQUID PISTON process. It can also mean, that the heat transfer at the wall is dominant. And since it is mainly the convection [START_REF] Odukomaiya | Thermal analysis of near-isothermal compressed gas energy storage system[END_REF], this heat transfer mode is the most important in determining the temperature evolution, in addition to the piston speed and the L/D ratio.

Compression Efficiency calculation

The efficiency measurement of the compression process using a liquid piston has been presented in details in Chapter 2, section 2.2.2. Three models have been tested in this study to compare them and to better estimate the compression efficiency η c . These models are available in the literature and use different compression paths to calculate and compare the energy used for the compression and the energy needed for an isothermal compression. The first model is Yan et al. [START_REF] Yan | Experimental study of heat transfer enhancement in a liquid piston compressor/expander using porous media inserts[END_REF], presented in Chapter 1 :

η c = (P f -P 0 )V iso + P 0 V 0 ln(CR) + 1 CR -1 V f V 0 (P (t) -P 0 )dV + 1 -1 CR (P f V f -P 0 V 0 ) + (P f -P 0 )V iso (4.22)
Where CR = P f /P 0 is the compression ratio, P f , P 0 and V f , V 0 are the pressure and volume at the beginning and the end of the compression. V iso is the volume for an isothermal compression.

The second one, Patil et al. [START_REF] Vikram | Efficiency Improvement Techniques For Liquid Piston based Ocean Compressed Air Energy Storage[END_REF] is also presented in chapter 1, where it is a simpler one, PISTON based only on the compression ratio CR and the polytropic index n [START_REF] Vikram | Efficiency Improvement Techniques For Liquid Piston based Ocean Compressed Air Energy Storage[END_REF] :

η c = ln(CR) + 1 CR -1 CR n-1 n -1 n-1 + CR -1 n -1 + (CR -1) CR -1 n -1 CR (4.23)
Later on, Neu proposed a new model [START_REF] Neu | Etude expérimentale et modélisation de la compression quasi isotherme d'air pour le stockage d'énergie en mer[END_REF]. Unlike the two previous models, its expression includes a term of isochoric and non-isobaric cooling. The efficiency it calculates therefore applies only between the isochoric lines V 0 and V f :

η c = P 0 V 0 ln( V 0 V f ) + V f V 0 -1 - V f V 0 (P (t) -P 0 )dV (4.24)
In order to compare the results, the three available models were used to calculate the compression efficiency of numerical results presented in the previous section and experimental ones realized by Neu [START_REF] Neu | Experimental study of convective heat transfer during liquid piston compressions applied to near isothermal underwater compressed-air energy storage[END_REF], using the same parameters, for the REMORA projects.

Effects of the velocity speed on the compression efficiency

The compression efficiencies calculated by using the three formulas and the parameters from the previous section (4.6.1) according to the piston speed are presented in figure 4. [START_REF] Gong | Numerical investigation of droplet motion and coalescence by an improved lattice Boltzmann model for phase transitions and multiphase flows[END_REF]. We can see that the efficiency increases when the speed of the piston decreases, which is an expected result regarding the thermodynamics of the compression process presented in chapter 2. The gap between the experimental and numerical results especially for the "fast" case can be due to experimental uncertainties that increase when the piston speed increases and mainly at the end of the compression where the temperature and velocity gradient are at their highest level. PISTON

Conclusion

The main objective of this work was to obtain a CFD model of the air flow in the liquid piston. Finite Volume Method resolution of Navier-Stokes equation coupled to VOF method for interface tracking achieved a numerical modeling of the flow and heat transfer inside the liquid piston. The numerical results were validated through comparison with an experimental study. The results obtained with the model enabled a better visualization and understanding of the evolution of the physics of liquid piston's compression process.

• 2D model is not appropriate for modeling the flow inside a liquid piston (different temperature evolution and velocity field). PISTON

• During the first seconds of compression the flow is axisymmetric and a structure where positive velocities are observed closer to the piston center and negative ones are closer to the piston walls. During this phase the temperature increases slowly.

• A transition occurs at t ∼ 10s, it totally destroys the flow structure while this one loses its axisymmetric property and the heat transfer is at its maximum. After that the flow becomes chaotic and the temperature increases rapidly to reach ∆T = 32.5K at the end of the compression confirming the near-isothermal nature of the process.

• The transition is first seen on the compression direction component of the velocity. It happens at the section where the velocity gradients are the highest. Shear stress resulting from frictions between positive and negative streams can be the main reason of this disruption. A 2D axisymmetric can not model appropriately the flow after the transition occurs.

• The Reynolds number based on the liquid column properties is not a relevant criteria to describe the flow evolution. Before the transition occurs air's velocity varies from [-5U pist , 12U pist ]. While Re air ∼ 100 based on U pist , which is relatively low, the flow undergoes several phases, and becomes chaotic after the transition. New criteria could be investigated to characterize the flow's regime.

• Effects of the temperature and piston speed can be seen on the air's temperature evolution and on the efficiency. A faster compression results in higher temperature evolution and thus lower efficiency. Same can be applied to wall temperature. Some differences between the numerical and experimental results are still present, especially the moment the transition occurs. While in this model the wall are modeled as isothermal, in real case, the walls are thick materials. This can have effects on the flow and heat transfer evolution through the compression process. A numerical model that can couple the multiphase flow to the heat transfer in the solid has been studied and is presented in chapter 5.

While the compression is now modeled, the reverse process (i.e Expansion) will be modeled in the next chapter in order to better understand the whole compression/expansion cycle using a liquid piston.

Chapter 5

Modeling of the expansion process using a liquid piston

Introduction

During the restitution of the stored energy, the reverse process as the one we presented in the previous chapter is performed, it is as important as the compression's one. The liquid piston mechanism is assumed to be perfectly reversible so the same organization applies to the expansion of the air allowing the production of electrical energy. It consists of three main processes (Figure 5.1). Water is first injected to fill the piston until a certain level, then the air injected at the storage pressure, the piston valves are then closed and the expansion process starts, where the air pushes the water to generate electricity from the reversible turbine/turbine.

As for the compression process we presented before, only the expansion process is studied here. The thermodynamics of the process has been presented in Chapter 2 with referent studies related to the liquid piston. Description of the approaches and key finding of studies related only to expansion are presented here.

Flow and heat transfer in a liquid piston during expansion process

Rodrigues presented a solid piston expansion process. His work involved a single expansion of air mixed with hydrogen to properly characterize the physics of expansion [START_REF] Jose | Heat Transfer During the Piston-Cylinder Expansion of a Gas[END_REF]. This work also highlights the importance of the Peclet number P e in the study of expansion. The experimentation of CO 2 expansion in a cycle with a compression was developed in the work of Baek [9]. This work was then validated by CFD in the paper that continued the research air compression and expansion machine by a liquid piston ( [START_REF] Török | Quasi-Izothermal Compressors And Expanders With Liquid Piston[END_REF] by comparing the advantages over an equivalent system with a solid piston. The modeling of a 2D expansion of a piston with and without different inserts (honeycomb-like metal foam structure) was studied by Zhang experimentally and through 2D CFD [START_REF] Zhang | Numerical modeling of three-dimensional heat transfer and fluid flowthrough interrupted plates using unit cell scale[END_REF]. The study was limited to short times, but proved the interest of the process. Their study showed that in the case of expansion without metal foam insert, it is advisable to have a slow piston speed at the beginning of the expansion because of the rapid volume evolution (and thus temperature evolution), implying a rapid temperature loss. On the other tested case, an expansion using a liquid piston with metal foam insert, it is recommended to have a fast piston speed at the beginning of the expansion process.

The experimental approach was studied by Yan et al. [START_REF] Yan | Experimental study of heat transfer enhancement in a liquid piston compressor/expander using porous media inserts[END_REF]. It focused on studying the efficiency of different porous media within a liquid piston compression/expansion chamber. Yan et al. study proved that, for a range of insert sizes, the impact of these inserts was very significant on the thermodynamic efficiency during expansion, bringing the the phenomenon closer to an isothermal case. Later on, Wieberdink improved the previous study by extending it to high pressure tests (210 × 10 5 P a).

Objectives and main contributions

The aim of this numerical simulation of the expansion process is to complete the compression process studied in the previous chapter in order to enable a complete study of the flow in the liquid piston (a complete cycle). As in the compression study, experimental results realized under the same project and from the literature allow us to validate the numerical model PISTON through a comparison analysis. Few studies exist on the expansion process, especially those dedicated to flow analysis. The present study will give a better illustration of the flow and heat transfer evolution inside a liquid piston during the expansion process.

Geometry and numerical parameters

The same model as the one we presented in Chapter 4, section 4.2 is used for modeling the expansion process using a liquid piston. Only a 3D model is studied for the expansion, since we showed that the flow in the compression process is a totally 3D one, it is thus expected a similar flow's nature for the expansion's one.

A similar geometry of the one studied for the compression process is used. For the 3D geometry, the liquid piston is modeled as a cylinder with L 0 length and D ∈ [-R, R] diameter. The water is ejected from the piston bottom at a velocity corresponding to the experimental one U pist = -0.033 m/s (Figure 5.2). L pist is the position of the water/air interface at time t. The piston walls are considered as isothermal at T walls and the domain is at an initial temperature T 0 . The initial domain filled by air from L pist,0 to the top of the piston while water is in the bottom and is ejected from at U pist . Numerical and physical parameters of the simulations are listed in table 5.1 : piston geometry, piston velocity (U pist ), compression time (t f ), Expansion Ratio (ER = P 0 P f ) and piston walls temperature (T walls ).

D (m) L 0 (m) L pist,0 (m) U piston (m/s) t f (s)
P 0 ER T walls (K) 0.0518 0.906 0.693 -0.033 20.6 4.4 × 10 5 4.4 300

Table 5.1: Physical and numerical parameters of the liquid piston expansion process

Mesh convergence study for the expansion process model

Grid independence study is conducted first to determine the correct mesh size and the mesh distribution in the numerical simulations. After proving its efficiency in the compression process modeling, only hexahedral type of meshes is used. Different parameters are used to determine the mesh convergence : average air temperature, local velocity at defined points. A structured O-Grid mesh is used, with mesh refinement at the wall and in the top half of the compression chamber where air is present at the beginning of the expansion and where the highest gradients were observed (Figure 4.6a). We showed that the most important evolution of the flow happens during the first 10s of the expansion. In order to reduce the computational time, the mesh convergence tests were computed for the first 10s only (t f = 10s). The walls temperature is T walls = 293K and the initial pressure is the one issued from the compression process studied before P 0 = 4. Table 5.2: Types of 3D meshes used in the mesh convergence study

The four studied meshes show close results until about t = 5s. From this time, the simulation with the least number of meshes diverges significantly from the other results. This is mainly due to the fact that this mesh is not enough refined at the walls. But as shown for the compression study, the mesh close to the wall is a key parameter in having a good meshing and accurate results. It is therefore estimated that the 830k mesh is the one that optimizes the computation time without constraining the accuracy of the results and is used for the rest of the study. The obtained results enabled us to validate the numerical model, which will be used for further analysis of the flow and the heat heat transfer of the expansion process.

Analysis of the flow and heat transfer in the liquid

piston during the expansion process

Description of the velocity and temperature fields

As in the compression process, no significant flow's property is observed in the water column, the main domain of interest is the air's one. The evolution of air's temperature during expansion can be described by four main evolutionary profiles (Figure 5. The evolution of the air's temperature during the expansion can be explained by comparing the influences of interface/near wall convection and surface exchanges. As the expansion starts, the air's density increases, due to the fact that the volume increases (and pressure decreases) and the heat exchange surface L/D increases. Therefore, there is a transition where the convection phenomena is the most dominant one, resulting in a higher heat exchange at the walls. The final temperature of the air is only 7K lower than the initial one. The temperature of the walls being considered constant at 300K and is greater than that inside the liquid piston, convection results in increasing the air's temperature.

The zones presented in the graph of air temperature's evolution can also be seen on the pressure's on (Figure 5.4). In fact, during the first second, the pressure decreases rapidly (P = 1/2 P 0 at t = 4.5s), it then takes another 10s for it to decrease by half again.

Investigation of velocity's local evolution

In addition to global evolution of the expansion parameters, numerical study gives us a view on the local variables (velocity, temperature, pressure...). The Longitudinal position In the second stage (B) which is located between 1.2s and 3.2s, the flow remains axisymmetric. The maximum positive velocity reached U z = 3 × U pist while the maximum negative velocity U z = -7 × U pist . The recirculation zone evolves in size and moves between 0.019m and 0.012m. The U z velocity component at R = ±0.019m shows strong gradients then the flow seems to recover like the recirculation zone. A second zone of low or zero velocities is also established for a short time near the walls. This is in agreement with the results obtained on the figure previously obtained in the section.

From t = 3.2s (Stage (C)), the flow is no longer axisymmetric and is agitated, which affects the position of the recirculation zone, as shown by the study of the points R = ±0.019m. The velocity of the other points seems to settle around half of the maximum velocities previously found. PISTON P-V diagram shows the paths that the process follows 5.14. We can see that the at the beginning of the expansion process, the actual trajectory is closer to the adiabatic one. This corresponds to the first zone (A) 5.5, where the temperature drops rapidly. As the process advances, the studied case becomes closer and closer to the isothermal. At the end of the expansion, the gap between the two is almost null. While a constant velocity is imposed at the outlet, P-V diagram shows that in order to get closer to isothermal expansion, we should consider implementing a variable piston speed, with very slow expansion at the beginning which can increase at the end of the expansion. It has been also proved by Yan et al. [START_REF] Yan | Experimental study of heat transfer enhancement in a liquid piston compressor/expander using porous media inserts[END_REF]. As for the compression process, we calculated the polytropic index n (Eq. 4.21). We found that n = 1.012. Which proves again that that the expansion using a LP is near-isothermal one.

Effects of wall temperature

Effects of the walls temperature (external medium) have been proven in the compression study. A similar study is conducted to measure the effects of the walls temperature on the evolution of the expansion process. The same geometry, boundary and initial conditions as in Table 5.1 are implemented. Only the walls temperatures have been changed (T walls = 293K and T walls = 300K) (Table 5.3).

Configuration

Cold Hot Walls temperature T walls (K) 293 300 Expansion ratio ER 2.5 2.5 12.5 × 10 5 P a Expansion ratio ER 5

Table 5.4: The parameters of the two studied velocity effects based on Yan et al configuration and keeps increasing until the end of the process where it reaches ∆ = 14K. Slower the expansion is, lower is the temperature evolution and closer it will be to a near-isothermal process. The results are concordant with the one obtained in the compression study. A slow compression/expansion enable to increase the heat exchanges between the air and the walls (the external medium) and thus reduce the temperature evolution caused by the pressure drop.

Figure 5.16: Comparison between air's temperature evolution for two different U pist

Expansion Efficiency

The efficiency measurement of the expansion process using a liquid piston has been presented in details in Chapter 2, section 2.2.2. Two models have been tested in this study to compare them and to better estimate the expansion efficiency η e . These models are available in the literature and use different expansion paths to calculate and is the ratio of the energy needed for an isothermal expansion and the one used for the expansion. The first model is Yan et al. [START_REF] Yan | Experimental study of heat transfer enhancement in a liquid piston compressor/expander using porous media inserts[END_REF], presented in Chapter 2 :

η e = W e actual E p = Ve V 0 (P (t) -P f )dV + (P 0 -P f )V 0 P 0 V 0 (ln(ER)) -(P f (V iso -V 0 )) + (P 0 -P f )V 0 (5.
2)

The second one, Patil et al. [START_REF] Vikram | Efficiency Improvement Techniques For Liquid Piston based Ocean Compressed Air Energy Storage[END_REF] is also presented in chapter 2, where it is a simpler one, based only on the compression ratio ER = P 0 /P f and the polytropic index n :

η e = 1 - ( 1 ER ) n-1 n n-1 -1 ER n-1 n + 1 ER ln(ER) + 1 ER -1 (5.3)
The above models are used to calculate the expansion efficiency of the main case and the one presented in the parametric study. As for the compression efficiency which follows the increasing temperature profile over time, the expansion efficiency also follows the temperature evolution. The closer it is to an isothermal one, the higher η e will be.

The final efficiencies, plotted in Figure 5.17 and Figure 5.18 show an increase in the expansion efficiency as the expansion time increases. This confirms the explanations given earlier in the previous section : a slow process increases the convective heat exchange between the piston wall and the air, thus compensating the temperature drop caused by the pressure drop leading to an increase in efficiency. On the other hand, for Patil's modeling (Figure 5.17) the efficiency value for the 8s Yan's 8s case is higher than that of the 10s expansion with the main configuration. This discrepancy is initially caused by a ratio of 3 between the lengths of the two configurations L 0 . It is also linked to Patil's model, which is based only on the expansion ratio at the end of the process. PISTON

Compression-Cooling-Expansion cycle

In a operating CAES system the storage/restitution cycle is composed of several steps. As shows Figure 5.19a, the first ones happen during the storage phase, where the air is compressed and and ejected to storage tanks. During the energy production (restitution) phase, the piston valve is opened to inject the compressed air into the piston (5.19b). Then, after closing the valves, the expansion process starts and the water pushed by the air goes through a turbine to produce electricity. For our studied case, we modeled a simpler cycle presented in Figure 5.20, the cycle is operated in a closer piston, where the air is first compressed. During the compression, the air's temperature increases. At the end of the this process, the air stays in the top of the piston to achieve an isochoric cooling. When the air returns to its initial temperature, the expansion process starts. At the end of the cycle, the air's final volume V f is equal to the initial one V 0 . The study of compression-cooling-expansion cycle numerically and experimentally shows a great similarity in the global evolution. But some disparities remain present.

Cycle's PV diagram

The PV diagram shows a path that goes from near-isothermal compression, isochoric cooling with small pressure drop and near-isothermal expansion that also includes a heating process (Figure 5.24). The cycle is then closed when P 0 ∼ P f and V 0 = V f . PISTON

Cycle Efficiency

The efficiency of the cycle η cycle can be calculated as [START_REF] Vikram | Modeling of liquid-piston based design for isothermal ocean compressed air energy storage system[END_REF][START_REF] Vikram | Energy and Exergy Analysis of Ocean Compressed Air Energy Storage Concepts[END_REF] :

η cycle = η c+cool × η e (5.4)
Yan and Neu's models allow us to calculate the isochoric cooling process efficiency combined to the compression one η c+cool [START_REF] Yan | Experimental study of heat transfer enhancement in a liquid piston compressor/expander using porous media inserts[END_REF][START_REF] Neu | Etude expérimentale et modélisation de la compression quasi isotherme d'air pour le stockage d'énergie en mer[END_REF]. Using Yan models (4.22, and 2.7), we found that :

η cycle = 0.95 * 0.931 = 88.5% (5.5)

Cooling process flow evolution

The investigation of the isochoric cooling stage has not yet been documented. Therefore it is interesting to describe and analyze the heat transfer and flow during this process. This would also better understand the impacts that can have this process on the expansion following it. At the end of compression, the air has a gradient of 70K between the interface and the top of the piston (figure 5.25).The air flow shows large agitations at the beginning of the cooling, which are resulted from the compression flow. This distribution is maintained until CHAPTER 5. MODELING OF THE EXPANSION PROCESS USING A LIQUID PISTON 4.5s after the beginning of the cooling (Figure 5.26). After that, the flow is at "rest", where the velocity field is almost null.

At the end of the cooling, a small temperature gradient is still present (Figure 5.27). This is a difference from the simulation of the expansion alone, but has a very small impact on it, since only 1K for initial temperature is not significant in comparison to T walls .

Conclusion

The main objective of this chapter is to provide a 3D study of the flow of air expansion by liquid piston. This objective was achieved by using the same numerical model as the one 165 CHAPTER 5. MODELING OF THE EXPANSION PROCESS USING A LIQUID PISTON used for the compression study. In parallel to this work, internal experimental tests were performed to compare the different results.

• As in the compression, during the first seconds of expansion, the flow is axisymmetric and a structure where negative velocities (in the direction of the piston displacement) are observed closer to the piston center and positive ones are closer to the piston walls. This stage lasts t ∼ 3.2s, which is three times shorter than in the compression. During this phase the temperature drops rapidly.

• A transition occurs at t ∼ 3.2s, it totally destroys the established flow structure. and the temperature reaches its lowest value T = 280K. After a short temperature free (1s); the temperature starts increasing due to the convection with the walls, and reaches T f = 293K, 7K lower than the initial temperature. The process can be qualified of a near-isothermal one.

• Effects of the temperature and piston speed can be seen on the air's temperature evolution and on the efficiency. A faster expansion results in higher temperature drop and thus lower efficiency. Same can be applied to wall temperature. The different stages that the flow undergoes are independent from the wall temperature or the piston, while the temperature evolution changes, the same pattern can be seen.

A thermodynamic cycle consisting of compression-isochoric cooling-expansion is modeled and the flow during the cooling is studied.

Chapter 6

Conjugate Heat transfer in the liquid piston 6.1 Introduction

Physical problem and technical barriers

During the compression process, the temperature of the air increases. Due to the slow compression time, the walls store a quantity of heat generated by the flow. In the previous chapter, the piston walls were modeled as a constant temperature boundary condition. While in real case, the interaction between the solid walls and the flow inside the piston can have an impact on the evolution of the process. In order to take into account this physical phenomenon and to study the influence of solid walls on the flow, the coupling between the solid domain and the fluid one is necessary (Figure 6.3). Thus two approaches have been developed, adapted and tested on OpenFoam in order to realize this modeling with thick walls. This constitutes a technical challenge as no similar model that couples multiphase flows solvers to heat transfer in solid one's.

Literature review of CHT modeling

Conjugate Heat Transfer (CHT) analysis is a simulation process that addresses the thermal interaction between a body and a fluid flowing over or through it. CHT problems occur whenever fluid convection and solid material conduction are taken into account simultaneously. This problem is classified as a multi-physics phenomenon which is characterized by at least two distinct computational sub-domains (fluid and solid regions), whose associated equations are coupled at shared interface.Transient thermal interactions between a fluid and a solid are difficult to perform because of their computational cost. Indeed, fluid and solid characteristic times differ by several orders of magnitude [START_REF] Moretti | Stability, convergence and optimization of interface treatments in weak and strong thermal fluid-structure interaction[END_REF].

The behavior of the interface conditions in CHT is also often studied using a normal mode analysis [START_REF] Dorfman | Conjugate problems in convective heat transfer[END_REF]. Roe et al. proposed an explicit coupling for CHT [START_REF] Roe | Combined interface boundary condition method for coupled thermal simulations[END_REF]. An interesting alternative is the steady-state approach described by Verstraete and Scholl [START_REF] Verstraete | Stability analysis of partitioned methods for predicting conjugate heat transfer[END_REF], it showed the utility of a simpler model to predict the CHT.

Solving the heat transfer in solid domain

In this section the equations used for solving heat transport in the solid are presented. The solid domain is modeled using unsteady heat transfer equation:

ρ s c p ∂T ∂t -∇.(λ∇T ) = Q (6.1)
where ρ s is the solid density, c p the specific heat capacity, λ the thermal conductivity and Q is a volumetric heat source. Note that for turbulent cases, λ is substituted by λ ef f , the effective thermal conductivity which take into account ν t the turbulent thermal viscosity, and where κ = λ ρcp is the thermal diffusivity [START_REF] Cheung | Conjugate Heat Transfer with the Multiphysics Coupling Library preCICE[END_REF].

λ ef f = κ ef f ρc p = (κ + κ t )ρc p = ( ν P r + ν t P r t )ρc p (6.2)
To simplify the interactions between the domains, instead of solving directly the temperature field, it is preferable to solve Equation 6.1 for total enthalpy, ∆H = c p ∆T , as shown in equation :

ρ s c p ∂H ∂t -∇.(λ∇T ) = Q (6.
3)

The thermal boundary conditions that can be applied to both the fluid and the solid domains are explained in the appendix. It can be a temperature (Dirichlet), a flux (Neumann) or both (Robin) 8.0.0.7.

Conjugate Boundary Condition

In conjugate heat transfer problems, boundary conditions of the fourth kind, or interface boundary conditions, state the continuity of the temperature and the heat flux at the fluidsolid interface Ω F S = Ω F ∩ Ω S . This reads :

T S = T f on Ω F S (6.4) k S ∂T s ∂n + Q = k F ∂T F ∂n on Ω F S (6.5)
In the Robin-Robin coupling, we assume that Newton's law of cooling applies:

q = h(T -T ext ) (6.6)

New implemented coupling approaches for solidfluid coupling

Numerical algorithms for CHT problems and multi-physics in general, are usually divided into two classes: partitioned and monolithic. While the first allows for direct use of a specialized solver in each region and consequently offering advantage of code reuse, the latter outperforms in terms of efficiency and stability. In partitioned strategy, fluid and solid solvers are called sequentially and accurate enforcement of coupling conditions necessitates a predictor-corrector loop which entails additional computational cost [START_REF] Sosnowski | Numerical model for thin liquid film with evaporation and condensation on solid surfaces in systems with conjugated heat transfer[END_REF]. In contrast, monolithic approach adopts simultaneous solution of identical physics, i.e. energy equation, throughout entire domain to accelerate convergence rate.

Model 1

To model the thin thermal walls, a temperature-dependent thermal-baffle boundary condition (BC) was developed. This was based on an existing type of BC called ThermalBaffle1D and was extended to 3D model with a coupling to the compressible multiphase fluid solver. Note that for a steady state heat transfer, the model can be applied directly through an interface boundary condition. The considered model is based on Cheung work [START_REF] Cintolesi | Large-eddy simulations of conjugate heat transfer with evaporationcondensation and thermal radiation[END_REF] and was extended to 3D unsteady heat transfer in the solid. The model has the advantage of being implemented directly in the fluid solver. The solid domain can be simply generated through an extrusion of the fluid domain. Which guaranties having conformal meshes. The heat transfer model in the solid is then solved through a multi-region' case. The non-linear diffusion coefficient κ ef f , which is additionally a function of the dependent variable T , is solved explicitly. Its value was calculated using the solution from the previous time step.

The main advantage of this method is that it generates instantly a conformal mesh (Figure 6.5). It also adopts the same time step for the fluid and solid solvers. This is a critical parameter, since, for the fluid, an adaptive time step is chosen due to the evolution of the velocity gradients during the compression/expansion processes.

Model 2

The second approach couples two already existent solvers in OpenFOAM. Two domains, with separate cases (mesh, solver...) are coupled through boundary conditions at each time step. For the fluid the model presented in previous chapters is used, while for the solid Eq. 6.1 is solved for unsteady heat transfer. The coupling works similarly as in the already existent solver "chtMultiRegionFoam", while the main contributions are the applications for a multiphase compressible solver for the fluid domain.

The idea of this approach is that each solution for one of the domains produces a boundary condition along the interface for the other. The process starts by solving the fluid domains in assigned order, with an initial guess of the temperature distribution at the interface (explicitly). The heat flux distribution obtained at the interfaces is then used to solve the energy transport in the solid domains and to obtain a new temperature distribution, and so on 6.4. If this process converges, the iterations continue until a desired accuracy is achieved. Further description of the coupling model can be found in [START_REF] Mousavi | Combination of reactingFoam and chtMultiRegionFoam as a first step toward creating a multiRegionReactingFoam, suitable for solid/gas phase[END_REF]. This model is more general than the first one, it could be used for complex solid geometries and can also model a radiation heat transfer model, which is not the case for the model 1. For 3D modeling a MPI parallel computation is necessary. Thus, the partition between the different should be done manually for this type of coupling. This method of decomposition decreases the efficiency of the parallel computation.

Results of the compression process

A set of simulations are conducted to test and compare the new models and to compare them to the isothermal walls model, and the experimental setup. In order to enable this comparison, the same material (Borosilicate glass) and wall thickness t w = 3.5mm as the one used experimentally is adopted.

Comparison of the two models

Comparison between the average temperature of the different models show no significant difference all over the compression process ∆T < 0.5K (Figure 6.6). As mentioned in the introduction, the temperature in the solid domain evolves very slowly at the beginning of the compression, this is due to the thermal inertia of the solid which have very large time scale compared to the fluid one [START_REF] Dorfman | Conjugate Problems in Convective Heat Transfer: Review[END_REF]. At the end of the process t = t f , the maximum solid temperature T s ∼ 297K locally, while T s,ave = 293.35K. The water domain cools the solid as it passes by, resulting in T s = 293K across the water domain, while in the solid zone in front of the air domain, the solid average temperature is T s,ave = 294.1K (Figure 6.7). The higher temperature are located at the piston top where the flow is almost null. It proves that convection between the solid and the fluid is the most dominant. For the flow evolution, the comparison of U z velocity field between the initial model (isothermal) and the two new proposed models shows similar results (Figure 6.8). An important difference, is that the structure disruption is slowly delayed by 0.4s as shown in Figure 6.9a. It can be seen that in the isothermal model, the transition stage has already began, while the structure is still established in the two CHT models. While this can be due to the effect of the walls temperature, it is more likely a numerical issue. Indeed, it was shown in the compression study that the effect of mesh refinement near the wall is a critical parameter, while for coupling reason, the mesh domain at the coupling interface is mapped to ensure connectivity with the solid domain. This can results in a local numerical perturbation that does not ensure an accurate no-slip boundary conditions especially for high gradients. Globally, the velocity, the temperature fields and their evolution, and through comparison, can enable us to validate the two implemented models.

Parametric study of the compression process

External Temperature

After modeling the CHT transfer between the solid and fluid domain, the interaction between the solid domain and the external medium has not yet been documented. To model it, three contrasting boundary conditions are tested (cf. Appendix 8.0.0.7). The first is the T ext = 293K called isothermal boundary conditions (Dirichlet) [START_REF] Jacobs | A Conservative Isothermal Wall Boundary Condition for the Compressible Navier-Stokes Equations[END_REF], which had been used previously where the temperature is imposed, the second one is the adiabatic boundary condition (Neumann), where the heat flux is null. Finally the third one is a BC that comes between the two, where a heat flux is imposed with h = 10 W m -2 K -1 (Neumann) (Equation 6.6). .9 shows the evolution of the temperature in the solid domain during the compression process. It can be seen that until t = 9s a very small gap of ∆T s = 0.02K can be seen. At the end of the compression this gap increases to ∆T s = 0.56K between the isothermal and the isoflux. The gap between the adiabatic and isoflux cases ∆T s = 0.17K is opposite to the expected results, but can be explained by the fact that an adiabatic case is more constraining one which can lead to numerical instabilities. In addition, the gap of local maximal temperature is very small and can be due to many other numerical, coupling reasons.

Materials

In addition to Borosilicate glass which was used in the experimental study and was presented for the model's validation, a parametric study including different materials with contrasting physical properties have been studied. The chosen materials have very different extreme properties in terms of conduction λ, density ρ s , heat capacity c p and thermal effusivity e = λρc p . These characteristic are the main ones that can influence the evolution of the temperature in the solid (Table 6

(Jm -2 K -1 s -1/2 )
1.4 × 10 3 3.5 × 10 4 3.5 × 10 2 Table 6.1: Physical parameters of the three studied materials Figure 6.10 shows the evolution of air's temperature during compression process with three different solid wall materials. While copper and Borosilicate glass show no significant difference, the less thermally diffusive material (Plexiglas) shows a significant impact. This is due to its low density which makes its thermal inertia (effusivity) small, thus it can not dissipate as much heat as the copper or the Borosilicate glass.

Expansion study with CHT

After completing the compression CHT investigation the expansion process is as much interesting, since the air's temperature evolution is higher and since the heat transfer is maximal at the end of the process as shown in the isothermal study. Same material and simulation properties as in the experimental and isothermal cases are imposed. Figure 6.11 shows a comparison of air's temperature evolution during expansion process, between the initial model (isothermal walls), the CHT models 1 and 2. The CHT model shows a small gap with the isothermal one at the transition moment t = 4s where a gap of ∆T air,ave = 1K can be seen. The different zones shown in previous chapter are also present and the final temperature is almost identical.

The solid wall effects can be most seen at the moment where the temperature is at its lowest and where its gradients are at their highest level. It is expected as has shown Dorfman [START_REF] Dorfman | Conjugate problems in convective heat transfer[END_REF] that at higher temperature differences, effects of a material solid wall can be more significant. While the CHT model results did not get closer the experimental ones, both the two numerical models remain inside the margin error.

Conclusion

In order to take into account the effect of the evolution of the wall temperature on the flow, a numerical study of CHT (Conjugate Heat Transfer) between the fluid domain and the piston walls is carried out. For this, two coupling approaches have been developed and implemented into the numerical model. Comparison of the two models show very few differences between the two CHT models. The model 1, which is more computationally efficient and easier to implement, was chosen for the rest of the wall effect study. The evolution of the temperature in the solid is well present, but is limited and occurs at the end of the compression. It has a small effect on the flow and heat transfer of the air's domain. A parametric study of different walls materials and external temperature boundary conditions is carried out in order to better quantify the effect of the walls on the liquid piston compression process. Three contrasting external boundary conditions (adiabatic, isothermal, and a specified heat flux) have been tested. Their effects can be seen on the wall's average temperature especially at the end of the compression, but it has small effects on the fluid domain. It can be seen that the most significant parameter is the materials conductivity. Comparative study of three various materials shows that Plexiglas walls, which is the least thermally diffusive studied one had a significant impact on the air's temperature (∆T = 3K). The expansion process was also studied using the CHT model 1 and results were compared to the initial model. We showed that, during the transition, where the temperature is at its lowest values, a gap of ∆T = 1K is seen, this can be due to dissipation of energy in the solid walls (thermal inertia). and tested. Firstly Shan-Chen modified approach for Single Component MultiPhase and Multi-Component MultiPhase flows have been used to model high density ratios are tested. The approach suffers from spurious currents especially for the MCMP flows resulting in unstable schemes. A second approach "Free-Energy" is implemented in a in-house code using two models "Lee-Lin" and "Inamuro", the first enabled to simulate a static droplet at 1000 density ratio but could not enable a model a dynamic case, while the second model, which has a very high computational cost enabled a simulation of two droplet collision with high Re and average density ratio [START_REF] Higuera | Boltzmann approach to lattice gas simulations[END_REF]. Finally, a novel VOF method for two-phase flow is tested. It allowed simulation of high density ratio/high viscosity ratios multiphase flows. The method suffers from lack of mass conservation. A decision was made to switch to a more classical approach which is presented in the later chapters. Chapter 4 concerns the modeling of LP air compression using Finite Volume Method for solving Navier-Stokes equations and Volume Of Fluid method for interface tracking. The numerical model was compared and validated through comparison with experimental results carried out inside REMORA project. The 3D nature of the flow is demonstrated. The flow undergoes several stages, the most important one is the establishment of typical axisymmetric flow structure. We found that during this period air's velocity varies from [-5U pist , 12U pist ], which makes the Reynolds number based on U pist not a reliable non-dimensional number for characterizing the flow evolution. Then a transition occurs which disrupts the structure and makes the flow totally agitated. The transition is first seen on the compression direction component of the velocity. It happens at the section where the velocity gradients are the highest. Shear stress resulting from frictions between positive and negative streams can be the main reason of this disruption.

A parametric study was conducted to measure the effects of compression parameters on the flow evolution and the efficiency. A faster compression results in higher temperature evolution and thus lower efficiency (drops from η c = 91% for U pist = 0.033 m/s to η c = 82% for U pist = 0.1 m/s ). The walls temperature (ambient external medium), is a key parameter for the air's temperature evolution. A higher walls temperature results in a higher air's temperature and thus a lower compression efficiency.

Chapter 5 studies the expansion process and same analysis as in Chapter 4 is conducted. Results show similarity of the flow evolution with the compression one. The main difference is the temperature evolution, which decreases at the beginning and then increases after the transition occurs. A parametric study is conducted and shows the important effect of expansion time in favoring the convection and thus temperature reduction and enhancing efficiency (up to η = 94%). A compression-isochoric cooling-expansion cycle is also studied. We showed that the near-isothermal expansion results also in a heating process, and thus the cycle is closer to an isothermal one since the polytropic index n = 1.012 for the expansion. The cycle's efficiency is also very high compared to an adiabatic one, as it achieves η cycle = 88.5%.

In chapter 6, two approaches are proposed to model the Conjugate Heat Transfer between solid walls and fluid (air and water). Both explicit approaches shows similar results. While few differences are observed with the previous model, effects of solid material (Plexiglas, copper and Borosilicate glass) were highlighted. We found that the material with the smallest effusivity has a significant impact on air's temperature evolution (∆T = 3K).

The liquid piston offers high compression/expansion efficiencies (91% to 95%), with its suitability for under-water applications, it constitutes a reliable technology for Compressed Air Energy Storage applications and mainly offshore ones. It also offers the possibility to use a parallel compression/expansion stage with several LP working in sequence. We showed in this work that we can achieve near-isothermal processes by optimizing their parameters (Piston speed profile, external medium temperature, solid wall materials...).

Perspectives

A better visualization of the flow inside a LP is now possible through the numerical model. While, in this thesis, a fine description is given for the flow structure establishment, the origins of it are yet to be determined. Its disruption is still a key element for further studies. A more in depth physical analysis is still needed to better understand the flow evolution. The effects of piston parameter (piston speed, wall temperature) can be seen on the temperature evolution. Heat exchanges that occur at the walls need to be measured through Nusselt number or heat transfer coefficient calculation. Further studies on this topic would help understand the effect of each parameter. As it has been shown in this thesis, the Reynolds number is not representative of the flow evolution. New dimensionless criteria for heat transfer estimation and transition prediction would help reduce and simplify the numerical model, which is a computationally heavy one.

The CHT developed models offer the possibility to study a more complex geometry (foam, porous media) and their effects on the flow. Other phenomena that have not been modeled such as mass transfer and phase change could be studied. Indeed, the local temperature observed at the end of the compression and especially in the expansion before the transition may result in a phase change of the air. Interaction between the gas domain and the fluid one, at high pressures could also provide guidelines to a better understanding of the flow.

A more on depth study of the effects of piston geometry (shape, L/D), solid wall materials and compression/expansion parameters are necessary to determine precise effects of each one. This can be realized through a topology optimization which would help increase the efficiencies of the compression/expansion processes using a LP for CAES applications.

Finally, simulations of LP at industrial scales and high pressure ratios, with multi-stages and successive cycles will ensure the optimized configuration of the LP for a CAES applications.

Molecular dynamics (Microscopic approach)

Realistic molecular level fluid flow simulation can be performed using the Molecular dynamics approach. The molecules are modelled individually allowing them to move and collide with respect to Newton's equation of motion. Simulation of this method is carried out using Potential functions to represent the molecular interaction and are solved deterministically. To represent a system with the big domain is difficult, as the simulation requires a large number of molecules which in-turn requires massive computation cost. Hence, molecular dynamics are used to represent small simulation which includes a feasible number of molecules to simulate. To simulate these dynamics, distribution functions are required which will be discussed in the next section.

Distribution Function

The distribution function is a statistical approach to represent the ensemble of many copies which gives the probability of finding the molecules with given position and momentum. It is enough to describe all the gas properties until it does not rely on the relative position of the molecules. If there are no collisions, the new position [START_REF] Michael | Lattice Boltzmann Modeling: An Introduction for Geoscientists and Engineers[END_REF] can be determined by:

f (x + dx, e + de, t + dt)dxde = f (x, e, t)dxde (8.1)
If there is a collision of the molecules, the molecules at the new position [START_REF] Michael | Lattice Boltzmann Modeling: An Introduction for Geoscientists and Engineers[END_REF] can be determined by:

f (x + dx, e + de, t + dt)dxde= f (x, e, t)dxde + Γ (+) -Γ (-) dxdedt (8.2)
Where, f represents the particle distribution function (PDF) which depends on the space(x), time(t) and the velocity of the particles (e). Γ (+) represents the particles added to the function due to the collision and Γ (-) represents the particle removed from the function due to collision.

LBM (Mesoscopic approach)

Lattice Boltzmann equations are derived from a method called Lattice gas automata which are simulated with the Boolean operations using the Boltzmann equations. LBM can be considered as one of the alternate methods to restore the Navier Stokes (N-S) equations, which can be demonstrated using Chapman-Enskog expansion [START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF].

There are different methods of obtaining the Lattice Boltzmann equation form the Boltzmann equation. One of them is by discretizing the Boltzmann equation using Hermite Polynomial series. The approximation of spatial and temporal discretizations are based on truncated Taylor expansions. Whereas, the molecular velocity is an independent variable in PDF.

LBM is mainly based on the discretized quantity of PDF which depends on the space, time and velocity parameters. The order of truncation of Taylor series was selected to recover the Navier Stokes equation based on Chapman-Enskog analysis.

The LBM methods are solved by discretization of the domain by lattice and the number points or directions in the lattice of the domain know as lattice structures. The type of lattice structure discretization is determined by the scheme that is used to solve the LBM. The grid spacing between each lattice point is assumed to be equal in all directions δx = δy = δz, which allows for easy implementation and numerical stability. The lattice provides the coupling between the time, space and velocity discretization to solve LBM. For implementing LBM, selection of structure of the grid, which is the lattice structure and basic principle of LBM, which is streaming and collision steps will be discussed in the coming sections. The lattice structure will be decided based on the problem that needs to be simulated and the accuracy of the simulation. Lattice structure is commonly represented as "DdQq" where 'd' represents the dimension of the simulation (2D or 3D) and 'q' represents the lattice direction allowed for the particle to travel, Higher the lattice direction, higher the accuracy of the simulation but also increases the computation cost. The commonly used lattice structures are D2Q5, D2Q9 lattice for 2-dimensional simulation and D3Q15, D3Q19, D3Q27 lattice which are used for 3-dimensional numerical simulation.

The lattice point has weights w i which are accounted for computation, as not all lattice vectors have the same length. Most of the lattice structures have 3 different weights where w 0 corresponds to 0 velocities, w s corresponds to short velocities and w l represents the long velocities. LBM is inspected using the truncated Taylor expansion and a multi-scale analysis to obtain the desired partial differential equations. The lattice must have a set of symmetry conditions to preserve the physics of flow [START_REF] Latt | Straight velocity boundaries in the lattice Boltzmann method[END_REF]. Where w i is the lattice weights, e is the velocity of the molecule, c s is the speed of sound in the lattice and δ αβ is the Kronecker delta.

Streaming and Collision

Streaming and collision are the two basic steps to simulate LBM, where streaming step includes the moment of molecules from one position to another and the collision step describes the position of molecules during the collision. Here, the basic model of Single Relaxation Time (SRT), Bhatnagar Gross and Krook (BGK) [START_REF] Bhatnagar | A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems[END_REF] will be discussed.

f a (x + e a ∆t, t + ∆t) = f a (x, t) - [f a (x, t) -f eq a (x, t)] τ (8.9)
The first part of the equation 8.9 f a (x + e a ∆t, t + ∆t) = f a (x, t) is the streaming step and the second half of the equation [fa(x,t)-f eq a (x,t)] τ is the collision step. Where f a is the PDF which depends on space and time, τ is the relaxation time which is determined based on the viscosity of the fluid.

The commonly used collision model is known as BGK collision model where the equilibrium distribution function f eq a can be obtained from using the macroscopic velocity u, macroscopic density ρ(x), molecular velocity e a and lattice weight w a as shown in the equation below.

f eq a (x, t) = w a ρ(x) 1 + 3 e a • u c 2 s + 9 2 (e a • u) 2 c 4 s - 3 2
u 2 c 2 s (8.10)
The lattice weights w a are introduced to maintain the symmetry of the lattice. c s is the speed of sound in lattice which is defined as c s = ∆x/∆t √ 3 . In the lattice, ∆x and ∆t are maintained as 1 in most of simulations, hence the value of c 2 s = 1/3. The velocity of the molecules in each lattice point (α) is given according to the lattice structure and for of "D2Q9"lattice structure the velocity of the molecules are given below. The lattice weights of different lattice points (α) are assigned for "D2Q9" lattice structure is as shown below.

w α =        4/9, α = 0 1/9, α = 1, 2, 3, 4 1/36, α = 5, 6, 7, 8 (8.12)
The relaxation parameter τ is related to the dynamic shear viscosity ν.

ν = c 2 s τ - 1 2 (8.13)
The macroscopic local density and the velocity of the fluid can be computed for each lattice using the PDF's f a and molecular velocity e a values as shown in the equations below.

ρ i = N α=0 f α,i (8.14) u i = 1 ρ N α=0
f α,i e α (8.15)

Boundary conditions

In most of the CFD simulations, implementation of boundary conditions is one of the difficult parameters. Setting up the boundary values at the inlet, outlet and the wall has a huge effect on stability and the accuracy of the solution. In LBM, treatment of the boundary conditions is carried out in different methodology as the equations are based on the PDF and in classical CFD methods (using N-S equations) the treatment is based on macroscopic variables. Proper treatment of boundary condition in the PDF must be carried out to obtain the macroscopic boundary conditions of the fluid. There are different methods of implementing the boundary conditions in which the commonly used boundary conditions are discussed.

Solid Boundary conditions

The computation of the PDF is based on the particle streaming from the neighbouring location of a lattice structure, but the values of PDF are unknown from the direction of the wall. These unknown values must be determined by the boundary condition. There are 2 major types of boundary conditions, wet node and link wise boundary conditions. Difference between these boundary conditions is the location of the boundary on the lattice. Bounce Back conditions (BBC) are commonly known link wise boundary condition used to simulate the no-slip conditions on the boundaries. As the name suggests, when the particle reaches the wall or boundary the particle will reflect or bounce back to the previous location with the directions of the reversed velocity. This phenomenon also helps in conserving mass and momentum of the flow. These boundary conditions are easy to implement regardless of the orientation of the wall.

There are 2 different types of the variance of Bounce Back conditions, Full way bounce Back condition (FWBBC) and Halfway bounce back condition (HWBBC).

In FWBBC, the particles will finish the complete cycle of the stream and collide to implement the condition. In this condition, the particle will enter from fluid node to the solid node in the streaming stem and during the collision step the direction of the velocity is reversed and in the next streaming step the particle will go back to the fluid node from the solid node completing the cycle. FWBBC is first-order accurate which leads to decrease in the accuracy, as LBM is second-order accurate. Figure 8.5: Full way bounce back condition [START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF] In HWBBC the particle just takes half the cycle to implement the condition. In the streaming step, as soon as there is the presence of solid node the particle will return to the original location and the direction of the velocity of the particle is reversed. This method is faster compared to FWBBC and increasing the accuracy of the solution. The periodic boundary conditions can be easily implemented by defining outlet values of PDF to the inlet PDF values.

f α (inlet) = f α (outlet) ∀ α ∈ [1, d] (8.16)
Implementations of zero derivative conditions and extrapolation boundary conditions are easy to implement as well. Where i is the location of the node. Here, i = 1 is the inlet location and i = nx is the outlet of the domain. Implementation of zero derivative conditions are shown below:

f α (i = 1, j) = f α (i = 2, j) (8.17)

f α (i = nx, j) = f α (i = nx -1, j) (8.18)
Extrapolation boundary conditions can be computed as shown below:

f α (i = 1, j) = 2f α (i = 2, j) -f α (i = 3, j) (8.19)
f α (i = nx, j) = 2f α (i = nx -1, j) -f α (i = nx -2, j) (8.20)

Velocity and Pressure BC

In most of the simulations, the problem has to be simulated with a required velocity or pressure (In the form of density in LBM) value at the open boundaries. The BC developed by Zou & He [START_REF] Zou | On pressure and velocity boundary conditions for the lattice Boltzmann BGK model[END_REF] will be discussed here to determine the unknown values at the lattice points. Assuming the velocity boundary condition is given on one side of the scheme, for example in D2Q9 scheme f 0 , f 2 , f 3 , f 4 , f 6 , f 7 are known. The values of f 1 , f 5 , f 8 and ρ has to be determined. We know from LBM that the density and the velocity can be determined by: ρ( x, t) = f i e i (8.22) Using the above equations, we can rearrange the equation as:

f 1 + f 5 + f 8 = ρ -(f 0 + f 2 + f 4 + f 3 + f 6 + f 7 ) f 1 + f 5 + f 8 = ρu i + (f 3 + f 6 + f 7 ) f 5 -f 8 = ρu j -f 2 + f 4 -f 6 + f 7 (8.23)
By using these values, we can determine:

ρ = 1 1 -u i [(f 0 + f 2 + f 4 + 2 (f 3 + f 6 + f 7 )] (8.24)
To find all the variables we need another equation, Zou & He [START_REF] Zou | On pressure and velocity boundary conditions for the lattice Boltzmann BGK model[END_REF] assumed that the bounce back conditions hold good for the non-equilibrium part of the particle normal to the boundary. It can be assumed as:

f 1 -f eq 1 = f 3 -f eq 3 (8.25)
Here, f 1 can be determined using the f eq and to determine the other values use the following equations: Computation of the pressure inlet condition can be carried out in a similar methodology, where we have the value of density instead of the velocity.

f 1 = f 3 + 2 3 ρu j f 5 = f 7 - 1

Forcing scheme

The external force acting on the simulation has to be included to obtain the required results. The interaction force of Shan & Chan method is also considered as external forces which will be discussed in later sections. Proper implementation of the forcing term is important to reduce the error and increase the numerical stability and for the simulation of multicomponent multi-phase, as the forcing scheme plays a major role in particle interaction.

There are different forcing schemes which can be used to implement the force term. Shifted velocity method, Exact Difference Method (EDM) and Guo's Forcing method are few of the scheme which will be discussed here. Shifted velocity method is widely used and easy to implement as it is directly included in the equilibrium distribution function in the form of macroscopic velocity and also used in Shan & Chen [START_REF] Shan | Lattice Boltzmann model for simulating flows with multiple phases and components[END_REF] original equation. Shifted velocity method is derived from the momentum equation where the force is included in the bulk velocity as shown in the equation below.

Where F represents the external forces, τ is the relaxation time and u is the macroscopic velocity which is computed using the PDF in equation 8.15 and ρ is computed using the equation 8.14.

u eq = u + τ F∆t/ρ (8.27) The velocity u in the equilibrium distribution function is replaced by the bulk velocity u eq to obtain the collision equation with forcing term. The external force can be implemented through the shift in PDF known as Exact Difference Method (EDM) [START_REF] Kupershtokh | On equations of state in a lattice Boltzmann method[END_REF] which maintains the second-order accuracy of LBM and is independent of the relaxation time τ . EDM is implemented in the LBM equation as:

f α (x + ∆x, t + ∆t) = f α (x, t) - 1 τ f α (x, t) -f (eq) α (ρ, u) + ∆F α (8.29)
The external forces F are include in the term ∆F α with the help of equilibrium distribution function. The EDM method removes the error caused due to the density dependence on relaxation time. ∆F α = f (eq) α (ρ, u + F∆t/ρ) -f (eq) α (ρ, u) (8.30)

The Guo's Forcing scheme [START_REF] Guo | A coupled lattice BGK model for the Boussinesq equations[END_REF] is based on direct addition of a discrete form of force into the equation of distribution function and taking into account of the force contribution due to momentum flux and lattice effect.

∆F α = 1 - 1 2τ
w a e αu eq c 2 s + e α u eq c 4 s

• F∆t (8.31) The implementation of this scheme is similar to the EDM scheme where ∆F α is replace in the LBM equation 8.29. The external force is included in F, w a is the lattice weight,u eq is the bulk velocity which is computed using 8.27.

The computation of interaction force for Shan & Chen method is carried out using the original shift velocity method in this report.

Palabos

Palabos is an academic general-purpose computational fluid dynamics (CFD) code with a kernel based on the LBM. The code was developed by the Scientific and Parallel Computing group (SPC) at the university of Geneva. The libraries of Palabos code are based on C++ language and Palabos code has open-source AGPLv3 license [palabos.unige.ch]. The code is highly parallelized using the multi-block and dynamic objects approach.

These code libraries are used to utilize the LB equation and the boundary conditions, as these equations are well parallelized in Palabos [START_REF] Latt | Palabos: Parallel Lattice Boltzmann Solver[END_REF]. Another comparison criteria is the local evolution of velocity. Figure 8.9 shows a comparison between experimental results and CFD ones (WALE and Smagorinsky LES SGS turbulence models) of local velocities at L z = 0.817m and R = 0. We can see that the transition is delayed of t = 1s between experimental and WALE results. Indeed, the Smagornskgy model does not capture the transition that occurs at t ∼ 9.5s correctly. Nicoud showed that this model has a disadtanagte of poor modeling of the transition regime [START_REF] Nicoud | Subgrid-Scale Stress Modelling Based on the Square of the Velocity Gradient Tensor[END_REF][START_REF] Menter | Review of the shear-stress transport turbulence model experience from an industrial perspective[END_REF]. For the above reasons, the WALE SGS for LES turbulence model has been chosen for our numerical model.

Yan's expansion process comparison

Another comparison of numerical results of expansion process has been conducted. This time the comparison was conducted with external study from Yan et al. experimental study [START_REF] Yan | Experimental study of heat transfer enhancement in a liquid piston compressor/expander using porous media inserts[END_REF]. The geometry (cylinder) and expansion parameters are presented in table 8.4 D (m) L 0 (m) L pist,0 (m) t f (s) T 0 , T walls (K) P 0 (×10 8.10.A small piston speed is imposed at the beginning of the process, which increases afterwards and then decreases or is stable. This 202 CHAPTER 8. APPENDICES choice was justified by Yan et al. by its effectiveness for liquid piston porous media inserts, where it showed a temperature reduction [START_REF] Yan | Experimental study of heat transfer enhancement in a liquid piston compressor/expander using porous media inserts[END_REF]. Figure 8.10: Experimental piston speed profiles of the two presented test cases [START_REF] Yan | Experimental study of heat transfer enhancement in a liquid piston compressor/expander using porous media inserts[END_REF] The simulation of the relaxation in 8s shows very encouraging results. (Figure 8.11) However, there are two areas where the results deviate from the model:

• Between 1.4-2.7 V/V0. This difference can be explained by the extraction of experimental diffuse points, the small margin of error can be accepted as viable

• From 4.2 V/V0. The numerical simulation shows no damping. The proposed explanation for this difference comes from the lack of information for the end speed of the simulation.

The simulation of the expansion in 3s has the same remarks as the slow expansion in 8s. In addition, the temperature evolution profile shows a relatively difference evolution. The flow has a higher outlet velocity which varies more which makes it more difficult to implement numerically with the available data. [START_REF] Yan | Experimental study of heat transfer enhancement in a liquid piston compressor/expander using porous media inserts[END_REF] and numerical ones of two expansion processes Similar air's temperature evolution profile as the one presented in chapter 5 can be seen here, especially for the slow (8s) expansion. While some minor gaps can still be seen in Figure 8.11 between numerical and experimental results which are mainly due to lack of some information of the experimental study, the comparison enables us to confirm the numerical model. Une thèse expérimentale sur le sujet a été réalisée entre Segula et le laboratoire GEPEA à Nantes [START_REF] Neu | Etude expérimentale et modélisation de la compression quasi isotherme d'air pour le stockage d'énergie en mer[END_REF]. Des études numériques ont suivi, et ont porté sur la modélisation du processus de compression dans un piston liquide avec la méthode VOF basée sur la résolution des équations de Navier-Stokes [START_REF] Mehdi | Méthode VOF pour la simulation numérique de l'écoulement de l'air comprimé et du transfert thermique associé dans un piston liquide[END_REF]. Ces études ont démontré qu'il était possible de modéliser l'écoulement et l'évolution de la température induits par la compression d'air dans un piston liquide. La méthode LBM (Lattice-Boltzmann Method) a été choisie initialement pour mener cette étude du fait de son caractère linéaire qui la rendait performante pour les gros calculs. Toutefois, la méthode présentait plusieurs verrous scientifiques qui se sont révélés très difficiles à résoudre. Il a donc été choisi de poursuivre cette étude avec une approche plus classique : la résolution des équations de Navier-Stokes avec la méthode des volumes finis sur OpenFoam.

Synthèse de la thèse

Objectifs de la thèse

Le résultat principal attendu de ce travail de recherche consiste à réaliser une étude globale des processus de compression/détente par piston liquide. Ceci consiste à réaliser des simulations numériques afin de comprendre l'évolution de l'écoulement, des champs de vitesse, de température et de pression. Des résultats théoriques sur l'efficacité énergétique et exergétique d'un piston liquide sur un cycle complet (compression-détente) sont notamment attendus, afin d'orienter de manière précise les études numériques. Par ailleurs, l'exploitation physique des résultats devra permettre une meilleure compréhension, une qualification et une quantification des phénomènes physiques dus à l'écoulement et l'échange thermique associé dans le piston liquide. Aussi, il est attendu l'exploration de solutions permettant l'amélioration du rendement énergétique du piston liquide. Il est aussi prévu, en perspective de ce travail, de réaliser une étude d'optimisation fine du rendement énergétique du système REMORA.

Verrous scientifiques

Á travers l'étude expérimentale précédente et l'étude bibliographique menée, il a été mis en évidence les différents phénomènes associés à la compression d'air par piston liquide. Ceuxci présentent des verrous scientifiques à lever. En fait, l'écoulement dans la chambre de compression est un écoulement multiphasique (gaz-liquide) et compressible auquel sont associés des transferts thermiques (convection mixte, conduction [START_REF] Odukomaiya | Thermal analysis of near-isothermal compressed gas energy storage system[END_REF]). D'autres phénomènes ont été observés lors de la compression d'air par piston liquide. Il a été constaté qu'au cours de la compression, un changement de régime convectif a lieu, ainsi l'occurrence du passage du régime convectif laminaire à une régime turbulent a été caractérisée par une formule empirique à partir des résultats expérimentaux [START_REF] Neu | Etude expérimentale et modélisation de la compression quasi isotherme d'air pour le stockage d'énergie en mer[END_REF] et a été confirmé numériquement dans l'étude menée précédemment [START_REF] Mehdi | Méthode VOF pour la simulation numérique de l'écoulement de l'air comprimé et du transfert thermique associé dans un piston liquide[END_REF]. Par rapport à la vitesse de la colonne d'eau, la vitesse d'écoulement de l'air varie de 1 à 10. Ceci peut induire un changement de régime d'écoulement. Expérimentalement, il a été difficile d'évaluer le champs de vitesse de l'air. Dans l'étude bibliographique, certaines études se sont basées sur la vitesse d'entrée de l'eau, ce qui ne permet pas de calculer correctement le champs de vitesse de l'air. Numériquement, la simulation d'un écoulement multiphasique à grand rapport de densités/grand Reynolds constitue un défi majeur. De plus la résolution de la thermique de l'écoulement rend cette approche inédite. D'autres approches possibles, notamment la simulation d'un piston solide avec la méthode de la frontière immergée (IBM) ou d'une frontière mobile nécessitent des adaptations spécifiques s'agissant d'un domaine évoluant dans le temps et d'un couplage avec le transfert de chaleur.

Chapitre 2 :

Ce chapitre présente une analyse de la technologie du Piston Liquide (LP) pour une application dans les systèmes CAES (Compressed Air for Energy Storage). Divers aspects sont abordés, notamment l'état des lieux des projets LP dans le monde et la tendance du développement. Elle met en évidence la complexité du couplage entre l'écoulement des fluides et le transfert de chaleur pendant les opérations de compression/détente, ainsi que les différentes actions proposées et mises en oeuvre pour améliorer le transfert de chaleur à l'intérieur de la colonne du piston. Alors que la plupart des études tentent de caractériser le transfert de chaleur à l'aide de corrélations empiriques (préexistantes et développées), l'évolution de l'écoulement n'est toujours pas bien comprise. Comme il est difficile d'avoir une visualisation complète de l'écoulement de manière expérimentale, la modélisation numérique est choisie pour étudier l'écoulement et le transfert de chaleur à l'intérieur d'un LP.

Chapitre 3 :

Dans ce chapitre la méthode Lattice Boltzmann (LBM) est présentée. Elle a été initialement choisie pour notre étude du fait de ses différents avantages liés à la rapidité de calcul et sa capacité à être couplée avec des algorithmes d'optimisation. Le principal défi de la modélisation de l'écoulement à l'intérieur d'un LP à l'aide de la LBM est la gestion de l'interface liquide/gaz pour un écoulement multiphasique avec un rapport de densité élevé. Bien que plusieurs méthodes de modélisation multiphasique existent déjà, elles sont limitées à la plage du rapport de densité (1-10), alors que ce rapport est de 1000 pour un écoulement eau/air. Différentes approches ont été mises en oeuvre et testées. Tout d'abord, l'approche modifiée de Shan-Chen, pour les écoulements multiphasiques à un composant et à plusieurs composants, a été utilisée pour modéliser des écoulements à rapports de densité élevés. L'approche souffre de courants parasites, en particulier pour les flux MCMP (MultiPhase Multi-Component), ce qui entraîne des schémas instables. Une deuxième approche "Free-Energy" est implémentée dans un code interne en utilisant deux modèles "Lee-Lin" [START_REF] Lee | Lattice Boltzmann simulations of micron-scale drop impact on dry surfaces[END_REF] et "Inamuro" [START_REF] Inamuro | Lattice Boltzmann methods for viscous fluid flows and for two-phase fluid flows[END_REF], le premier a permis de simuler une goutte statique à un rapport de densité de 1000 mais n'a pas permis de modéliser un cas dynamique, tandis que le deuxième modèle, qui a un coût de calcul très élevé a permis une simulation de collision de deux gouttes avec un Re élevé et un rapport de densité moyen [START_REF] Higuera | Boltzmann approach to lattice gas simulations[END_REF]. Enfin, un couplage de la LBM avec la méthode VOF appliqué à un écoulement diphasique est testé. Cette nouvelle approche inédite a permis la simulation d'écoulements multiphasiques à grands rapports de densité/viscosité et à hauts régimes d'écoulement. Toutefois, elle souffre d'un manque de conservation de la masse. Une décision a ensuite été prise de choisir une approche plus classique qui est présentée dans les chapitres suivants.

Chapitre 4 :

La modélisation de la compression d'air dans un LP est présentée dans ce chapitre. Elle consiste à la résolution des équations de Navier-Stokes via la méthode des volumes finis (FVM) et la méthode de volume du fluide (VOF) pour le suivi des interfaces.

Les résultats obtenus via le modèle numérique ont été comparés et ensuite validés par comparaison avec les résultats expérimentaux réalisés dans le cadre du projet REMORA. La nature 3D de l'écoulement est démontrée. En effet, celui-ci passe par quatre phases, la plus importante étant l'établissement d'écoulement axisymétrique structuré. Nous avons constaté que pendant cette période la vitesse de l'air varie de [-5U pist , 12U pist ], ce qui fait que le nombre de Reynolds calculé à partir des paramètres initiaux du gaz et basé sur U pist n'est pas un nombre adimensionnel fiable pour caractériser l'évolution de l'écoulement. Ensuite, une transition se produit. Elle perturbe celui-ci et le rend totalement agité. La transition est d'abord observée sur la composante de la vitesse dans la direction de la compression. Elle se produit là où les gradients de vitesse sont les plus importants. Le cisaillement entre les courants positifs et négatifs pourrait être la raison principale de cette perturbation.

Une étude paramétrique a été menée pour mesurer les effets des différents paramètres de compression, notamment la vitesse du piston et la température des parois sur l'évolution de l'écoulement et l'efficacité. Une compression plus rapide entraîne une évolution de la température plus importante et donc un rendement plus faible (chute de η c = 91% pour U pist = 0.033 m/s à η c = 82% pour U pist = 0.1 m/s ). La température des parois (milieu extérieur ambiant), est un paramètre clé pour l'évolution de la température de l'air. Une température de paroi plus élevée entraîne une température de l'air plus élevée et donc une efficacité de compression plus faible.

Chapitre 5 :

Dans le chapitre 5, le processus de détente est étudié, la même analyse que celle du chapitre 4 est menée. Les résultats montrent une similitude entre l'évolution de l'écoulement en détente et celle observée au cours de la compression. La principale différence est l'évolution de la température, qui diminue au début puis augmente après la transition. Une étude paramétrique est menée et montre l'effet important du temps de détente qui pourrait favoriser la convection et donc la réduction de la température et l'amélioration du rendement (jusqu'à η = 94%). Un cycle compression-refroidissement isochore-détente est également étudié. Nous avons montré que la détente polytropique entraîne également un processus de réchauffement de l'air, et donc que le cycle est plus proche d'un cycle quasi-isotherme puisque l'indice polytropique est de n = 1.012 pour la détente. L'efficacité du cycle est également très élevée par rapport à un cycle adiabatique, puisqu'elle atteint η cycle = 88.5%.

Chapitre 6 :

Dans le chapitre 6, deux approches sont proposées pour modéliser le transfert de chaleur conjugué entre la paroi solide et le domaine fluide. Les deux approches explicites donnent des résultats similaires. Si peu de différences sont observées avec le modèle précédent, les effets du matériau solide (Plexiglass, cuivre et verre borosilicate) ont été mis en évidence. Nous avons constaté que le matériau ayant la plus faible inertie thermique (effusivité) a un impact significatif sur l'évolution de la température de l'air (∆T = 3K).

Conclusion et perspectives

Le piston liquide offre des efficacités de compression/détente élevées (91% à 95 %), avec une meilleure convenance aux applications sous-marines, il constitue une technologie fiable pour l'implémentation dans systèmes de stockage d'énergie par air comprimé. Nous avons montré dans ce travail que nous pouvons obtenir des processus polytropiques qui peuvent se 209 CHAPTER 8. APPENDICES rapprocher de quasi-isothermes en optimisant les paramètres du processus (profil de vitesse du piston, température du milieu externe, matériaux de paroi...).

Une meilleure visualisation de l'écoulement à l'intérieur d'un LP est maintenant possible grâce au modèle numérique développé. Tandis que, dans cette thèse, une description fine est donnée pour l'établissement de la structure de l'écoulement, les origines de celle-ci restent encore à déterminer. Sa perturbation reste un élément clé pour des études ultérieures. Une analyse physique plus approfondie est encore nécessaire pour mieux comprendre l'évolution de l'écoulement. Les échanges de chaleur qui se produisent à la paroi doivent être quantifiés par le calcul du nombre de Nusselt ou du coefficient de transfert de chaleur. Des études complémentaires sur ce sujet permettraient de comprendre l'effet de chaque paramètre. Comme il a été montré dans cette thèse, le nombre de Reynolds ne permet pas de caractériser l'écoulement. De nouvelles lois permettant l'estimation du transfert thermique et la prédiction de la transition permettraient de réduire et de simplifier le modèle numérique, qui reste lourd en termes du coût de calcul.

Les modèles développés pour la modélisation de la paroi solide offrent la possibilité d'étudier des géométries plus complexes (forme ou différentes insertions) et leurs effets sur l'écoulement. D'autres phénomènes qui n'ont pas été modélisés comme le transfert de masse ou le changement de phase peuvent aussi être étudiés. En effet, la température locale observée en fin de compression et surtout lors de la détente avant la transition peut entraîner un changement de phase de l'air.

Une étude plus approfondie des effets de la géométrie du piston (forme, L/D), des matériaux de la paroi solide et des paramètres de compression/expansion est nécessaire pour déterminer les effets précis de chacun. Ceci peut être réalisé par une optimisation de la topologie qui permettrait d'augmenter l'efficacité des processus de compression/détente en utilisant un LP pour les applications CAES.

Enfin, des simulations de LP à l'échelle industrielle et à des rapports de pression élevés, avec des cycles multi-étages et successifs, permettront d'optimiser la configuration du LP pour les applications CAES. Le principe de piston liquide présente plusieurs avantages par rapport à un piston solide pour une utilisation dans un système de stockage d'énergie sous forme d'air comprimé. Dans ce contexte, le projet REMORA développé par Segula Technologies où le piston liquide est utilisé pour réaliser les processus compression (stockage) et détente (restitution) d'air comprimé dans l'objectif de stocker l'énergie électrique produite. Les travaux sur la compréhension des phénomènes physiques liés à la compression/détente par piston liquide sont absents de la littérature scientifique. L'objectif de cette thèse est de modéliser numériquement la compression / détente par piston liquide. Initialement, la méthode LBM (Lattice Boltzmann Method) est utilisée mais n'a pu donner les résultats souhaités en raison de verrous scientifiques majeurs.

Le choix a ensuite porté sur la résolution des équations de Navier-Stokes et la méthode VOF (Volume of Fluid) pour le suivi de l'interface. Les résultats du modèle numérique ont pu être comparés avec des études expérimentales ce qui a permis de les valider. La nature 3D de l'écoulement et la transition (du régime d'écoulement et de transferts thermiques) qui se produit au cours de la compression et de la détente ont pu être étudiés. Un modèle a ensuite développé pour prendre en compte l'effet de l'évolution de la température dans les parois sur l'écoulement et la température de l'air. Une étude paramétrique de ces processus a aussi été réalisée et a montré l'effet de certains paramètres (vitesse du piston, température et matériau de la paroi) sur l'évolution de la température et le rendement.

Title : Investigation of a compressed air energy storage system: Flow and heat transfer numerical modeling in a liquid piston Keywords : Liquid piston, CFD, VOF Method, Heat transfer, compression/expansion Abstract : The liquid piston principle has several advantages over the solid for a compressed air energy storage application. In this context, the REMORA project developed by Segula Technologies where the liquid piston is used to carry out the processes of compression (storage) and expansion (restitution) of compressed air with the objective of storing eletrical energy. Proper understanding of the physical phenomena related to the compression/expansion using a liquid piston are missing from the scientific literature. The objective of this thesis is to model numerically the compression / expansion processes using a liquid piston. Initially, the LBM (Lattice Boltzmann Method) is used but could not give the expected results because of major scientific barriers. The choice was then made to solve the Navier-Stokes equations coupled to VOF (Volume of Fluid) method to track the water/air interface. The results of the numerical model have been compared with experimental ones which enabled to validate them. The 3D nature of the flow and the transition (of the flow and heat transfer regime) that occurs during compression and expansion have been highlighted. A model was then developed to model the effect of the temperature evolution of the solid walls on the flow. A parametric study of these processes was also carried out and showed the effect of some parameters (piston speed, ambient temperature, solid material) on temperature evolution and process efficiency.

Figure 2 . 1 :

 21 Figure 2.1: Trend of published documents (journal articles, conference proceedings, book chapters, patents) with "Liquid Piston" in title or abstract for in Dimensions [27] (August 2021)

Figure 2 . 2 :

 22 Figure 2.2: Operation of LP during the compression (a) and the expansion (b) processes

Figure 2 . 3 :

 23 Figure 2.3: P -V and T -S diagrams for air compression and expansion processes using different cycles (Otto, Joule, Polytropic and Isothermal).Adapted from[START_REF] Lemofouet | Investigation and optimisation of hybrid electricity storage systems based on compressed air and supercapacitors[END_REF]

Figure 2 . 4 :

 24 Figure 2.4: P -V curves showing compression (a) and expansion (b) trajectories for a LPin an open-accumulator system CAES[START_REF] Yan | Experimental study of heat transfer enhancement in a liquid piston compressor/expander using porous media inserts[END_REF]

Figure 2 . 5 :

 25 Figure 2.5: Evolution of air flow pattern along with the advancement of the water during the compression operation [102]

Figure 2 . 6 :

 26 Figure 2.6: Heat transfer modeling for LP compression/expansion, including the convectionbetween the gas and liquid through the interface, the convection at the inner and outer walls of the piston column and the conduction through the wall. Adapted from[START_REF] Odukomaiya | Thermal analysis of near-isothermal compressed gas energy storage system[END_REF]

Figure 2 . 7 :

 27 Figure 2.7: Investigated alternative liquid spray configurations of GLIDES system [109]

Figure 2 . 8 :

 28 Figure 2.8: Efficiency increase of Liquid Piston by the liquid spray injection. (a) Compression Ratio (CR)=10 and variable Mass Loading (ML); (b) variable CR, MeanDroplet Diameter (MDD)=20 µm). Adapted from[START_REF] Qin | Liquid piston compression efficiency with droplet heat transfer[END_REF]

Figure 2 . 9 :

 29 Figure 2.9: Five types of porous inserts used in compression/expansion experiments [174], including 3 Acrylonitrile Butadiene Styrene (ABS) interrupted-plates with 2.5 mm, 5 mm and 10 mm spacing between plates and 2 aluminium foams with 10 and 40 pores per inch (PPI)

Figure 2 . 10 :

 210 Figure 2.10: Comparison of compression/expansion efficiencies with or without the porous media inserts [174]. (a). Compression efficiency vs. power density (CR=10); (b). Expansion efficiency vs. power density (expansion ratio=6)

Figure 2 . 12 :

 212 Figure 2.12: Honeycomb insert configuration for HTE in LP. (a). thin tubes [179]; (b).honeycomb matrix[START_REF] Zhang | Numerical analysis of heat exchangers used in a liquid piston compressor using a one-dimensional model with an embedded two-dimensional submodel[END_REF]

Figure 2 . 13 :

 213 Figure 2.13: Hollow spheres concept for LP. Full LP domain ans Hollow spheres layer at the water/air interface

Figure 2 . 14 :

 214 Figure 2.14: Experimental results on the effect of hollow spheres with different materialsand sizes for CR=2. Adapted from[START_REF] Ranganath Ramakrishnan | Temperature abatement using hollow spheres in liquid piston compressor for Ocean Compressed Air Energy Storage system[END_REF]

Figure 2 . 16 :

 216 Figure 2.16: Wire mesh insert for LP [116]. (a). schematic model; (b) photo of aluminum and copper wire mesh in Archimedean spiral form

Figure 2 . 17 :Figure 2 . 18 :

 217218 Figure 2.17: Wire mesh insert on the temperature evolution of LP [116]. (a) Effect of mesh materials; (b) Effect of compression speed

Figure 2 .

 2 Figure 2.19: CFD-based shape optimization of the LP column [177]. (a). Effect of single shape parameter on the efficiency improvement; (b). Optimal geometry of the LP column and the temperature field at the end of compression

22 )uFigure 3 . 1 :

 2231 Figure 3.1: Channel flow boundary condition

Figure 3 . 2 :

 32 Figure 3.2: Comparison of analytical solution and the LB numerical simulation

  (a) Velocity (top) and pressure (bottom) contours at initial time step (b) velocity (top) and pressure (bottom) contours after few time steps

Figure 3 . 3 :

 33 Figure 3.3: Velocity and pressure contours at different time steps

  (a) velocity (top) and pressure (bottom) contours at time step where the flow reaches the end of the channel (b) Velocity (top) and pressure (bottom) contours as the flow develops

Figure 3 . 4 :

 34 Figure 3.4: Velocity and pressure contours at different time steps

Figure 3 . 6 :

 36 Figure 3.6: Static drop SCMP simulation schematic

Figure 3 . 7 :

 37 Figure 3.7: Density distribution contour of the simulation for T r value of 0.8

7 Figure 3 . 8 :

 738 Figure 3.8: Density distribution for different value of T r

  Figure 3.9: Spurious Velocity at the interface of 2 phases at reduction temperature T r = 0.8

7 Figure 3 . 10 :

 7310 Figure 3.10: Spurious Velocity at the interface of 2 phases at different value of reduction temperature T r

Figure 3 . 11 :

 311 Figure 3.11: Density variation for different value of β

Figure 3 . 12 :

 312 Figure 3.12: Effect of β on velocity and pressure.

Figure 3 . 13 :

 313 Figure 3.13: Static drop MCMP simulation boundary condition

Figure 3 . 14 :

 314 Figure 3.14: Density distribution of air in MCMP

Figure 3 . 16 :

 316 Figure 3.16: Plot of density variation of water and air.

  (a) Velocity variation of air component. (b) Velocity variation of water.

Figure 3 . 17 :

 317 Figure 3.17: Velocity variation at the end of the simulation.

Figure 3 .

 3 Figure 3.18: Density distribution of air with condensation of water droplets in MCMP

  (a) Negative density of water component at the beginning. (b) Negative density of water component at the end of the simulation

Figure 3 . 19 :

 319 Figure 3.19: Negative density of water component.

Figure 3 .

 3 Figure 3.20: Total density of the flow simulation.



  Figure 3.21: Channel flow of MCMP simulation boundary condition where air is pushed by water

  (a) Density distribution before applying the inlet and outlet condition. (b) Density distribution after applying the inlet and outlet condition.

Figure 3 . 22 :

 322 Figure 3.22: Density distribution before and after applying the actual condition of channel flow

  (a) Velocity distribution of air. (b) Velocity distribution of water.

Figure 3 . 23 :

 323 Figure 3.23: Velocity distribution of air and water after applying inlet and outlet BC.

  (a) Density distribution when the flow reaches middle of channel. (b) Density distribution when the flow reaches almost end of channel.

Figure 3 . 24 :Figure 3 . 25 :

 324325 Figure 3.24: Density distribution as the flow advances in the channel

Figure 3 . 26 :

 326 Figure 3.26: Flowchart of the different steps of the numerical implementation of Inamuro's model

Figure 3 .

 3 Figure 3.27: Caption

Figure 3 . 28 :Figure 3 . 29 :

 328329 Figure 3.28: Time evolution of droplet shape

5 Figure 3 . 30 :

 5330 Figure 3.30: U x and u y velocity fields of liquid droplet inside a gas domain

•Figure 3 . 31 :Figure 3 . 32 :

 331332 Figure 3.31: Flowchart representation of the key steps of the LBM-VOF coupling in Palabos

Figure 3 . 33 :

 333 Figure 3.33: Geometry of the tested cases two-phase flow and compression flow

  (a) U z component of the velocity field after flow establishment in a multiphase flow without compression (Open-Top) (b) Volume Fraction distribution in the computational domain, representing the two phases. α = 1is the liquid phase and α = 0 is the gas phase

Figure 3 . 34 :

 334 Figure 3.34: Velocity and pressure fields of the "Open-top" test

Figure 3 . 36 :Figure 3 . 37 :

 336337 Figure 3.36: Velocity and pressure fields of the compression1 test

Figure 4 . 1 :

 41 Figure 4.1: Energy storage processes (compression, compressed air ejection to storage tank and then water ejection)

Figure 4 . 2 :

 42 Figure 4.2: Classification of numerical methods for two-phase flows.Methods that use a sharp-interface approach have a white background while diffuse-interface approaches have a gray background[START_REF] Mirjalili | Interface-capturing methods for two-phase flows : An overview and recent developments[END_REF]

(4.6) 103 CHAPTER 4 .

 1034 MODELING OF THE COMPRESSION PROCESS USING A LIQUID PISTON where e k = |U| 2 2

Figure 4 . 3 : 4 :

 434 Figure 4.3: Experiment setup, as realized by Neu and Subrenat [104]

Figure 4 . 5 :

 45 Figure 4.5: A comparison of air's average temperature evolution through time with three 2D meshes

  (a) Radial view of a hexahedral Mesh (O-grid) with local refining near the wall (b) Radial view of prisms mesh

Figure 4 . 6 :

 46 Figure 4.6: Radial view the two types of meshes used for the mesh convergence study

Figure 4 . 7 :

 47 Figure 4.7: A comparison of air's average temperature evolution through time with several used meshes

  (a) 2D model U z velocity field at t = 2.24s (b) 3D model U z velocity field at t = 2.24s

Figure 4 . 8 :

 48 Figure 4.8: Comparison between the U z velocity fields of the 2D and 3D models at t = 2.24s

Figure 4 . 9 :

 49 Figure 4.9: Comparison between the U z velocity fields of the 2D and 3D models at t = 6.7s

  (a) 2D model temperature field at t = 8.92s (b) 3D model temperature field at t = 8.92s

Figure 4 . 10 :

 410 Figure 4.10: Comparison between the temperature fields of the 2D and 3D models at t = 8.92s

Figure 4 . 11 :

 411 Figure 4.11: Comparison on the evolution of average air temperature over compression time between numerical results and experimental data [104].

Figure 4 . 12 :

 412 Figure 4.12: Comparison on the evolution of average air pressure over compression time between numerical results and experimental data [104].

( a )

 a Experimental axial velocity field[START_REF] Neu | Experimental investigation of internal air flow during slow piston compression into isothermal compressed air energy storage[END_REF] (b) Numerical axial velocity field

Figure 4 . 13 :Figure 4 . 14 :

 413414 Figure 4.13: Comparison of the axial velocity field between experimental [104] and CFD results at t = 3.15s

  Figure 4.15: The U z (m/s) velocity component profile between t = 0s and t = 5s

Figure 4 . 16 :

 416 Figure 4.16: The U z (m/s) velocity component profile and streamlines between t = 0s and t = 5s

, 4 .

 4 20a). It then moves along the walls until it returns to the area near the piston's top, where it is again carried into the central flow.

Figure 4 . 17 :

 417 Figure 4.17: The U z (m/s) velocity component profile at t = 9.5s, while the transition occurs

  , 4.20b).

Figure 4 . 18 :

 418 Figure 4.18: The U z (m/s) velocity component profile after the transition (at t = 12.5s)where the flow is chaotic

Figure 4 . 19 :

 419 Figure 4.19: The U z (m/s) velocity component profile at the end of the compression (at t = 20.5s

( a )

 a Temperature field at t = 4.25s (b) Temperature field at t = 9.4s (c) Temperature field at t = 14.25s (d) Temperature field at t = 20s

Figure 4 . 20 :

 420 Figure 4.20: Evolution of temperature field through compression time

Figure 4 . 21 :

 421 Figure 4.21: Evolution of velocity U z profile at L z = 0.650m

Figure 4 . 22 :

 422 Figure 4.22: Evolution of temperature profile at L z = 0.650m

Figure 4 . 23 :

 423 Figure 4.23: Longitudinal L z and radial R local positions presented

Figure 4 . 24 : 4 . 5 . 4 . 2

 4244542 Figure 4.24: Local velocity component U z through time at L z = 0.216m

Figure 4 . 25 : 4 . 5 . 4 . 3

 4254543 Figure 4.25: Local velocity component U z through time at L z =0.443m

Figure 4 . 26 :

 426 Figure 4.26: Local velocity U z through time at L z =0.650m

Figure 4 . 27 :

 427 Figure 4.27: Local velocity U x through time at L z =0.650m

Figure 4 . 29 :

 429 Figure 4.29: Comparison between experimental results and CFD results of local velocities at L z = 0.817m

Figure 4 . 30 :

 430 Figure 4.30: Radial view of velocity U z during the first apparition of the structure disruption at t = 8.2s. With adapted scale (m/s) for each position

Figure 4 . 31 :

 431 Figure 4.31: Radial sections (L z = 0.4m, L z = 0.6m, L z = 0.7m, L z = 0.8m) of the U z velocity at different moments

Figure 4 .

 4 [START_REF]EIA projects nearly 50% increase in world energy usage by 2050, led by growth[END_REF] shows that for the L z = 0.216m and for 133CHAPTER 4. MODELING OF THE COMPRESSION PROCESS USING A LIQUID PISTONL z = 0.443m the flow rate increases linearly with the piston headway even when the transition occurs to meet the water flow rate when the interface reached L z . At L z = 0.650m and L z = 0.817m, in stage (B) (Figure4.11, the numerically calculated flow rate is different than the one calculated by the assumption. Although both results show a similar profile. Closer the the cylinder head L z = 0.883m, the flow undergoes a recirculation and a dispersion in the radial directions (Figure4.15). The proposed correlation does not take into account the negative currents, especially after the transition occurs t ∼ 10s and close to the piston top. The model evolves linearly, while the transition changes the flow structure.

Figure 4 . 32 :

 432 Figure 4.32: Flow rate evolution during time at different positions

Figure 4 .

 4 Figure 4.33: P-V diagram of the modeled compression process using a liquid piston

Figure 4 . 34 :

 434 Figure 4.34: Comparison between air's temperature evolution for two different U pist = 0.033 m/s; 0.1 m/s

  4 are used, while only T walls changes. While starting from the same initial temperature, the air's Configuration cold reference hot Walls temperature T walls (

Figure 4 . 36 :

 436 Figure 4.36: Comparison between air's temperature evolution for three different T walls = 288K; 293K; 298K

Figure 4 . 38 :

 438 Figure 4.38: Evolution of compression efficiency depending on the wall temperature using the three different models

Figure 5 . 1 :

 51 Figure 5.1: Energy production or recovering processes (water injection, air injection then expansion)

8 ×Figure 5 . 2 :

 852 Figure 5.2: Liquid Piston geometry and boundary conditions of the expansion process

Figure 5 . 4 :

 54 Figure 5.4: Comparison between numerical results and experimental data of air's pressure evolution during expansion process

Figure 5 . 5 :

 55 Figure 5.5: Comparison of the evolution of average air temperature over expansion time between numerical results and experimental data

  5):• (A) High quasi-linear decrease of air's average temperature from the initial time and until t ∼ 1.2s.• (B) Slower decrease of air's average temperature 1.2 < t < 3.2s • (C) Constant air's average temperature from 3.2 < t < 4.2s • (D) Temperature increases until the end of the expansion At the beginning of the expansion (Zone A), a symmetrical structure is set up. The latter consists of a fast central flow, surrounded by a faster recirculation zone close to the wall. The temperature drops rapidly, except in the near-wall zones (Figure 5.6).

Figure 5 . 6 :

 56 Figure 5.6: Temperature (top) and velocity U z (bottom) fields at t = 1s

Figure 5 . 7 :

 57 Figure 5.7: Temperature (top) and velocity U z (bottom) fields at t = 2.4s

Figure 5 . 8 :Figure 5 . 9 :

 5859 Figure 5.8: Temperature (top) and velocity U z (bottom) fields at t = 4s

Figure 5 . 11 :Figure 5 . 12 :Figure 5 . 13 :

 511512513 Figure 5.11: Evolution of local velocity U z through time at L z = 0.817m

Figure 5 .

 5 Figure 5.14: P-V diagram of the modeled expansion process using a liquid piston

Figure 5 . 18 :

 518 Figure 5.18: Comparison of Expansion Efficiencies for the studied configurations using Yan model

Figure 5 .

 5 Figure 5.19: Compression-storage-Expansion cycle for CAES applications

Figure 5 . 23 :

 523 Figure 5.23: Comparison between experimental and numerical results of the pressure evolution during a compression-cooling-expansion cycle

Figure 5 .

 5 Figure 5.24: P-V diagram of the compression-cooling-expansion cycle

Figure 5 .Figure 5 . 27 :

 5527 Figure 5.25: Temperature (top) and velocity U z (bottom) fields at t cool = 0.2s

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

 616263 Figure 6.1: Longitudinal view of a liquid piston with solid walls

Figure 6 . 5 :

 65 Figure 6.5: Conformal meshing of the two domains Solid (in red) and fluid (in green)

Figure 6 . 6 :

 66 Figure 6.6: Evolution of air's temperature during compression process with CHT using the two models and isothermal walls

Figure 6 . 7 :

 67 Figure 6.7: Comparison of temperature fields in the solid domain between the two proposed models during the compression process

Figure 6 . 8 :

 68 Figure 6.8: Comparison of U z velocity fields between the initial model (isothermal) and the two new proposed models during the compression process

  (a) t = 9s (b) t = t f = 25s

Figure 6 . 9 : 1 Figure 6

 6916 Figure 6.9: Comparison of temperature fields in the solid domain between external wall Boundary conditions (isothermal T ext = 293K, adiabatic, and convective heat flux). All the simulations used the model 1

Figure 6 . 10 :

 610 Figure 6.10: Evolution of air's temperature during compression process with three different solid wall materials

Figure 6 . 11 :

 611 Figure 6.11: Evolution of air's temperature during expansion process with the initial model (isothermal walls), the CHT model 1 and the experimental results

Figure 8 . 1 :

 81 Figure 8.1: Molecular streaming and collision representation

Figure 8 . 2 :

 82 Figure 8.2: Lattice grid discritization for simulation using LBM

we

  i e iα e iβ = e 2 s δ αβ (8.5) (d) i w i e iα e iβ e iγ = iα e iβ e iγ e iδ = c 4 s (δ αβ δ γδ + δ αγ δ βδ + δ αδ δ βγ ) (8.7) (f ) i w i e iα e iβ e iγ e iδ e i = 0 (8.8) There are models which require additional properties like thermal LBM which needs the velocities interacting with next to nearest neighbours. The figure 8.3 represents the lattice structure of D2Q9 where f represents the PDFs at different points of the lattice.

Figure 8 . 3 :

 83 Figure 8.3: D2Q9 Lattice structure

8.0.0. 3

 3 Link wise Boundary conditionsLink wise Boundary conditions are computed considering the location of the boundary next to the last fluid node, which is on the solid part of the lattice as shown in the figure8.4, hence the flow cannot be evaluated by PDF at that location. Different methodology has to be chosen to compute to values at he boundary.

Figure 8 . 4 :

 84 Figure 8.4: Link wise boundary condition[START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF]

Figure 8 . 6 :

 86 Figure 8.6: Half way bounce back condition[START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF]

Figure 8 . 7 :

 87 Figure 8.7: Wet node boundary condition[START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF]

3 e a • u eq c 2

 32 f eq a (x, t) = w a ρ(x) 1 +

Figure 8 . 8 :

 88 Figure 8.8: Air's average temperature evolution : A comparison between Comparison between experimental results and CFD ones (WALE and Smagorinsky LES SGS turbulence models)

Figure 8 . 9 :

 89 Figure 8.9: Comparison between experimental results and CFD ones (WALE and Smagorinsky LES SGS turbulence models) of local velocities at L z = 0.817m

( a )Figure 8 . 11 :

 a811 Figure 8.11: Normalized air's average temperature evolution comparison betweenexperimental results[START_REF] Yan | Experimental study of heat transfer enhancement in a liquid piston compressor/expander using porous media inserts[END_REF] and numerical ones of two expansion processes

  Ce travail de recherche s'intègre dans le cadre d'une thèse CIFRE entre Segula Technologies et le LTeN (Laboratoire de Thermique et Énergie de Nantes). Cette thèse s'inscrit dans le cadre du projet REMORA (Système de Stockage d'énergies sous forme d'air comprimé) et vise à étudier numériquement avec la méthode du Volume du Fluide (VOF) l'écoulement et les transferts thermiques engendrés par la compression quasi-isotherme d'air par piston liquide (LP).

Titre:

  Etude d'un système de stockage d'énergie par compression d'air : Modélisation numérique de l'écoulement et du transfert thermique dans un piston liquide Mots clés : Piston Liquide, CFD, Méthode VOF, transferts thermiques, compression/détente Résumé :

  defined compression/expansion formulas (Eqs. 2.6 CHAPTER 2. REVIEW ON LIQUID PISTON TECHNOLOGY FOR COMPRESSED AIR ENERGY STORAGE and 2.7) only related to CR and n the polytropic index :

  2.7: without liquid spray (base configuration), liquid spray with or without heat regeneration. The tests conducted under the base configuration showed a temperature rise up to 40 K during CHAPTER 2. REVIEW ON LIQUID PISTON TECHNOLOGY FOR COMPRESSED AIR ENERGY STORAGE

		100%							100%
	(%)	95%						(%)	95%
	Efficiency η C	80% 85% 90%						Efficiency η C	80% 85% 90%
	Compression	65% 70% 75%					Pre-mixed Injection Direct Injection	Compression	65% 70% 75%					Pre-mixed Injection Direct Injection
		60%							60%
		0	1	2	3	4	5	6	6	7	8	9	10	11	12
				Mass Load (ML)					Compression Ratio (CR)
				(a)							(b)

Compression efficiency η c Storage Power Density (kW/m 3) Sinusoidal Speed Constant Speed APIA AIA Optimal

	100%
	95%
	90%
	85%
	80%
	75%
	70%
	100	200	300	400	500	600	700	800	900	1 000
	Figure 2.11: Comparison on the efficiency as a function of storage power density for
		different compression trajectories (CR=10). Adapted from [141]

Table 2 .

 2 2: Summary table of the HTE technologies from the literature with main results and conditions

	Figure 2.20: RTE efficiency improvement of a LP-based CAES system by implementing
	various HTE concepts (uncertainty bars represent 95% confidence interval), adapted from
	[114]

2.3.10 Summary of empirical heat transfer correlations

	Nusselt number N u has been shown to be the most used non-dimensional number for esti-mating the heat transfer inside the LP during compression/expansion processes. Different correlations have been proposed and used based on their own applicable conditions. The generic form used to write N u as a function of Re and P r numbers is defined in Eq. 2.20 Reference Year A a b c Dittus & Boelter [28] 0.023 0.8 0.3 0.0 Zhang et al. [183] 2014 8.456 + 0.325Re 0.625 P r 1/3 Sieder& Tate [150] 0.023 0.8 0.3 0.14 Zhang et al. [180] 2014 N u = 2f r lc 1/2 L d Re lc P r Ramakrishnan et al. [132] 2016 2 + 0.589Ra 0.25 D 1+( 0.469 9/16 4/9 P r ) Annand [4] 0.70 0.7 0.7 0.0 Neu[103] 2020 6.67 Re D P r D L 0.36 and summarized in Table 2.3 together with other proposed correlations. Skelland [152] N u x = 1.41 3k+1 4k 1/3 D h ReP r x 1/3 Neu [103] 2020 6.17 Re D P r D L 0.48	Remarks Porous media; L/D > 1; obtained numerically Fully developed turbulent flow in circular tubes; valid only for L/D > 60 Fully developed flow Honeycomb geometry Used for hollow spheres [132] 0.7 < P r < 160; Re > 10000; L/D > 60 Internal combustion engines Laminar regime; 60 < L/D < 200 Can be applied to foams in a pipe [119] Turbulent regime; 60 < L/D < 200
	Adair et al.[1]	0.053	0.8	0.6	0.0	Reciprocating compressors
	Hamilton[47]	0.0245	0.8	0.6	0.0	Inner duct turbulent flow
	Wakao & Kaguei [154]			1.748 + 1.19Re 0.6 D P r 1/3	Porous media; Re based on mean pore diameter
	Liu & Zhou [86]	0.75	0.8	0.6	0.0	Reciprocating compressors
	Zukauskas [188]			9.48 + 0.29Re 0.6 D P r 1/3	Porous media Re based on mean pore diameter
	Hsieh et al. [51]	0.16	1.1	0.0	0.1	Reciprocating compressors
	Cengel [19]	0.664	0.5	1/3		Laminar flow in a pipe [166]
	Cengel [19]	0.023	0.8	0.3		Turbulent flow in a pipe [166]
	Kamiuto and Yee [63]	0.124 0.791 0.791		Open-cellular porous materials based on equivalent strut diameter
						Porous media
	Lu et al. [89]	0.76	0.4	0.37	0	1 ≤ Re D ≤ 40; D = (1 -e -1-0.04 D f )
						where D f is the filament diameter
						Porous media
	Lu et al. [89]	0.76	0.4	0.37	0	1 ≤ Re D ≤ 40; D = (1 -e -1-0.04 D f )
						where D f is the filament diameter
	Lu et al. [89]	0.52	0.5	0.37	0	40 ≤ Re D ≤ 10 3 ; Porous media
	Lu et al. [89]	0.25	0.4	0.37	0	10 3 ≤ Re D ≤ 2 × 10 5 ; Porous media
	Zhang et al. [178]		N u = 9.700 + 0.0879Re 0.792 P r 1/3	Porous media; L/D > 1; Interrupted-plate

Table 2 .

 2

	3: Summary of some empirical heat transfer correlations used for solid pistons and liquid pistons

Table 3 .

 3 .1 we can notice as T r reduces the density ratio between the phases increases. LATTICE BOLTZMANN METHOD 1: Vapor and liquid density of water at different reduction temperature T r

	Reduction Temperature(T r ) Vapor density Liquid density Density ratio
	0.6	0.001	8.7	8700.00
	0.65	0.0024	8.43	3512.50
	0.7	0.0126	8.07	640.48
	0.75	0.0494	7.68	155.47
	0.8	0.132	7.22	54.70
	0.85	0.283	6.63	23.43
	0.9	0.543	5.922	10.91
	0.95	1.0074	4.952	4.92

  26: LATTICE BOLTZMANN METHOD Start Compute density Using Eqs. 8.9 and 3.30, compute f i (x, t) and g i (x, t), then compute φ * (x, t + ∆t) and u * (x, t + ∆t) Calculate density Density ρ is calculated using equation 3.39 Compute pressure Using Eqs 3.46 and 3.48, compute p(x, t + ∆t). The iteration is repeated until |p n+1 -p n |/ρ <

	Compute velocity
	u

* (x, t + ∆t) Advance Advance one time step and return to Step 1 Stop

Table 3 .

 3

	Method	Advantages	Challenges
	Shan-Chen	Easy to implement Highly parallelizable	Thick interface Spurious currents esp. for MCMP Low Reynolds number
	Phase-Field (free-energy)	Lee-Lin Inamuro	High density ratio (1000) Very small spurious current Adapted for high Reynolds Relatively high density ratios (120)	Thick interface Very small Reynolds Heavy computational load Poisson equation for pressure correction
	VOF		Innovative method Sharp interface High density ratios (1000)	Mass conservation problems Small literature available

6: Tested methods comparison, advantages and remaining challenges

Table 4 .

 4 2: Evolution of air's viscosity using Sutherland's model According to liquid piston parameters and applying the equation 4.9, the evolution of the kinematic viscosity µ air is listed in the table 4.2. While c p,air is considered as constant.

	.1):
	C	120	-	Sutherland Coefficient
	T 0	291, 15	K	Reference Temperature
	µ 0 18, 27 • 10 -6 Pa.s Reference Viscosity
	Table 4.1: Sutherland equation's parameters
				104

Table 4

 4

	𝐷
	𝑨𝒊𝒓

.

[START_REF] Arjomand | Heat transfer analysis of liquid piston compressor for hydrogen applications[END_REF]

shows the Reynolds number for each phase and for each configuration. The Reynolds number of the water and air (Eqs. 4.12 and 4.13 ), calculated based on the liquid piston velocity U p ist, are estimated to be 1706 and 104, respectively.

Table 4 .

 4

5: Calculation of the Reynolds Number based on compression parameters

Table 5 .

 5 3: The parameters of the two studied walls temperature cases Figure 5.15 shows the evolution of the temperature of the two studied cases. The temperature gap ∆T = 7K at the beginning of the compression process. The gap increases as PISTON

	Configuration	Fast	Slow
	Diameter D (m)	0.0508
	Length L 0 (m)		0.353
	Walls temperature T walls (K)		300
	Expansion time t f (s)	3	8
	Initial Pressure P 0 (P a)

depending on expansion time -Patil model

	0.7									Config2 10s (num) Config2 20s (exp)
	0.6
	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
						t / t f
	0.84 0.86 0.8 0.9 Rendement 0.88 0.9 0.92 0.94 0.96 1 Efficiency 1.1	Yan 8s Numerical 10s Numerical 20s Expansion Efficiency Evolution Yan 3s Evolution rendement
	2	4	6	8	10	12	14	16	18	20	22
	0.7					t f
	t / t f model 0.5 Figure 5.17: Comparison of Expansion Efficiencies for the studied configurations using Patil 0 0.1 0.2 0.3 0.4 0.6 0.7 0.8 0.9 1 0.6

de detente pour differents temps de detente -Formule Yan

											BoYan 3s
											BoYan 8s
											Config2 10s (num)
											Config2 20s (exp)
				Expansion Efficiency
		0.85
		0.8
	Efficiency	0.75
		0.7
		0.65
		2	4	6	8	10	12	14	16	18	20	22
							t f

depending on expansion time -Yan model

	Evolution
	Yan 3s
	Yan 8s
	Numerical 10s
	Numerical 20s

  .1).

	Material	Borosilicate glass Copper Plexiglas (Poly(methyl methacrylate) (PMMA))
	c p	746	385	1200
	ρ (kg/m 3 )	2230	8400	1050
	λ W/m/K	1.25	380	0.1
	κ(×10 -6 m 2 /s)	1.2	117	0.107
	e

5

  Pa) P 0 /P f U pist (m/s)

	Fast (3s) 0.0508 0.353 Slow (8s)	0.3037	∼ 3 ∼ 8	300	12.5	∼ 6	Variable Variable

Table 8 .

 8 4: Parameters of expansion process of experimental Yan et al. study Yan et al. imposed a variable piston speed at the outlet during expansion process. The implemented profiles are presented in Figure

(a) (b)

(a) Density field at t * = 1.5 (b) Density field at t * = 5.5

(a) Velocity profile in the gas phase (b) Velocity profile in the liquid phase Figure 3.35: Velocity profile of the two phases in the "outlet" test this explain a part of the rapid pressure growth, it still does not explain the pressure field discontinuities. All the pressure correction algorithms tested and implemented have not corrected this problem. One other issue that can explain the numerical instabilities faced, is the low relaxation time due to the density ratios, it can go down to τ L = 0.505. Using a regularised model reduced to velocity field discontinuities at the beginning of the compression, but could not solve the other issues.

(a) Velocity field (U Z ) of the compression1 test (b) Pressure field showing discontinuities

classical approach consisting of Navier-Stokes equation solving coupled to VOF method for interface tracking is chosen for the rest of the study.

(a) 2D model U z velocity field at t = 6.7s (b) 3D model U z velocity field at t = 6.7s

Acknowledgments

Assumptions

The following assumptions are made in the numerical model:

• The compressed gas (air) is considered a perfect gas;

• The liquid phase (water) is considered as incompressible flow; • No mass transfer occurs between water and air.

Geometry and numerical parameters 4.2.2.1 Experimental setup

To enable experimental comparison and validation of numerical results, the physical parameters of the study are the same as the ones used by Neu et al. [START_REF] Neu | Experimental study of convective heat transfer during liquid piston compressions applied to near isothermal underwater compressed-air energy storage[END_REF].

In the experimental setup (Figure 4.3), a solid piston is used instead of a liquid one to better handle the Particle Image Velocimetry (PIV) measurements which are used to analyze the velocity fields.

While the pressure is measured by a pressure sensor and the piston position (which serves for the calculation of the volume (Eq. 4.10)) is measured by a position sensor of the solid piston, the air's average temperature is calculated using ideal gas law (Eq. 4.11).

V =

(L 0 -L pist πD 2 ) 4 (4.10) T = T air,ave = P T 0 P 0 V 0 = P V mr (4.11)

Numerical parameters

For the 3D geometry, the liquid piston is modeled as a cylinder with L 0 length and D diameter. While the 2D geometry is considered as a plane with L 0 length and D ∈ [-R, R] width. The water is injected at the piston bottom at a velocity corresponding to the experimental one, its average through time is U pist = 0.033 m/s (Figure 4.4). L pist is the position of the water/air interface at time t. The piston walls are considered as isothermal at T walls and the domain is at an initial temperature T 0 . The initial domain is only filled by air L pist = 0m and water is injected from the bottom at a constant velocity U pist . PISTON

Effects of the wall temperature on the compression efficiency

The compression efficiencies calculated for the three different numerical cases presented in previous section (4.6.2) are presented in figure 4.38. They show a linear evolution of η c as a function of the walls temperature. The higher the temperature of the walls, the lower the efficiency. This result is consistent with the fact that a high wall temperature indicates a loss of energy lost during the compression process and therefore a lower efficiency. Thus, the compression efficiency evolves in a linear way for the three models and follows the following equations :

Validation and comparison to experimental results

The obtained numerical results are compared to experimental ones obtained in a parallel study under the REMORA project. The only data available for comparison are the air's pressure and air's average temperature evolution through compression time. Data on the velocity fields are still lacking in the literature. The comparison criteria will therefore be the pressure and air's average temperature. Figure 5.4 shows the pressure evolution during the expansion process, it evolves from P 0 = 4, 14 × 10 5 P a to P f = 10 5 P a.

As in the compression process, the evolution of the pressure in the numerical and experimental results is similar, with a small gap (max 1500 Pa) (Figure 5.4). This is mainly due to the direct link between the pressure evolution and the air's volume one. The pressure decreases rapidly in the first 4s, where it is divided by 2. Then it decreases slowly in the last 16s to achieve the desired pressure P f = 10 5 P a.

The air's temperature evolution is a key parameter that define the efficiency of the process. Its evolution shows an interesting form which can be visible in the figure 5.5. The temperature decreases at the beginning of the compression. As in the compression, the temperature remains constant during the transition. The highest gap between numerical and experimental results is noted (∆T = 3K at t = 3.2s). The temperature rises then until the end of the expansion at t = 20.6s where T air,ave = 293K with less than 1K gap between the two curves.

While no precise experimental uncertainty analysis is available for the expansion setup, we used the same one as in the compression σ = ±1.6K, but with a widen range which gave PISTON the expansion progresses, it reaches 8K at t = 3s where the lowest temperature is observed. The gap is then constant at ∆T = 8K until the end of the process. It shows, as in the compression, that the walls temperature has a huge impact on the final temperature. We can consider that a ∆T difference on T walls will result in a ∆T gap in T f .

If

Effects of piston Speed

In order to investigate the effect of the piston speed on the expansion process two cases have been studied. The configuration is based on Yan et al. study [START_REF] Yan | Experimental study of heat transfer enhancement in a liquid piston compressor/expander using porous media inserts[END_REF], which has similar parameters than those studied previously (Table 5. Table 5.5: The parameters of the studied cycle Numerical results were compared to experimental ones. The piston speed is almost constant at 0.033m/s, except at the beginning and the end where a ramp of t = 0.1s is imposed. The final numerical position is approximately identical to the experiment with the negligible difference of 1mm. The evolution of the pressure is approximately the same. However, there is a notable difference at the end of the compression ∆P c = ∆P cool = 1 × 10 4 P a until the end of PISTON the cooling and the gap slightly increases at the end of the expansion ∆P cool = 1.5 × 10 4 P a. The reasons behind this gap are still not clearly identified, but could be due to leakage in the experimental setup especially at high pressure.

The pattern of the air's temperature during the compression is very similar until t = 5s, from there a difference starts to appear (Figure 5.22). This gap is amplified and at the end of compression is evaluated at ∆T = 6.5K. It was shown in the parametric compression study that a difference of ∆ on T walls would result in a ∆ gap of T f . Thus in comparison with the main case studied in the chapter 4, a T f = 334K is expected, which is the case for the numerical result but not for the experimental one. As for the pressure the main source of this gap could be the amount of substance's leakage.

The temperature then drops during the isochoric cooling, from T c to T cool . Numerically T cool = 301K which is 1K higher than T 0 . But experimentally the temperature decreases to lower values than the ambient temperature, which seems not to be possible in term of thermodynamics as shown in chapter 2 section 2.2. As in Neu et Subrenat [START_REF] Neu | Experimental investigation of internal air flow during slow piston compression into isothermal compressed air energy storage[END_REF], an invert method is used to calculate the temperature. Thus, a loss in the amount of substance would result in a temperature drop. To correct this loss, the temperature is readjusted to the T walls = T cool = 300K (Figure 5.22). General conclusions and perspectives

Overview

Compressed air energy storage (CAES) systems have demonstrated the potential for the energy storage of power plants. One of the key factors to improve the efficiency of CAES is the efficient thermal management to achieve near isothermal air compression/expansion processes. In this context, REMORA project is developed and proposes a CAES using a Liquid Piston. It has been found that LP is a promising concept for isothermal CAES. However, the complex and transient fluid flow and heat transfer behaviors inside the LP are yet to be characterized and mastered. The present study aims to model numerically the flow and heat transfer inside a LP during the compression/expansion processes.

To enhance the heat transfer and increase the efficiency of the compression/expansion processes many approaches have been tested including liquid spray, wire mesh, porous media, optimal trajectory, hollow spheres and optimal geometry of the piston column.

Chapter 2 presents a review on the Liquid Piston (LP) technology for CAES. Various aspects are discussed including the state-of-the-art on LP projects all over the world and the trend of development. It highlighted the complexity of the coupled fluid flow and heat transfer during the compression/expansion operations, and different actions proposed and implemented to enhance the heat transfer inside the piston column. While most of the studies try to characterize the heat transfer using empirical (pre-existent and developed) correlations, the flow evolution still lacks proper understanding. As it is difficult to have a complete visualizations of the flow experimentally, numerical modeling is chosen to study the flow and heat transfer inside a LP.

Chapter 3 presents the Lattice Boltzmann Method (LBM) which was initially chosen for our study. The main challenge in modeling the flow inside a LP using the LBM is the handling of liquid/gas interface for a multiphase flow with high density ratio. While, several methods for multiphase modeling already exist, they are limited in range of density ratio (1-10), while this ratio is 1000 for a water/air flow. Different approaches were implemented The LBM (Lattice Boltzmann Method) is originally based on the lattice gas automata [START_REF] Higuera | Boltzmann approach to lattice gas simulations[END_REF] and only being able to capture the hydrodynamics effects [START_REF] Shi | Theory of the lattice Boltzmann method: Lattice Boltzmann models for nonideal gases[END_REF]. The LBM closes the gap between the micro-scale (Molecular Dynamics) and the macro-scale (Navier-Stokes), thus, is is based on microscopic models and mesoscopic kinetic equations. The variables of concern in CFD are the macroscopic, observable quantities. In the incompressible Navier-Stokes equation, the degrees of freedom are defined by the space and time dependent velocity e = e(x, t) and the pressure p = p(x, t). In classical CFD, these variables are often discretized and directly calculated on the each cell of the grid. The LBM simulated the dynamics of each particle through the particle distribution function, a quantity derived from the kinetic theory of gases. It gives a probability density for the presence of particles in phase space. It is actually not always clear how a given configuration of the macroscopic variables can be used to initialize the distribution function. Reason why LBM is ms becoming more and more popular, is the fact that LBM is solved locally. It has high degree of parallelization, hence it is ideal for parallel machines (computational clusters) and can be solved in graphics Process Units (GPU).

List of publications

From Kinetic Theory to LBM

The fluid is made up of molecules moving at very high speeds of around 300m/s and interact with elastic collision [START_REF] Michael | Lattice Boltzmann Modeling: An Introduction for Geoscientists and Engineers[END_REF]. By knowing the position vector and the momentum of the individual particles at a given time, it is possible to predict the state of the system in future by utilizing the distribution function f N (x N , p N , t) (where N is the number of particles) with the help of classical mechanics [START_REF] Michael | Lattice Boltzmann Modeling: An Introduction for Geoscientists and Engineers[END_REF]. This methodology is directly used to simulate molecular dynamics which is discussed in the next section

Appendix B Numerical schemes

First realized tests showed a divergence of the simulations. Thus, the choice of numerical parameters, especially discretiezation schemes is a key element. After several tries, the chosen schemes are presented in tables 8.1, 8.2, 8.3. More details about each discretization paramater can be found in [START_REF] Christopher | {OpenFOAM} {Programmers} {Guide}[END_REF].

Variable Linear Solver

Comparison of turbulence models

The LES turbulence model was chosen for its capacity to model both the laminar flow regime, the transition and the turbulent regime. Two LES SGS model were chosen : WALE and Smagorinksy-Lilly. The two SGS models are presented in section 4.2. In order to choose the best one between these two SGS, a comparison study is conducted. In addition to air's average temperature evolution, local data comparison has also been used.

Comparison between experimental results and the two LES SGS turbulence models, show that no significant difference can be seen at air's average temperature evolution (Figure 8.8). One ∆T = 1K can be seen at the end of the compression process, which is an insignificant gap. where q is the heat flux, h is the heat transfer coefficient and T ext is the external media temperature. These three types of boundary conditions are widely employed as approximations at the fluid-solid interfaces in an uncoupled heat transfer analysis. A uniform temperature (isothermal walls) or heat flux, for example, could be assumed. Convective boundary conditions are also frequently utilized, in which the combined effects of the fluid flow are summed up in the heat transfer coefficient h, which must be calculated empirically. h is not a material feature, it is determined by a variety of parameters including shape, surface temperature distribution, flow conditions, and the fluid's thermophysical properties.

Appendix C Boundary conditions in heat transfer modeling