Isabelle Cheylan

Song Zhao

Pierre Boivin

Pierre Sagaut

Compressible Pressure-Based Lattice-Boltzmann applied to humid air with phase change

A new compressible pressure-based Lattice Boltzmann Method is proposed to simulate humid air flows with phase change. The variable density and compressible effects are fully resolved, effectively lifting the Boussinesq approximation commonly used, e.g. for meteorological flows. Previous studies indicate that the Boussinesq assumption can lead to errors up to 25%, but the model remains common, for compressible models often suffer from a lack of stability. In order to overcome this issue, a new pressure-based solver is proposed, exhibiting excellent stability properties. Mass and momentum conservation equations are solved by a hybrid recursive regularized Lattice-Boltzmann approach, whereas the enthalpy and species conservation equations are solved using a finite volume method. The solver is based on a pressure-based method coupled with a predictor-corrector algorithm, and incorporates a humid equation of state, as well as a specific boundary condition treatment for phase change. In particular, boundary conditions that handle mass leakage are also proposed and validated. Three test cases are investigated in order to validate this new approach: the Rayleigh-Bénard instability applied to humid air, the atmospheric rising of a condensing moist bubble, and finally the evaporation of a thin liquid film in a vertical channel. Results indicate that the proposed pressure-based Lattice-Boltzmann model is stable and accurate on all cases.

Introduction

Many application fields require to account for humid air physics when modeling and simulating flows. Examples of important industrial application fields are the simulation of high pressure water vapor flows in steam turbines and the simulation of cooling systems. This is also illustrated by atmospheric flows, in which atmospheric moisture plays a crucial role in the formation of clouds, storms, weather fronts, precipitation and other phenomena [START_REF] Haltiner | Dynamical and physical meteorology[END_REF]. It governs most aspects of our weather and climate in the form of vapor, liquid water, and ice. All these mechanisms can be solved numerically by invoking the Boussinesq approximation [START_REF] Zerroukat | A moist Boussinesq shallow water equations set for testing atmospheric models[END_REF], which consists in neglecting density variations except in the gravity term in the Navier-Stokes equations. This incompressibility condition simplifies the equations, and is widely adopted in weather forecast tools, but comes with a restriction on the range of flow physics. In fact, it is restricted to small variations of temperature [START_REF] Gray | The validity of the boussinesq approximation for liquids and gases[END_REF], which can be very limiting in the scope of meteorological applications, or phase change mechanisms. Cherkasov et al. [START_REF] Cherkasov | Limitations of the Boussinesq Model on the Example of Laminary Natural Convection of Gas between Vertical Isothermal Walls[END_REF] showed that the Boussinesq assumption can lead to an error of 25% in a problem of laminar natural convection in a gas between two vertical isothermal plates heated up to different temperatures, if the temperature of the hot wall is no more than three times the temperature of the cold wall. Schatzmann [START_REF] Schatzmann | Effects of the Boussinesq Approximation on the Results of Strongly-Buoyant Plume Calculations[END_REF] also studied a variety of industrial stack plumes with different initial density ratios, and made comparisons of predictions with and without the Boussinesq approximation. It was found that the difference in jet-property predictions between the two models increases with increasing initial density defect ratio. In the case of the largest ratio tested, the maximum height of rise predicted with the Boussinesq model was 20% larger than with the non-Boussinesq model. This shows the importance of a non-Boussinesq model to study cases with large differences of temperature, or more generally cases where the thermodynamic variables vary a lot. Other less restrictive assumptions, such as using pseudo-incompressible equations [START_REF] Wp. O'neill | A moist pseudo-incompressible model[END_REF], or anelastic equations [START_REF] Bannon | Theoretical foundations for models of moist convection[END_REF], can also be used, which all rely on simplifications of the fully compressible Navier-Stokes equations to tackle the compressibility issue. Previous studies have also demonstrated that it is possible to use a full compressible model: Feng and Tao [START_REF] Feng | A Compressible Thermal Lattice Boltzmann Model with Factorization Symmetry[END_REF] developed a full compressible thermal model using a multispeed lattice Boltzmann method, but the cost of this method is prohibitive for an industrial use as will be explained in section 3.2. Another solution to reduce the cost of a compressible solver is to couple a finite volume method to a Lattice-Boltzmann model as was shown by Luan et al [START_REF] Luan | Coupling of finite volume method and thermal lattice Boltzmann method and its application to natural convection[END_REF]. They applied their model to a case of natural convection and obtained stable results, but their study was restricted to single component flows. The main disadvantage of the compressible model is its lack of stability as shown by Guiho [10]. This stability issue of the compressible model will also be demonstrated with the third test case presented here, because a stable solution with the density-based solver could not be obtained. In order to overcome this problem, the objective of our study is to broaden the scope of application of the compressible model thanks to an improved solver stability. The novelty of our method is to propose a stable multi component compressible solver that can handle humid air for meteorological applications, and to achieve a better accuracy without using the common Boussinesq assumption. New boundary conditions will also be proposed in order to handle the mass leakage at boundaries. The rest of the introduction will focus on the features of our pressure-based compressible hybrid Lattice-Boltzmann method.

Most previous studies dealing with humid air and moist cases are based on the Navier-Stokes equations, e.g. [START_REF] Grabowski | Monotone Finite-Difference Approximations to the Advection-Condensation Problem[END_REF][START_REF] Smolarkiewicz | Simulation of all-scale atmospheric dynamics on unstructured meshes[END_REF][START_REF] Deardorff | Usefulness of Liquid-Water Potential Temperature in a Shallow-Cloud Model[END_REF]. Here, we chose to use the Lattice-Boltzmann method (LBM) which is a mesoscopic approach, based on the Boltzmann equation, that recovers the Navier-Stokes equations on the macroscopic level. The advantages of the LBM are its automatically generated cubic mesh which removes the need of a body-fitted mesh, as well as its very high parallel efficiency and very high numerical accuracy. The LBM has been gaining more and more interest in the industry, with practical applications in the automotive industry, the aerospace industry, but also in civil engineering, among other applications. Feng et al. [START_REF] Feng | Hybrid recursive regularized lattice Boltzmann simulation of humid air with application to meteorological flows[END_REF] have also developed the LBM for meteorological applications, and their study will serve as a guideline for this article. Nevertheless, they use the Boussinesq approximation which is not the case here as explained previously.

The novelty of the present work is that a new fully compressible pressure-based Lattice-Boltzmann solver with a hybrid recursive regularization is used and developed in order to incorporate specific features: the humid air equation of state, the phase change mechanisms and the evaporation/condensation boundary conditions. Because of the presence of a small fraction of liquid in the humid air, the perfect gas law cannot be used. Instead, a humid air equation of state is used, e.g. see [START_REF] Haltiner | Dynamical and physical meteorology[END_REF][START_REF] Sommeria | Three-Dimensional Simulation of Turbulent Processes in an Undisturbed Trade Wind Boundary Layer[END_REF]. This equation of state, which originates in a linearization of a full multispecies multiphase equation of state makes it possible to study cases with small fractions of water vapor and liquid. Phase change between the vapor and the liquid is also taken into account and evaporation/condensation processes can therefore be simulated. In order to do so, the instantaneous equilibrium phase change model given in [START_REF] Sommeria | Three-Dimensional Simulation of Turbulent Processes in an Undisturbed Trade Wind Boundary Layer[END_REF] is used to simulate the condensation of a bubble rising in a moist atmosphere, but also the evaporation of a thin liquid film at the walls of a heated channel [START_REF] Laaroussi | Effects of variable density for film evaporation on laminar mixed convection in a vertical channel[END_REF]. Special attention is thus given in order to incorporate the phase change mechanisms in the boundary conditions. In order to account for phase change mechanisms, an enthalpy equation is used in the present work instead of the entropy equation used in [START_REF] Farag | A pressure-based regularized lattice-Boltzmann method for the simulation of compressible flows[END_REF]. This hybrid method is based on the coupling of the mass and momentum equations solved by the LBM, with the enthalpy equation solved by a finite-difference method. Doing so, one does not need (when simulating fully compressible flows) to expand the equilibrium function up to an order too expensive for an industrial use. The regular D3Q19 lattice can thus be kept, and the general cost of the solver, which is fully compressible, stays reasonable.

Another aspect of the solution is the HRR collision model, which damps non-physical modes which can lead to instabilities at high Reynolds number flows. The HRR collision model stabilizes the solution and makes it possible to use the solver in a more general framework [START_REF] Jacob | A new hybrid recursive regularised Bhatnagar-Gross-Krook collision model for Lattice Boltzmann method-based large eddy simulation[END_REF][START_REF] Coreixas | Recursive regularization step for high-order lattice Boltzmann methods[END_REF]. The last novel aspect of this LBM algorithm is the development of a new segregated, pressure-based approach, coupled with a predictor-corrector method [START_REF] Farag | A pressure-based regularized lattice-Boltzmann method for the simulation of compressible flows[END_REF]. The predictor-corrector method is an algorithm in which the first step (predictor) consists in finding the weakly compressible solution. Then, the second step (corrector) is to correct the first step in order to recover a fully compressible solution. This type of algorithm is quite common in Navier-Stokes methods, and its efficiency depends on the robustness and capability of the first step. Here, only explicit schemes are used for both steps and the predictor step is based on the artificial compressibility method pioneered by Chorin [START_REF] Chorin | A numerical method for solving incompressible viscous flow problems[END_REF]. The pressure-based collision model is used for the first step (predictor) of the algorithm. Pressure-based methods present an advantage in some cases which is the following: by controlling explicitly the pressure oscillations in flows with a large density variation, pressure-based solvers have proven to be more robust and accurate than classical single relaxation time density-based LBM solvers. Indeed, computing the pressure directly from an evolution equation, and not from the equation of state, reduces spurious pressure oscillations in multifluid Navier-Stokes methods, see [START_REF] Karni | Hybrid Multifluid Algorithms[END_REF][START_REF] Abgrall | Computations of compressible fluids[END_REF] and in Lattice-Boltzmann Methods [START_REF] Filippova | A novel lattice BGK approach for low Mach number combustion[END_REF][START_REF] Filippova | A novel numerical scheme for reactive flows at low Mach numbers[END_REF][START_REF] Lee | Pressure evolution lattice-Boltzmann-equation method for two-phase flow with phase change[END_REF]. Nevertheless, the predictor step gives a pressure evolution equation which is not valid for high-Mach number flows, hence the need for a correction step to recover the full compressible solution.

This article is organized as follows: the first section presents the governing equations of the method, beginning with the macroscopic equations and the specific features needed for the humid air applications : the humid equation of state, the phase change mechanisms, and the special boundary conditions for evaporation/condensation. The numerical method is then presented, with the pressure-based LBM coupled with a predictor-corrector algorithm and the finite volume method to solve the species and the enthalpy conservation equations. In the second section, the test cases are presented in increasing order of complexity. First, the Rayleigh-Benard instability is studied, in a single component configuration (dry air), and then in a multiple components configuration (dry air with water vapor and water liquid). This first test case will prove that the pressure-based code gives consistent results in a simple configuration, and is compared to the density-based code. The second test case deals with a moist bubble rising in the atmosphere, which is condensing during its rise. The purpose of this second test case is to demonstrate that the full compressible effects are taken into account in a large domain, without using the Boussinesq assumption, which is most commonly used in the meteorological context as shown by previous studies [START_REF] Bannon | Theoretical Foundations for Models of Moist Convection[END_REF][START_REF] Maronga | The Parallelized Large-Eddy Simulation Model (PALM) version 4.0 for atmospheric and oceanic flows: Model formulation, recent developments, and future perspectives[END_REF]. Finally, the last test case is about a heated vertical channel with a thin liquid film on its walls. Dry air is entering the channel from the top, and produces the evaporation of the thin liquid film, thus creating vapor at the wall. This last case shows the superiority of our model compared to the density-based code, because the latter did not give usable results contrary to the former. All the test cases are performed with the pressure-based version and compared to the density-based version of the code when possible, and/or to the reference when available, in order to validate the accuracy and robustness of the pressure-based method.

Governing equations 2.1 Macroscopic equations

The humid air flow is assumed to follow the compressible Navier-Stokes equations, linking the local mass volume ρ, velocity vector u α and energy.

Under this formalism, mass and momentum conservation read

∂ρ ∂t + ∂ρu β ∂x β = 0 ( 1
)
∂ρu α ∂t + ∂ρu α u β + δ αβ p -Π αβ ∂x β -ρg = 0, ( 2
)
where p is the pressure, to be linked to the local state in Eq. ( 7), g is the gravity acceleration, and Π αβ is the stress tensor

Π αβ = µ ∂u α ∂x β + ∂u β ∂x α -δ αβ 2 3 ∂u γ ∂x γ . ( 3
)
Energy conservation is considered in its enthalpy h form,

∂h ∂t + u α ∂h ∂x α = 1 ρ (Π αβ ∂u α ∂x β - ∂q α ∂x α + Dp Dt ), (4)
where the pressure derivative Dp Dt is assumed to be negligible, and q α is the heat flux, assumed to depend on the temperature T as

q α = -λ ∂T ∂x α + ρ N k h k Y k V k,α , (5)
with λ the thermal conductivity. It is defined as:

λ = µ P r k Y k C p,k (6)
V k,α is the diffusion velocity and it is evaluated through a Fickian approximation.

Moist equation of state and phase change

To close the system of governing equations, a thermodynamic closure is required, linking the pressure with the mass volume and energy (or enthalpy / temperature). Considering a mixture of dry air, water vapor, and liquid water with respective mass fractions Y d , Y v and Y l , and assuming that mass fractions of vapor and liquid water remain small, one obtains the following equation of state [START_REF] Haltiner | Dynamical and physical meteorology[END_REF][START_REF] Sommeria | Three-Dimensional Simulation of Turbulent Processes in an Undisturbed Trade Wind Boundary Layer[END_REF],

ρrT = p(1 -0.61Y v )(1 + Y l ), ( 7
)
where r is the mass gas constant of dry air.

The enthalpy h is linked to the temperature T through

h = 3 k=1 h k Y k , (8)
with

h k = C p,k T + q k . (9)
C p,k is the constant pressure mass heat capacity of species k, and q k its reference enthalpy.

The mean heat capacity C p therefore reads

C p = k Y k C p,k (10)
For the rest of the demonstration, we introduce the following thermodynamic quantity:

θ = p ρc 2 s = rT c 2 s ( 11
)
as the ratio of thermodynamic perfect gas pressure to the classical athermal LBM pressure.

In considering three constituents, two additional conservation equations need to be appended to the system as

∂Y v ∂t + u α ∂Y v ∂x α = 1 ρ ∂ ∂x α (ρD y ∂Y v ∂x α ) -Q (12) ∂Y l ∂t + u α ∂Y l ∂x α = 1 ρ ∂ ∂x α (ρD y ∂Y l ∂x α ) + Q, ( 13
)
where D y is the diffusion coefficient, given by the Schmidt number of each specie:

D y = µ ρSc y ( 14
)
and Q is a source term, accounting for potential phase change. The dry air mass fraction is deduced as

Y d = 1 -Y v -Y l .
Phase change, i.e. condensation or evaporation, is formulated as in the NCAR general atmospheric model [START_REF] Washington | A January simulation experiment with the two-layer version of the NCAR global circulation model[END_REF]: condensation takes place when saturation is reached at a given grid point, and the specific humidity is bounded by the saturation value Y sat v :

Y sat v = p sat p -(1 -p sat ) , (15)
where = 0.622 is the ratio of molecular masses of water vapor and air. The saturation pressure is assumed to follow an Antoine law

p sat (T ) = 610.78e (17.269 T -273.16 T -35.86 ) . (16)
The phase change source term is then computed through a saturation adjustment [START_REF] Sommeria | Three-Dimensional Simulation of Turbulent Processes in an Undisturbed Trade Wind Boundary Layer[END_REF]:

Q = -Y l /δ t if Y v < Y sat v and Y l < ∆Y v ∆Y v /δ t otherwise (17)
where

∆Y v = C p rT 2 Y sat v L v (T ) 2 + C p rT 2 (Y v -Y sat v ), (18)
with the mass latent heat L v defined as

L v (T ) = h v (T ) -h l (T ). ( 19
)
3 Numerical method

We now present the hybrid numerical method used to solve the governing equations. Mass and momentum conservation equations ( 1) and ( 2) are solved via a pressure-based Lattice-Boltzmann solver, coupled with a finite volume solver for energy and species equations ( 4), [START_REF] Smolarkiewicz | Simulation of all-scale atmospheric dynamics on unstructured meshes[END_REF] and [START_REF] Deardorff | Usefulness of Liquid-Water Potential Temperature in a Shallow-Cloud Model[END_REF]. The novelty of the method relies on the fact that the full compressible effects are taken into account, contrary to previous studies on the matter which use the Boussinesq assumption as explained in the introduction. Another novel aspect of our approach is the mass conserving boundary conditions, which have been developed specifically to handle mass leakage at the boundaries when the simulation domain is big, which is typical of meteorological applications.

Pressure-based Lattice-Boltzmann solver coupled to a predictorcorrector approach

Lattice-Boltzmann methods rely on a discretization on space, time, and velocity of the Boltzmann equation [START_REF] Krüger | The Lattice Boltzmann Method. Principles and Practice[END_REF][START_REF] Guo | The Lattice Boltzmann Method and its applications in engineering[END_REF]. f i represents the variable solved by the LBM, and corresponds to the density distribution of particles at space and time (x, t) with velocity c i . It is a mesoscopic approach which recovers the macroscopic quantities through the moments of the distribution function f i :

     i f eq i = ρθ = p/c 2 s i c i,α f eq i = ρu α i c i,α c i,β f eq i = ρu α u β + ρθc 2 s δ αβ (20)
The Lattice-Boltzmann equation, discretized in space, time, and velocities, is written:

f i (x + c i δ t , t + δ t ) = f eq i (x, t) + (1 - 1 τ )f neq i (x, t) + 1 2 F E i (x, t) (21)
The term F E i (x, t) accounts for the gravity term F g i and for correction terms [START_REF] Farag | A pressure-based regularized lattice-Boltzmann method for the simulation of compressible flows[END_REF]. There are several ways to include an external force into the LBM. Here, we use the approach from [START_REF] Mohamad | A critical evaluation of force term in lattice boltzmann method, natural convection problem[END_REF] in which the gravity term is written as follows:

F g i = ρω i A α c iα c 2 s + (u α A β + u β A α )(c iα c iβ -c 2 s δ αβ ) 2c 4 s (22)
In this approach, when an external force is taken into account, such as gravity, the macroscopic velocity is also modified as:

u = 1 ρ i c i,α f eq i + A 2 (23)
The other correction terms included in F E i are detailed in Farag et al [START_REF] Farag | A pressure-based regularized lattice-Boltzmann method for the simulation of compressible flows[END_REF]. The equilibrium function is written:

f eq i = ω i [ρθ + c iα ρu α c 2 s + a (2),eq αβ H (2) iαβ 2c 4 s + a (3),eq αβγ H (3) iαβγ 6c 6 s ] ( 24
)
where H iαβ and H iαβγ are respectively the second and third order Hermite polynomials. The non equilibrium distribution function is also expanded in Hermite polynomials and reconstructed by a Recursive Regularized BGK LBM [START_REF] Jacob | A new hybrid recursive regularised Bhatnagar-Gross-Krook collision model for Lattice Boltzmann method-based large eddy simulation[END_REF][START_REF] Coreixas | Recursive regularization step for high-order lattice Boltzmann methods[END_REF] :

f neq i = ω i [ a (2),neq αβ H (2) iαβ 2c 4 s + a (3),neq αβγ H (3) iαβγ 6c 6 s ] (25)
In the present approach, the non-equilibrium coefficients a

(2),neq αβ and a

(3),neq αβγ are evaluated with the hybrid recursive collision model introduced in [START_REF] Jacob | A new hybrid recursive regularised Bhatnagar-Gross-Krook collision model for Lattice Boltzmann method-based large eddy simulation[END_REF]. This model is used to enhance numerical stability, as well as enforce a traceless lattice-Boltzmann tensor by removing its trace. It is observed to yield accurate and robust results for a broad range of applications and flow physics, e.g. [START_REF] Jacob | Wind comfort assessment by means of large eddy simulation with lattice Boltzmann method in full scale city area[END_REF][START_REF] Merlier | Lattice-Boltzmann large-eddy simulation of pollutant dispersion in complex urban environment with dens gas effects: Model evaluation and flow analysis[END_REF][START_REF] Feng | Hybrid recursive regularized lattice Boltzmann simulation of humid air with application to meteorological flows[END_REF][START_REF] Feng | Hybrid recursive regularized thermal lattice Boltzmann model for high subsonic compressible flows[END_REF][START_REF] Guo | An efficient lattice Boltzmann method for compressible aerodynamics on D3Q19 lattice[END_REF].

The pressure-based solver is very similar to the density-based solver, except that the macroscopic equations which are recovered are a slightly different due to the change in the zeroth order moment ρθ. In the pressure-based solver, the macroscopic equation on mass which is recovered is the following:

∂ρθ ∂t + ∂ρu β ∂x β = 0 (26)
This equation represents the pressure evolution equation in the low-Mach nearly incompressible limit, but it is not valid in fully compressible flows with strong dilatational effects. A few corrections are necessary in order to recover the full compressible limit (which correspond to the corrector step mentioned above). The corrections, done on the zeroth order moment and on the second order moment, are the following:

ρ(t + δt, x) = i f * i (t + δt, x) + ρ(t, x)[1 -θ(t, x)] (27)
a cor αβ = c 2 s δ αβ [ρ(t + δt, x)(1 -θ(t + δt, x)) -ρ(t, x)(1 -θ(t, x))] (28)
More details about the pressure-based LBM algorithm and on the correction terms can be found in Farag et al. [START_REF] Farag | A pressure-based regularized lattice-Boltzmann method for the simulation of compressible flows[END_REF].

Finite volume solver for species and enthalpy

A first-order explicit Euler scheme is used for the time integration :

φ n+1 = φ n + δ t [RHS C (φ n ) + RHS D (φ n ) + Φ] ( 29
)
where RHS C , RHS D represent the convection term and the diffusion term, respectively. The convective flux is constructed using the MUSCL scheme [START_REF] Kim | Methods for the Accurate Computations of Hypersonic Flows: I. AUSMPW+Scheme[END_REF], with the van Albada limiter [START_REF] Hirsch | Numerical computation of internal and external flows: The fundamentals of computational fluid dynamics[END_REF]. A classical second order central difference scheme is adopted for the diffusion term. The last term Φ is related to remaining source terms, which are treated explicitly in time. It is also possible to solve the scalar transport and the enthalpy equations by using the LBM, with the multiple distribution function model for example [START_REF] Feng | A three dimensional lattice model for thermal compressible flow on standard lattices[END_REF], or a multispeed approach [START_REF] Feng | A Compressible Thermal Lattice Boltzmann Model with Factorization Symmetry[END_REF]. Nevertheless, these methods are more cumbersome numerically than the hybrid method proposed here as noted by [START_REF] Feng | Regularized thermal lattice Boltzmann method for natural convection with large temperature differences[END_REF], which is why we decide to couple the LBM with a finite volume solver for the species and enthalpy equations.

Boundary conditions

In the LBM, several types of boundary conditions can be used, which can be classified mainly in two types : the conditions which deal with distribution functions directly (such as Bounce-Back and interpolated Bounce Back conditions [START_REF] Bouzidi | Momentum transfer of a lattice fluid with boundaries[END_REF][START_REF] Chun | Interpolated boundary condition for lattice Boltzmann simulations of flows in narrow gaps[END_REF]), and the conditions which deal with macroscopic quantities from which distribution functions are then deduced, or reconstructed, hence the term "reconstruction" method [START_REF] Verschaeve | A curved no-slip boundary condition for the lattice Boltzmann method[END_REF][START_REF] Ginzburg | Local second-order boundary methods for lattice Boltzmann models[END_REF] for this second approach. The reconstruction method can be summarized with the following steps:

• Depending on the desired boundary condition, we apply the Dirichlet or Neumann condition to macroscopic quantities : velocity, density, enthalpy, and species

• Reconstruction of f eq i from macroscopic values on boundaries (eq 24)

• Reconstruction of f neq i from macroscopic values on boundaries (eq 25)

• Computation of f i using f eq i and f neq i Details of the present implementation are available in [START_REF] Feng | Solid wall and open boundary conditions in hybrid recursive regularized lattice Boltzmann method for compressible flows[END_REF][START_REF] Wilhelm | An explicit power-law-based wall model for lattice Boltzmann method-Reynolds-averaged numerical simulations of the flow around airfoils[END_REF][START_REF] Wilhelm | A new explicit algebraic wall model for LES of turbulent flows under adverse pressure gradient[END_REF] and will not be repeated here for the sake of brevity.

For enthalpy and species which are solved via the finite volume approach, the macroscopic values of h, Y l , Y v are simply imposed (Dirichlet condition) or extrapolated (Neumann condition) on all boundaries via a first order interpolation.

This approach is most commonly used in industrial applications, because it is possible to use some already existing conditions (wall models for turbulent applications for example) which have been developed on macroscopic quantities, often for Navier-Stokes solvers. On the other hand, the mass conservation on solid impermeable boundaries is not necessarily respected in the reconstruction method. For many types of engineering problems for which the numerical domain is small and gravity does not play a big role, it is not an important concern. Nevertheless, the mass conservation issue with the reconstruction method can cause a problem for cases in which gravity plays an important role because the numerical domain is large (atmospheric test cases for example). In that case, the leakage of mass at the boundaries can lead to inaccuracy and stability issues. This type of problem was encountered here for the bubble ascending in a moist atmosphere case. It was therefore mandatory to respect the mass conservation at the boundaries. To do so, it is easier to use the first approach: handling distribution functions directly. For compressible problems, respecting the mass conservation at the boundary can be mathematically written:

dm dt = V dρ dt dV = - V ∇(ρu)dV = - S ρu n dS = 0 ( 30
)
with n the normal vector of the boundary S. Knowing that in the LBM, ρu n = i c i f i , this leads to a condition on the functions coming out and coming in the domain at the boundary, which resembles the bounce-back method.

outgoing i c i f i (x k , t n ) + incoming i c i f i (x k , t n ) = 0 (31)
If we use the bounce back condition for the incoming unknown direction:

outgoing i c i f col i (x k-1 , t n-1 ) + incoming i c i f col i (x k , t n ) = 0 (32)
If we suppose that the boundary is flat (which is the case in the bubble test case), and knowing that c i = -c opp(i) with opp(i) the opposite direction of i, the condition is written:

f col i (x k , t n ) = f col opp(i) (x k-1 , t n-1 ) (33)
with i the incoming direction.

Validation of the method

The proposed new pressure-based fully compressible Lattice Boltzmann Method for humid air will now be validated considering several test cases.

Rayleigh-Bénard instability

The Rayleigh-Benard instability is a natural convection phenomenon, in which a fluid is enclosed between two horizontal plates, the one below being heated and the one on top being cooled, so that the fluid develops a regular pattern of convection cells known as Benard cells. It is an important mechanism of mass and heat transfer in nature (meteorology, astrophysics) and in many industrial applications (heat exchangers). The gravity accelerates the flow downwards whereas the viscosity counteracts the motion so that a convection pattern appears. The Prandtl number P r and Rayleigh number Ra are defined as:

P r = ν D (34)
Ra = P r gH 3 ∆T ν 2 T 0 (35)
∆T is the temperature difference between the top and bottom walls, T 0 is the reference temperature set to 300.0 K, H is the domain height set to 1m, and ν is the kinematic viscosity given by the value of Ra in the simulations. For the test cases presented in this section, the Rayleigh number is equal to 10 4 and the Prandlt number is equal to 0.71. Another important parameter is the reference velocity used to scale the velocity on figure 3: Only dry air is simulated in this first test case, so that the pressure-based code can be validated and compared to results from Feng et al. [START_REF] Feng | Hybrid recursive regularized lattice Boltzmann simulation of humid air with application to meteorological flows[END_REF], and Ouertatani et al [START_REF] Ouertatani | Numerical simulation of two-dimensional Rayleigh-Bénard convection in an enclosure[END_REF] in which a nearly incompressible model was used along with a Boussinesq approximation. The square simulation domain is 1m × 1m, on a 50 × 50 grid. Figure 2 shows the temperature field once the steady state has been reached. It is in good agreement with reference [START_REF] Feng | Hybrid recursive regularized lattice Boltzmann simulation of humid air with application to meteorological flows[END_REF]. The results are in very satisfactory agreement with the reference data, meaning that the pressure-based LBM solver can be accurately used to investigate flows in the Boussinesq regime. The slight differences observed may be due to the fact that our pressure-based solver does not rely on the Boussinesq assumption, contrary to the references [START_REF] Feng | Hybrid recursive regularized lattice Boltzmann simulation of humid air with application to meteorological flows[END_REF][START_REF] Ouertatani | Numerical simulation of two-dimensional Rayleigh-Bénard convection in an enclosure[END_REF]. Indeed, we can assume that without the Boussinesq assumption, in the momentum and mass equations, the density variations (others than in the buoyancy force) are not neglected, hence leading to small differences between the reference and the results presented here.

u 0 = gβH∆T (

Multiple component case : moist air (dry air with water vapor and water liquid)

This test case presents a multiple component Rayleigh-Bénard instability, which is composed of dry air, water vapor, and water liquid. The water vapor and liquid fractions are set respectively to 0.002 and 0.001 on the top wall, and to 0 on the bottom wall. The bottom wall is heated at 309.5K and the top wall is cooled at 299.5K. Periodic boundary conditions are applied at the left and right walls. The domain is extended to observe several convection cells. The temperature field is displayed on figure 5. Because the reference article does not present results with the gravity applied to all three components, but only to water components, the comparison is made between the density-based version of the ProLB solver, and the present pressure-based version. Indeed, the densitybased version of ProLB has already been validated in many cases [START_REF] Guo | An efficient lattice Boltzmann method for compressible aerodynamics on D3Q19 lattice[END_REF][START_REF] Feng | Solid wall and open boundary conditions in hybrid recursive regularized lattice Boltzmann method for compressible flows[END_REF] and will serve as a reference guideline when no data is available in the literature. Figure 5 presents the temperature field: several convection cells can be observed in the extended domain. A very good agreement is obtained between the density-based and pressure-based methods, which means that the new pressure-based code can be accurately used to simulate moist convection cases, with several components. Here, both density-based and pressure-based codes give good results and demonstrate good accuracy and robustness.

Rising bubble

A rising moist bubble is studied in a two dimensional domain, with a height of 2.4 km and a width of 3.6 km, on a 720 x 480 grid with a 5m mesh size. The base state has a constant relative humidity of 20% with a pressure of 850 hPa and a surface temperature of 283 K. The bubble is initially located at the center of the domain, at z = 0.8 km. The kinematic viscosity ν is set to 1 m 2 /s and the Prandlt number to 1. The relative humidity RH is set to 100% inside the bubble, for r < 200m. A transition layer is assumed for 200m < r < 300m, defined as follows :

RH = 20% + 80%cos 2 ( π 2 r -200 100 ) (37)
The vapor mass fraction can then be deduced as:

Y v = .p sat RH p 0 -RH(1 -)p sat (38)
Concerning other initial conditions, we cannot use the reference initial conditions from [START_REF] Grabowski | Cloud-environment interface instability: rising thermal calculations in two spatial dimensions[END_REF][START_REF] Feng | Hybrid recursive regularized lattice Boltzmann simulation of humid air with application to meteorological flows[END_REF], which were based on a nearly incompressible approach. Indeed, in these references a constant density ρ 0 is used as a reference state. This is possible when invoking the Boussinesq approximation, because density variations are only taken into account in the gravity term via a modified and special formulation of the gravity. Here, we use a fully compressible solver, which means that setting a constant initial density is not correct nor representative of the hydrostatic atmospheric state. The correct way to initialize the simulation is to find the initial state (p, ρ, T ) which respects the perfect gas law, the hydro static equilibrium state, and another constraint on the temperature. Several assumptions on the atmospheric temperature field can be made: isothermal atmosphere, or linear variation of the temperature for example. Here, we chose to use the latest, which leads to the following set of equations:

     p = ρRT dp dz = -ρg dT dz = -K (39)
Solving the previous system of equations leads to the following expression of the initial state:

     p = p 0 (1 -Kz T 0 ) (g/RK) ρ = P RT T = T 0 -Kz (40)
The coefficient K is set to 0.0065 here, so that the temperature varies linearly from 283 K (ground surface) to 267.4 K (top surface).

Regarding boundary conditions, both boundary conditions explained in section 3.3 were tested on a simplified 2.4km high air column with identical initial conditions as the bubble test case. Because the fluid is initialized to the hydrostatic state, there should be no movement of fluid in the air column and the velocity should be equal to zero. When using the reconstruction method, we encountered mass conservation issues at the boundaries due to the presence of gravity in a very high domain. Nonphysical mass flux is coming in the domain because the top and bottom walls do not conserve exactly the mass with the reconstruction method. To solve this problem, we used the second boundary condition [START_REF] Merlier | Lattice-Boltzmann large-eddy simulation of pollutant dispersion in complex urban environment with dens gas effects: Model evaluation and flow analysis[END_REF] explained in section 3.3 which conserves the mass. The result with both boundary conditions is displayed on Fig. 8 The mass coming in the domain with the BB boundary condition is better controlled with incoming velocity values of the order of 1e -3 m/s compared to 1e -1 m/s with the reconstruction method. We will hence use the bounce back condition for the bubble test case.

The contour of water vapor is displayed on figure 9. The bubble rises and expands over time, condensing during its rise as predicted by the reference. Nevertheless, a quantitative comparison with the reference article is not possible because the initial states are different, as it is not possible to impose a constant reference density initially since the Boussinesq assumption is not used here. In order to illustrate this problem, figure 10 shows a comparison of the initial densiy and temperature, in the Boussinesq reference case and in our pressure-based compressible solver. Imposing a constant density field initially led to a non physical result in our case, because we do not use the Boussinesq assumption due to the compressible formulation of our solver, so it was mandatory to impose a varying density initially, which agrees with the hydrostatic solution. Because of this initial difference, the results quickly differ: at two minutes, the highest position of the bubble in our case is 1152m, 5% above the measured position 1100m in the reference. At three minutes, the order of error is 10%, the bubble being located at 1352m in our case, against 1194m in the reference. Both density-based and pressure-based codes give almost identical profiles. The results are in good agreement with [START_REF] Grabowski | Cloud-environment interface instability: rising thermal calculations in two spatial dimensions[END_REF], with the bubble expanding and condensing during its rise.

Heated vertical channel with evaporating liquid film

The last test case deals with a heated vertical channel with evaporating thin liquid films at the walls in which a steady, downward, laminar flow is streaming. This test case shows that the pressure-based LBM solver also includes condensation/evaporation boundary conditions, which are stable and accurate. In order to validate the results, a reference case from [START_REF] Laaroussi | Effects of variable density for film evaporation on laminar mixed convection in a vertical channel[END_REF] is selected. First, it is assumed that the liquid films are extremely thin so that the boundary conditions for heat and mass transfer from [START_REF] Lin | Analysis of combined buoyancy effects of thermal and mass diffusion on laminar forced convection heat transfer in a vertical tube[END_REF] can be used. At the inlet section, the fluid temperature, velocity, and vapor mass fraction are constant and uniform, and the pressure is interpolated. At the outlet section, the pressure is constant and equal to the atmospheric pressure, and the fluid velocity, temperature and vapor mass fraction are interpolated. At the walls, we assume that the gas and liquid film interface are in thermodynamic equilibrium, which leads to the following expression for the vapor mass fraction at the wall under the extremely thin liquid film assumption:

Y v,w = M v M a p sat (T w ) p -p sat (T w )(1 -Mv Ma ) (41)
Liquid evaporation produces a normal velocity at the wall, which is given by [START_REF] Lin | Analysis of combined buoyancy effects of thermal and mass diffusion on laminar forced convection heat transfer in a vertical tube[END_REF]: The inlet Reynolds number is 300 and the height and diameter of the channel are respectively 2m and 2cm. The mesh size is 0.5 mm, which gives a 40 x 4000 grid. The pressure p is assumed to be equal to the ambient pressure p a for heights of the order of 1m. This leads to a constant value of Y v,w , because the wall is also heated at a constant temperature T w . Here, we chose to simulate a case from [START_REF] Laaroussi | Effects of variable density for film evaporation on laminar mixed convection in a vertical channel[END_REF] with T w = T 0 = 327.5K and Y v,w = 0.1. The kinematic viscosity ν is equal to 1.74e -5 m 2 /s and the Schmidt number is equal to 0.529.

u t = 0 u w = -Dv 1-w i ( ∂ω ∂n ) i (42)
Results are presented on figures 16 and 17. The comparison is made between our results and those given in [START_REF] Laaroussi | Effects of variable density for film evaporation on laminar mixed convection in a vertical channel[END_REF] on downwards velocity along the channel, at 0.2 m and at 1m. The results are in good agreement with the reference data, both in magnitude and tendency. Nevertheless, the discrepancies between the reference and our results is bigger than in the first test case. Several reasons can be invoked to explain these discrepancies: first, the boundary conditions are not applied in the same manner. In [START_REF] Laaroussi | Effects of variable density for film evaporation on laminar mixed convection in a vertical channel[END_REF], source terms were introduced in the cells adjacent to the walls because direct implementation of special boundary conditions is not possible in the commercial software they use, which is a black box on this matter. Here, we directly impose the correct boundary conditions on the walls, without introducing the notion of adjacent cells. We prefer to use this second approach for the boundary conditions because it is identical to the original study [START_REF] Lin | Analysis of combined buoyancy effects of thermal and mass diffusion on laminar forced convection heat transfer in a vertical tube[END_REF], and also the article [START_REF] Azizi | Buoyancy effects on upward and downward laminar mixed convection heat and mass transfer in a vertical channel[END_REF] on which the authors of the reference [START_REF] Laaroussi | Effects of variable density for film evaporation on laminar mixed convection in a vertical channel[END_REF] rely. In fact, we assume that the liquid films are extremely thin so that boundary conditions for heat and mass transfer can be introduced as was suggested by Lin et al. [START_REF] Lin | Analysis of combined buoyancy effects of thermal and mass diffusion on laminar forced convection heat transfer in a vertical tube[END_REF]. Secondly, they use a special formulation of the compressible flow equations: the density does not fully depend on the pressure, but only on the ambient pressure. In Ansys Fluent which is the commercial software that they are using, this special formulation is called the "Incompressible Ideal Gas Law", for which the density depends only on the operating pressure p op and not on the local relative pressure field :

ρ = p op R/W T (43)
This not the case in our solver because it is fully compressible, which means that the density fully depends on the local relative pressure and temperature. This could explain the discrepancies observed because we are not using the same formulation for the density. Indeed, the authors are already observing differences between their Boussinesq results and their Incompressible Ideal Gas Law results, so it seems logical that differences are also observed for our fully compressible model.

Regarding the density-based version of the Lattice Boltzmann Method, stability issues were encountered in this test case, which led to numerical instability before steady state was reached. According to [START_REF] Farag | A pressure-based regularized lattice-Boltzmann method for the simulation of compressible flows[END_REF], pressure-based solvers were initially designed to handle low-Mach flows with large density variations, which is typically the case here (multiphase interfacial flow). This comes from the fact that pressure-based solvers explicitly control pressure oscillations because the pressure is directly computed using an evolution equation instead of deducing it from an equation of state. In this case, we can assume that the density-based solver is unstable because it does not handle well the pressure oscillations caused by the large density variations at the walls, due to the presence of a thin liquid film. In order to stabilize the density-based simulations in this test case, we tried to refine the mesh in order to better capture and control the pressure oscillations at the walls. The mesh was refined 10 times, leading to a mesh composed of 16 million cells, leading to a dramatic growth of the computational cost, therefore demonstrating that the pressure-based method is more stable and robust and efficient on coarse meshes than the density-based method on this type of test cases. The coefficient σ in the HRR model [START_REF] Jacob | A new hybrid recursive regularised Bhatnagar-Gross-Krook collision model for Lattice Boltzmann method-based large eddy simulation[END_REF] also had to be modified in order to stabilize the computation.

Conclusion

A new Lattice Boltzmann method for humid air able to fully account for compressibility effect was proposed and successfully assessed on several test cases. This method is an extension of the segregated pressure-based method proposed in [START_REF] Farag | A pressure-based regularized lattice-Boltzmann method for the simulation of compressible flows[END_REF] to the case of moist air with phase change. To this end, a modified equation of state is implemented, and the entropy equation used in the original scheme is replaced by an enthalpy equation, which is better suited to account for condensation and evaporation phenomena. New boundary conditions are also introduced in order to tackle the mass leakage issue at the boundaries, often encountered in large simulation domains. The present method is observed to be more accurate and more robust than the density-based Lattice Boltzmann method for humid air previously proposed in [START_REF] Feng | Hybrid recursive regularized lattice Boltzmann simulation of humid air with application to meteorological flows[END_REF]. Indeed, the third test case about the heated vertical channel was only studied and validated with our compressible pressure-based model, because the compressible density-based code did not produce stable results and the solution exploded. The pressure-based code is hence more general, since the densitybased method was restricted to the Boussinesq approximation and weak dilatational effects. In fact, the compressible model is more accurate than the Boussinesq assumption but suffers from stability issues. These were handled thanks to the pressure-based code which is more stable than the density-based code, as shown in this paper thanks to the different test cases.

36 ) 4 . 1 . 1

 36411 Single component case : dry airThis test case is studied in a square domain, the top wall being cooled at T C = 299.5 K and the bottom wall being heated at T H = 300.5 K. Adiabatic conditions are set for the temperature on the right and left walls. The temperature difference T H -T C is chosen to be small enough to allow the comparison with results based on the Boussinesq approximation.

Figure 1 :

 1 Figure 1: Schematic diagram of the single component Rayleigh-Benard case

Figure 2 :Figure 3 :Figure 4 :

 234 Figure 2: Temperature for the single component case : reference [14] (left) and pressure-based solver (right)

Figure 5 :

 5 Figure 5: Temperature for the three component Rayleigh-Benard case

Figure 6 :

 6 Figure 6: Velocity (cross) and density (square) on a horizontal line at the middle of the domain. Symbols = reference data (density-based method), solid line = pressure-based solver

Figure 7 :

 7 Figure 7: Velocity (cross) and density (square) on a vertical line at the middle of the domain. Symbols = reference data, solid line = pressure-based solver

Figure 8 :

 8 Figure 8: Vertical velocity in the air column for two types of boundary condition implementations

Figure 9 :

 9 Figure 9: Vapor contour at initial time (left) and at t = 150s (right)

Figure 10 :

 10 Figure 10: Initial condition for the density (left) and for the temperature (right)

Figure 11 :Figure 12 :Figure 13 :-

 111213 Figure 11: Profiles of density along a vertical line x=1800m at t=100s

Figure 14 :

 14 Figure 14: Profiles of vertical velocity along a vertical line x=1800m at t=100s

Figure 15 :

 15 Figure 15: Profiles of vertical velocity along a vertical line x=1800m at t=100s

Figure 16 :-

 16 Figure 16: Downwards velocity at 0.2m

Figure 17 :

 17 Figure 17: Downwards velocity at 0.2m

Acknowledgement

This work was performed using the ProLB software. Centre de Calcul Intensif d'Aix-Marseille is acknowledged for granting access to its high performance computing resources.