A Leroy

D Violeau

M Ferrand

L Fratter

A Joly

A new open boundary formulation for incompressible SPH

Keywords: SPH, incompressible, open boundary conditions

In this work a new formulation for inflow/outflow boundary conditions in an incompressible Smoothed Particles Hydrodynamics (ISPH) model is proposed. It relies on the technique of unified semi-analytical boundary conditions that was first proposed for wall boundary conditions in 2013, then extended to open boundaries in the framework of weakly-compressible SPH (WCSPH). An ISPH model relying on that formulation for solid boundaries was then proposed, which is the one considered here. It includes a buoyancy model for temperature effects and a k -turbulence closure. There are two main requirements for the imposition of open boundaries in ISPH: an algorithm to let particles enter and leave the domain, and the correct imposition of open boundary conditions on the fields. Regarding the algorithm for particles creation/destruction, it relies on the variation of mass of the particles located at the open boundaries. When the mass of such a particle reaches a threshold, a new particle is released. On the other hand, the imposition of open boundary conditions on the fields is done by prescribing the value of the boundary terms appearing in the semi-analytical formulation. The formulation was first validated in 2-D on a cut dam-break, a case of propagation of a solitary wave and a Creager weir. It was then extended to 3-D and tested on a 3-D circular pipe. A preliminary application case consisting of two connected pipes at different temperatures was then simulated. The results are promising since in all cases the fluid enters and leaves the domain as prescribed and generating none or very few reflected waves.

Introduction

The high computational cost of the Smoothed Particle Hydrodynamics (SPH) method restrains the range of application in fluid dynamics to rather short scale problems, so that an efficient formulation for open boundaries is required for its use on an industrial level. Such a formulation would make it possible to accurately represent flows in truncated domains and is a first step towards the coupling of SPH with other numerical methods. However, the prescription of boundary conditions in SPH is problematic, especially when it comes to open boundary conditions.

Regarding wall boundary conditions, several methods have been proposed until now, the most classical ones in the SPH community being i ) the imposition of repulsive forces at the boundary [START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF][START_REF] Monaghan | SPH particle boundary forces for arbitrary boundaries[END_REF], ii ) the use of ghost particles across the boundary [START_REF] Shao | Incompressible SPH method for simulating Newtonian and non-Newtonian flows with a free-surface[END_REF][START_REF] Issa | Numerical assessment of the smoothed particle hydrodynamics gridless method for incompressible flows and its extension to turbulent flows[END_REF][START_REF] Lee | Comparisons of weakly compressible and truly incompressible algorithms for the SPH mesh free particle method[END_REF] and iii ) the use of mirror particles across the boundary [START_REF] Oger | Specific pre/post treatments for 3-D SPH applications through massive HPC simulations[END_REF].

Another class of methods is also starting to gain importance, where the SPH interpolation is renormalised so as to restore its consistency close to boundaries [START_REF] Kulasegaram | A variational formulation based contact algorithm for rigid boundaries in two-dimensional SPH applications[END_REF]. The accuracy of the renormalisation process and the treatment of the boundary terms are the two key-issues of this technique. They led to several works [START_REF] Kulasegaram | A variational formulation based contact algorithm for rigid boundaries in two-dimensional SPH applications[END_REF][START_REF] Feldman | Dynamic refinement and boundary contact forces in SPH with applications in fluid flow problems[END_REF][START_REF] Ferrand | Improved time scheme integration approach for dealing with semi analytical boundary conditions in SPARTACUS2D[END_REF][START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF][START_REF] Mayrhofer | Investigation of wall bounded flows using SPH and the unified semi-analytical wall boundary conditions[END_REF][START_REF] Leroy | Unified semi-analytical wall boundary conditions applied to 2-D incompressible SPH[END_REF][START_REF] Violeau | Exact computation of SPH wall renormalising integrals in 3-D[END_REF] which made it possible to accurately represent flows close to solid boundaries and to prescribe arbitrary Neumann or Dirichlet boundary conditions. Ferrand et al. [START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF] obtained good results with a k -turbulence model where Neumann boundary conditions could be prescribed exactly on k and for the first time in SPH, the condition on the turbulent dissipation rate being non-homogeneous.

On the other hand, three types of methods have been proposed regarding open boundary conditions. First, methods where particles are let to enter or leave the domain through the use of buffer layers (see [START_REF] Vacondio | SPH modeling of shallow flow with open boundaries for practical flood simulation[END_REF] for example). The open boundary conditions are then prescribed by setting the physical quantities of the particles in the buffer zones. The second type of method relies on the semianalytical boundary conditions [START_REF] Kassiotis | Semi-analytical conditions for open boundaries in smoothed particle hydrodynamics[END_REF]. The open boundaries are discretised like walls: through a mesh with particles at the vertices. The mass of these particles is let to evolve so as to inject or remove mass from the computational domain and the boundary conditions are accurately prescribed through the boundary terms of the differential operators, which run over the boundary mesh elements (segments in 2-D, triangles in 3-D). This kind of method is more complex than the ones based on buffer layers but makes it possible to accurately prescribe open boundary conditions on the fields, as was shown in [START_REF] Kassiotis | Semi-analytical conditions for open boundaries in smoothed particle hydrodynamics[END_REF]. A third kind of method was recently proposed in [START_REF] Hirschler | Open boundary conditions for ISPH and their application to micro-flow[END_REF], based on the mirror particles technique.

In this work the mirror axis of each particle moves at the open boundary, and new particles are created behind that axis when it reaches a limit position.

Regardless of the method chosen to discretise open boundaries, with weakly compressible SPH (WCSPH) the prescription of the fields at inflow or outflow boundaries should be done based on the theory of Riemann invariants, since a hyperbolic system of equations is then solved [START_REF] Mahmood | Absorbing inlet/outlet boundary conditions for 2D SPH turbulent freesurface flows[END_REF][START_REF] Kassiotis | Semi-analytical conditions for open boundaries in smoothed particle hydrodynamics[END_REF]. Otherwise, prescribing both the pressure and the velocity at an open boundary may result in spurious waves. This is not the case when it comes to the incompressible SPH method (ISPH), where a pressure Poisson equation is solved. ISPH was developed as an alternative to WCSPH in order to improve the pressure estimation [START_REF] Lee | Comparisons of weakly compressible and truly incompressible algorithms for the SPH mesh free particle method[END_REF][START_REF] Xu | Accuracy and stability in incompressible SPH (ISPH) based on the projection method and a new approach[END_REF]. In this approach the correct prescription of boundary conditions is a key-point for the accuracy of the scheme. The first ISPH models were based on ghost or mirror particles for solid walls, which do not make it possible to accurately prescribe arbitrary boundary conditions on the pressure. With the projection method employed in these models, a non-homogeneous boundary condition should be applied on the pressure, which was made possible with the work by Leroy et al. [START_REF] Leroy | Unified semi-analytical wall boundary conditions applied to 2-D incompressible SPH[END_REF]. Indeed, they proposed an ISPH model based on the semi-analytical boundaries, with correct pressure wall boundary conditions. The accuracy of the model was shown to reach that of mesh-based methods, especially on confined flows. The formulation employed for the free-surface treatment in that work remains problematic, since it involves the detection of the free-surface particles and the prescription of their pressure to zero. A new free-surface formulation was proposed in [START_REF] Nair | An improved free surface modeling for incompressible SPH[END_REF], that does not involve any tracking of the interface anymore, but it would have to be adapted to the USAW framework. In the present work, we thus use the same formulation for the free-surface as in [START_REF] Leroy | Unified semi-analytical wall boundary conditions applied to 2-D incompressible SPH[END_REF].

ISPH is nevertheless a promising method for industrial applications, even regarding confined flows where mesh-based methods find it hard to simulate flows around moving bodies that come to touch each other for instance. In order to apply the method to real cases, a formulation for open boundaries is necessary.

Few works regarding open boundary conditions for ISPH have been proposed until now, some based on the buffer layers technique (see [START_REF] Yildiz | SPH with the multiple boundary tangent method[END_REF] for example) and the one cited earlier ( [START_REF] Hirschler | Open boundary conditions for ISPH and their application to micro-flow[END_REF]) based on mirror particles.

The present work is thus an extension of the work proposed in [START_REF] Leroy | Unified semi-analytical wall boundary conditions applied to 2-D incompressible SPH[END_REF] to the treatment of open boundaries in ISPH, based on the semi-analytical boundary conditions. The ISPH model considered here includes a k -turbulence closure (see also [START_REF] Violeau | Two attempts of Turbulence Modelling in Smoothed Particle Hydrodynamics[END_REF][START_REF] Violeau | One and two-equations turbulent closures for Smoothed Particle Hydrodynamics[END_REF]) and a buoyancy model. First, the model will be described with emphasis on the prescription of the pressure boundary condition at open boundaries. Then, the results obtained on several 2-D validation cases will be presented, before showing first results in 3-D.

Governing equations and modelling choices

The system of equations to be solved is composed of the incompressible Reynolds-Averaged Navier-Stokes (RANS) equations coupled to a heat equation and to the k -turbulence closure. The Boussinesq approximation is used to account for density variations so that the system reads [START_REF] Hervouet | Hydrodynamics of Free Surface Flows: Modelling with the Finite Element Method[END_REF]:

                                     ∇ • v = 0 dk dt = P + G -+ 1 ρ ∇ • (µ k ∇k) d dt = k (C 1 P + C 3 G -C 2 ) + 1 ρ ∇ • (µ ∇ ) dv dt = - 1 ρ ∇p + 1 ρ ∇ • (µ E ∇u) + g (1 -β(T -T 0 )) dr dt = v dT dt = K E ∇ 2 T (1)
In this system, r is the particle position, v is the Lagrangian velocity and u is the Eulerian velocity. Both velocities are equal in our SPH model except for boundary particles. When the k -turbulence closure is used, the Eulerian velocity of wall boundary particles is prescribed with a wall function (see [START_REF] Leroy | Unified semi-analytical wall boundary conditions applied to 2-D incompressible SPH[END_REF]).

On the other hand, at open boundaries the particles' Eulerian velocity is set as the ingoing/outgoing velocity, as we will see in Section 5.2. For all open boundaries, the Lagrangian velocity is equal to zero, whereas it is possible to represent moving solid walls [START_REF] Leroy | Unified semi-analytical wall boundary conditions applied to 2-D incompressible SPH[END_REF]. t is the time, ρ is the density (which is assumed constant in the framework of the Boussinesq approximation for weakly buoyant flows), k is the turbulent kinetic energy field, is its dissipation rate, µ is the dynamic molecular viscosity (which may be a function of the temperature), µ T is the dynamic eddy viscosity, µ E = µ + µ T is the effective (i.e. total) dynamic viscosity. We also define the two variables

µ k = µ + µ T σ k and µ = µ + µ T σ .
σ k , C 1 , C 2 and σ are model constants described in Table 1. The kinematic molecular viscosity is denoted by ν = µ ρ and the eddy viscosity by ν T = µ T ρ . ν T is modelled as a function of k and of as usual [START_REF] Launder | Mathematical models of turbulence[END_REF]:

ν T = C µ k 2 (2)
where C µ is a constant defined in Table 1. Note that the eddy viscosity is imposed equal to zero at the walls. In case of a laminar flow, ν T is set to zero and the k and equations are not solved.

g is the gravity field (of magnitude g = 9.81ms -2 ), p = p + 2 3 ρk with p the pressure, β is the coefficient of thermal expansion, T is the temperature field, T 0 is the mean temperature. K E = K +K T is the effective thermal diffusivity, with K the molecular thermal diffusivity and K T = ν T P r T , P r T being the turbulent Prandtl number, taken as 0.85 (see e.g. [START_REF] Rodi | Turbulence models and their applications in hydraulics[END_REF]). P is the production of turbulent kinetic energy and G is a buoyancy production/destruction term. C 3 is set to 1 if G ≤ 0 and 0 otherwise. P is calculated according to a mixed linear-quadratic model [START_REF] Guimet | A linearised turbulent production in the kmodel for engineering applications[END_REF]:

P = min C µ kS, ν T S 2 (3)
where S = √ 2S : S is the scalar mean rate-of-strain. G is a buoyancy production/destruction term modelled through [START_REF] Rodi | Turbulence models and their applications in hydraulics[END_REF]:

G = βK T ∇T • g (4)
3. SPH interpolation in the frame of unified semi-analytical wall boundary conditions

In this section we summarise the unified semi-analytical wall (USAW) boundary conditions used herein. In this work, fluid particles which do not belong to a boundary are called free particles a ∈ F. Solid and liquid boundaries are In this framework, the discrete SPH differential operators are different from the classical ones [START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF]. The antisymmetric form of the discrete SPH gradient of a field A at particle a reads:

v 1 s v 2
(∇A) a ≈ G γ,+ a {A b } = 1 γ a b∈P V b (A a + A b ) ∇w ab - 1 γ a s∈S (A a + A s ) ∇γ as (5)
where V b is the volume of particle b defined as V b = m b ρ , with m b the mass of the particles. The density being kept constant in ISPH, we will omit the particle subscript in its notation. The particles' mass is also constant, except for vertex particles located at open boundaries. On the other hand, w ab = w h (r a -r b ), with r the position and w h the SPH kernel: in this work the 5th order Wendland kernel [START_REF] Wendland | Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree[END_REF]. γ a is the wall renormalisation factor mentioned in the introduction, defined as in [START_REF] Kulasegaram | A variational formulation based contact algorithm for rigid boundaries in two-dimensional SPH applications[END_REF] and [START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF]:

γ a = Ω∩Ωa w(r a -r )d n r ( 6
)
where Ω is the fluid domain, Ω a is the compact support of the kernel at particle a and n is the space dimension. Note that γ a is equal to 1 far from boundaries, due to the normalisation property of the kernel. On the other hand, γ a is inferior to 1 when the kernel support intersects the wall, as in Figure 1. γ a is computed through the analytical formula proposed in [START_REF] Leroy | Unified semi-analytical wall boundary conditions applied to 2-D incompressible SPH[END_REF] in 2-D and through a dynamic governing equation in 3-D, as in [START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF]. ∇γ as is the contribution of segment s to the gradient of γ a , defined as:

∇γ as = ∂Ωs∩Ωa w(r a -r )n s d n-1 r ( 7
)
∂Ω s is the boundary area spanned by segment s and n s is the inward unit normal to the wall on s (see Figure 1). The following property holds [START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF]:

∇γ a = s∈S ∇γ as (8)
The terms ∇γ as are computed through an analytical formula in 2D [START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF] and in 3D [START_REF] Violeau | Exact computation of SPH wall renormalising integrals in 3-D[END_REF][START_REF] Mayrhofer | Unified semi-analytical wall boundary conditions in SPH: analytical extension to 3-d[END_REF].

It is also possible to define a discrete symmetric gradient:

(∇A) a ≈ G γ,- a {A b } = - 1 γ a b∈P V b A ab ∇w ab + 1 γ a s∈S A as ∇γ as (9)
where A ab = A a -A b and A as = A a -A s . In case the discrete gradient of a vector field is calculated, the formulae ( 5) and ( 9) remain unchanged except that A ab ∇w ab and A as ∇γ as are replaced by A ab ⊗∇w ab and A ab ⊗∇γ as respectively.

In ISPH, it is important to use skew-adjoint gradient and divergence operators since the projection method is based on this property at the continuous level [START_REF] Cummins | An SPH projection method[END_REF]. In this work we use an antisymmetric SPH gradient and a symmetric SPH divergence operator. Their skew-adjointness properties are reviewed in [START_REF] Ferrand | Improved time scheme integration approach for dealing with semi analytical boundary conditions in SPARTACUS2D[END_REF][START_REF] Mayrhofer | Investigation of wall bounded flows using SPH and the unified semi-analytical wall boundary conditions[END_REF].

The symmetric form of the SPH divergence operator reads:

(∇ • A) a ≈ D γ a {A b } = - 1 γ a b∈P V b A ab • ∇w ab + 1 γ a s∈S A as • ∇γ as (10)
Finally, the discrete Laplacian operator proposed by Ferrand et al. [START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF] reads:

[∇ • (B∇A)] a ≈ L γ a {B b , A b } = 2 γ a b∈P V b B ab A ab r 2 ab r ab • ∇w ab - 1 γ a s∈S [B s (∇A) s + B a (∇A) a ] • ∇γ as ( 11
)
where B is a diffusion coefficient for the field A, B ab is a mean value between a and b (here the arithmetic mean) r ab = r a -r b and r ab = |r ab |. In case A is a vector, the Laplacian will be denoted by L γ a {B b , A b }. In case B = 1, it will be denoted by L γ a {A b }. The boundary term of the Laplacian operator is simplified so that the SPH Laplacian operator used here reads:

[∇ • (B∇A)] a ≈ L γ a {B b , A b } = 1 γ a b∈P V b (B a + B b ) A ab r 2 ab r ab • ∇w ab - 2 γ a s∈S B s (∇A) s • ∇γ as (12)

Incompressible SPH with the unified semi-analytical boundary conditions

In the ISPH model considered here the Chorin projection method [START_REF] Chorin | On convergence of discrete approximations to the Navier-Stokes equations[END_REF] is used to solve the Navier-Stokes equations, as in [START_REF] Lee | Comparisons of weakly compressible and truly incompressible algorithms for the SPH mesh free particle method[END_REF][START_REF] Leroy | Unified semi-analytical wall boundary conditions applied to 2-D incompressible SPH[END_REF]. The space-time discretisation of system (1) thus reads, ∀a ∈ F, ∀b ∈ P ∪ S:

                                                 k n+1 a -k n a δt = P n a + G a -n a k n+1 a k n a + 1 ρ L γ a {µ k,b , k n b } n+1 a -n a δt = n a k n a C 1 P n a + C 3 G a -C 2 n+1 a + 1 ρ L γ a {µ ,b , n b } ṽn+1 a -v n a δt = L γ a {ν n+1 T,b , u n b } -[β(T n a -T 0 ) -1] g L γ a {p n+1 b } = ρ δt D γ a {ṽ n+1 b } v n+1 a -ṽn+1 a δt = - 1 ρ G γ,+ a {p n+1 b } T n+1 a -T n a δt = L γ a {K n+1 T,b , T n b } r n+1 a = r n a + δtv n+1 a ( 13
)
where ṽn+1 a is a predicted velocity field. A particle shift is also used for stabilisation reasons. It consists in an additional particle motion based on a Fickian diffusion law, adapted to the USAW boundary conditions. For a more detailed description of this particle shift see [START_REF] Xu | Accuracy and stability in incompressible SPH (ISPH) based on the projection method and a new approach[END_REF][START_REF] Lind | Incompressible smoothed particle hydrodynamics for free-surface flows: A generalised diffusion-based algorithm for stability and validations for impulsive flows and propagating waves[END_REF][START_REF] Leroy | Unified semi-analytical wall boundary conditions applied to 2-D incompressible SPH[END_REF].

The time-step size δt must satisfy several conditions to ensure the stability of the scheme (see [START_REF] Violeau | Optimal time-step for incompressible SPH[END_REF]). Here an empirical criterion on the time step size is used, as in [START_REF] Leroy | Unified semi-analytical wall boundary conditions applied to 2-D incompressible SPH[END_REF].

The wall boundary conditions on the pressure are obtained by projecting the third line of (13) onto the normal to the wall, and prescribing the wall boundary condition on the velocity: v n+1 • n = 0. For all the segments s belonging to wall boundaries, this yields:

∂p ∂n n+1 s = ρ δt ṽn+1 s • n s = ρg + µ∇ 2 u n s • n s (14)
As in [START_REF] Lee | Comparisons of weakly compressible and truly incompressible algorithms for the SPH mesh free particle method[END_REF][START_REF] Lind | Incompressible smoothed particle hydrodynamics for free-surface flows: A generalised diffusion-based algorithm for stability and validations for impulsive flows and propagating waves[END_REF][START_REF] Leroy | Unified semi-analytical wall boundary conditions applied to 2-D incompressible SPH[END_REF], the free-surface particles are detected through a criterion on the value of the position divergence and a zero pressure is prescribed at freesurface particles. This ISPH model with USAW boundary conditions for walls was widely verified in [START_REF] Leroy | Unified semi-analytical wall boundary conditions applied to 2-D incompressible SPH[END_REF][START_REF] Leroy | Buoyancy modelling with incompressible sph for laminar and turbulent flows[END_REF].

Open boundaries for ISPH with the semi-analytical boundary conditions

There with N sv the set of segments s directly connected to v, N sv its size, S s the surface of segment s (or length in 2-D). At inflow boundaries, the mass flux [START_REF] Kassiotis | Semi-analytical conditions for open boundaries in smoothed particle hydrodynamics[END_REF] is positive and the mass of each vertex v increases until it reaches a higher threshold set to +0.5m ref , with m ref the mass of a free particle. Then, a free particle is created at the same location while m ref is subtracted to m v , so that m v goes down to -0.5m ref .

∀v ∈ V i/o , ṁn v = 1 N sv s∈Nsv ρS s (u s -v s ) • n s (15) s v (a) (b) (c) (d)
Figure 2 illustrates that process. In this way the mass variation is smooth with respect to space and time. At outflow boundaries, the mass flux is negative and when a free particle crosses a segment to get out of the domain it is deleted and its mass is distributed onto the vertices directly linked to the segment, a weight β a,v being associated to each of these vertices.

The weights are computed as in [START_REF] Kassiotis | Semi-analytical conditions for open boundaries in smoothed particle hydrodynamics[END_REF]:

• in 2-D, for v 0 and v 1 connected to s:

β a,v0 = p 1 • r v0v1 |r v0v1 | 2 β a,v1 = p 0 • r v1v0 |r v0v1 | 2 = 1 -β av0 (16)
• in 3-D, for v 0 , v 1 and v 2 connected to s:

β a,v0 = 1 2 [p 2 × r v2v1 ] • n s 1 2 [r v0v1 × r v0v2 ] • n s β a,v1 = 1 2 [p 0 × r v0v2 ] • n s 1 2 [r v0v1 × r v0v2 ] • n s β a,v2 = 1 2 [p 1 × r v1v0 ] • n s 1 2 [r v0v1 × r v0v2 ] • n s (17)
An illustration of the notations and of the fraction of segment area β a,v attributed to a vertex is provided on Figure 3. p i is the projection of r avi on s:

(a) (b) The coefficient associated to the vertex v 0 is proportional to the red area so that the largest amount of mass goes to the closest vertex particle.

p i = r avi -(r avi • n s )n s (v i
being one of the vertices linked to s). In this way, the authors of [START_REF] Kassiotis | Semi-analytical conditions for open boundaries in smoothed particle hydrodynamics[END_REF] took care that the largest amount of mass is attributed to the closest vertex particle to the point where a fluid particle is destroyed. This aimed at making the mass displacement as small as possible.

The following mass evolution equation is thus solved at the end of each time-step ∀v ∈ V i/o :

m n+1 v = m n v + δt ṁn v + δm n v ( 18
)
with δm n v the mass variation due to particle creation/destruction and ṁn v the mass flux corresponding to the imposed velocity at the open boundary.

Imposition of the inflow/outflow boundary conditions

We consider two options for the prescription of inflow/outflow boundary conditions: either the velocity or the pressure is set. Generally, boundaries where the velocity is imposed are inflow boundaries whereas boundaries where the pressure is imposed are outflow boundaries, but this is not exhaustive. Indeed, the algorithm makes it possible for particles to leave the domain through

Boundaries with prescribed velocity

The Dirichlet condition on the Eulerian velocity u is imposed at the vertex particles in the correction step of the projection method (fifth line of ( 13)). The velocity gradient at the segments is then deduced through a linear interpolation involving the particles close to the boundary. The Dirichlet condition on T is imposed at the vertex particles and the normal temperature gradient is computed through a linear interpolation between the segments and the surrounding free particles. The homogeneous Neumann condition on the pressure is imposed through the pressure Laplacian in the pressure Poisson equation (fourth line of ( 13)). Moreover, the pressure of vertex particles belonging to an inlet is extrapolated from the surrounding fluid particles so that a homogeneous Neumann condition is imposed. This ensures that the pressure at the open boundary is compatible with the imposed velocity. Thus, the boundary conditions imposed during the time-scheme when the velocity is prescribed read:

             u n+1 v = u pv v T n+1 v = T pv v ∂p ∂n n+1 s = 0 (19)
The compatible conditions on the complementary entities read:

                   ∂u ∂n n+1 s = u as δr as • n s ∂T ∂n n+1 s = 0 p n+1 v = p pv v = 1 α v b∈F V b p n b w vb (20)
where the superscript pv denotes the value of a field at an open boundary with prescribed velocity. In addition, δr as is defined by

δr as = max(r as • n s , δr) (21)
and α v is the Shepard filter:

α v = b∈P V b w vb (22)
The values of the fields at the segments of the boundary are then deduced from a mean of the directly linked vertex particles. The Neumann conditions in equations ( 19) and ( 20) are imposed in the boundary terms of the SPH Laplacian of T , p and u (second term in the right-hand side of ( 12)).

Boundaries with prescribed pressure

The Dirichlet condition on the pressure is imposed at the vertex particles. It can be either a fixed pressure value (like the hydrostatic pressure) or a radiative condition such as the one proposed by Orlanski [START_REF] Orlanski | A simple boundary condition for unbounded hyperbolic flows[END_REF], so as to let waves leave the domain:

∂p ∂t + C ∂p ∂n ∂Ωo = 0 (23)
with C a celerity usually taken as √ gd, d being the initial elevation of the free-surface above the bed at the outlet which is given by the user. Then, the pressure condition reads:

         p pp s = p n s - Cδt α s b∈F V b p n b -p n s δr sb w sb p pp v = 1 N sv s∈Nsv p pp s ( 24
)
where the superscript pp denotes the value of a field at an open boundary with prescribed pressure. The compatible Neumann condition imposed on the pressure is obtained through a linear interpolation of the surrounding free particles pressure. On the other hand, the homogeneous Neumann condition on T and u is imposed at the segments in the boundary term of the Laplacian operator [START_REF] Leroy | Unified semi-analytical wall boundary conditions applied to 2-D incompressible SPH[END_REF]. A compatible Dirichlet condition is deduced through an interpolation on the surrounding free particles. Thus, the boundary conditions imposed during the time-scheme when the pressure is prescribed read:

                 ∂u ∂n n+1 s = 0 ∂T ∂n n+1 s = 0 p n+1 v = p pp v ( 25
)
and for the complementary entities the compatible conditions read:

                   u n+1 v = 1 α v b∈F V b u n+1 b w vb T n+1 v = 1 α v b∈F V b T n+1 b w vb ∂p ∂n n+1 s = p n+1 a -p pp s δr as (26)

Solving the pressure Poisson equation

In the framework of the USAW boundary conditions, the pressure Poisson equation (second line of ( 13)) reads:

     2 γ a b∈P V b p n+1 ab r 2 ab r ab • ∇w ab - 2 γ a s∈S ∇p n+1 s • ∇γ as      = ρ δt D γ,- a {ṽ n+1 b } (27)
where the unknowns are the set of pressures p n+1 (recall that p ab = p a -p b ).

Taking the boundary conditions described in sections 4 and 5 into account (equation ( 14) and 3rd lines of ( 19) and ( 26)) and keeping in the left-hand side only the terms involving the unknown pressures, this equation becomes: Besides, the product of the columns corresponding to these particles with the unknown pressure vector is known (either through a Dirichlet condition or an approximation of a Neumann condition) and passed to the right-hand side.

2 γ a       b∈P V b p n+1 ab
In the end the system to be solved does not involve the inlet/outlet particles anymore and they appear in the right-hand side:

2 γ a             b∈P\V i/o V b p n+1 ab r 2 ab r ab • ∇w ab + b∈V i/o V b p n+1 a r 2 ab r ab • ∇w ab - s∈Spp p n+1 a δr as |∇γ as |             = ρ δt D γ,- a {ṽ n+1 b } + 2ρ γ a s∈S\S i/o ṽn+1 s -v wall s δt •∇γ as - 2 γ a s∈Spp p pp s δr as |∇γ as | + 2 γ a b∈V i/o V b p i/o b r 2 ab r ab • ∇w ab (29)
where p i/o b denotes either p pv b or p pp b . This equation corresponds to a linear system:

Ap = B ( 30
)
where p is the unknown vector of all particles pressures p a , B is the vector of right-hand side values at all particles and A is a sparse matrix corresponding to the discrete Laplacian operator. The Laplacian matrix is non-symmetric because of the γ a factor, of the term involving outlet segments and because V b is not constant with the USAW boundary conditions: the volumes of the vertex particles are lower than that of the free particles. To solve system [START_REF] Chorin | On convergence of discrete approximations to the Navier-Stokes equations[END_REF], the Bi-CGSTAB linear solver [START_REF] Vorst | Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems[END_REF] is used.

2-D validation

Schematic dam-break with an outflow condition

This case consists of a schematic dam-break on a flat bottom, which was cut so as to test the outlet formulation, and check that the fluid leaves the domain without reflections. The outlet boundary is the left-wall (at all times)

in Figure 4, where the pressure is prescribed according to the Bernoulli equation:

p pp = -1 2 ρv pp2
. All the lengths are made dimensionless by the height of the fluid column at the initial time, d. The dimensionless width of the domain is equal to 1. The initial dimensionless interparticular space was taken equal to 6 × 10 -3 .

The viscosity of the fluid was set to a relatively large value, 10 -2 m 2 s -1 , because the ISPH model used here still suffers from particle leaks for free-surface flows with low viscosity (see [START_REF] Leroy | Unified semi-analytical wall boundary conditions applied to 2-D incompressible SPH[END_REF]).

The dimensionless time is defined by:

t + = t √ gd (31)
Figure 4 shows the velocity field shape at several dimensionless times. On the left boundary, some particles adhere to the wall due to the relatively high value of viscosity used for this case. The fluid correctly leaves the domain without visible reflections at the outlet. The free-surface shape of the same non-cut dam-break simulated with ISPH-USAW is provided and appears in black in the Figure . The agreement is good between the two simulations, which shows that the proposed formulation is reliable.

2-D solitary wave

A case of propagation of a solitary wave on a slope is presented here, where the wave breaks before leaving the domain. the case. The incoming free-surface elevation is prescribed as a solitary wave (solution to the Korteveg-De Vries equation) [START_REF] Dean | Water wave mechanics for engineers and scientists[END_REF]:

η(x, t) = Asech 2 [k(x -Ct -x 0 )] ( 32
)
where η is the free-surface elevation compared to a reference water level d, A

is the wave amplitude, taken as A = d 2 , k = 3A 4d 3 is the wave number and C = g(A + d) is the wave celerity. x 0 is the initial position of the wave, equal to x 0 = x inlet -4 k here. At the inlet (left boundary in Figure 5), the time-dependent water height d t is used to impose the following velocity profile:

             d t (t) = d + η(x inlet , t) u x (z, t) = C η(x inlet , t) d(t) u z (z, t) = z d t (t) ∂η ∂t (x inlet , t) (33)
with x inlet the horizontal coordinate of the inlet. At the outlet (right boundary in Figure 5), the pressure is imposed through the Orlanski radiative boundary condition [START_REF] Launder | Mathematical models of turbulence[END_REF] with a celerity of √ gd.

Figure 6 shows the propagation of the solitary wave on the slope with a dynamic molecular viscosity of 10 -6 m 2 s -1 (no turbulence model was used in the simulation). The dimensionless time t + is defined as in the previous section, through [START_REF] Lind | Incompressible smoothed particle hydrodynamics for free-surface flows: A generalised diffusion-based algorithm for stability and validations for impulsive flows and propagating waves[END_REF]. The colours correspond to the pressure field obtained with ISPH-USAW. This simulation was run with 30315 particles (δr = 0.01m). The wave correctly enters the domain and leaves it after breaking, apparently without reflections, although there was no attempt at quantifying the amount of reflection in the domain.

This shows that the proposed open boundary formulation makes it possible

Figure 6: Propagation of a solitary wave on a slope with ν = 10 -6 m 2 s -1 : pressure field obtained with ISPH-USAW ranging from 0P a to 5837P a (blue to red).

to generate waves and let them leave the computational domain quite smoothly, even when the free-surface shape is complex due to wave breaking.

Creager Weir

A Creager weir is a spillway with a geometry such that for a given upstream head (called the dimensioning head H d ), the pressure on the weir is equal to the atmospheric pressure. The weir has a width L, which will not be considered here, as we are in a 2D simulation. The geometry is described on Figure 7.

The upstream face, between points M 1 and M 2 , is defined by three arcs of circles having centres with coordinates x i and z i and radii R i . Their values and the arc boundaries are given in [START_REF] Violeau | Fluid Mechanics and the SPH method[END_REF]. The downstream face is given by:

z H d = - 1 2 x H d 1.85 (34)
The hydraulic head is defined with:

H = 1 2g U 2 + d + P atm ρg ( 35
)
where U is the velocity, d the water depth and P atm the atmospheric pressure. Dimensional analysis gives the flow rate per meter q as a function of g, H,

a function C q H H d , H
Hp and H p the weir height, see Figure 7:

q = 2gH 3 C q H H d , H H p (36)
The function C q is also called the discharged coefficient, it is determined by experiments [START_REF] Chow | Open-Channel Hydraulics[END_REF]:

C q = C p H H d 0.12 (37)
From experiments C p can be set as C p ≈ 0.485. Here, the value of the dimensional head is chosen as H d = 0.3m, and a weir height H p = 0.4m. We choose to simulate two steady cases of imposed upstream head: H/H d = 0.5

and H/H d = 1.33. With the relations ( 36) and [START_REF] Violeau | Fluid Mechanics and the SPH method[END_REF] and knowing that q = hU , the water depth d and the bulk velocity U at inlet can be deduced. The inlet condition is imposed at a distance 4H p before the weir. At the beginning of the simulation, the fluid is at rest and the pressure is set as hydrostatic:

p(z) = ρg(d -z) (z
is the vertical elevation above the bed). For this case the k -turbulence closure was used. At the inlet, a logarithmic horizontal velocity profile is imposed, with the corresponding k and profiles [START_REF] Launder | Mathematical models of turbulence[END_REF]: Here, u * is the shear velocity upstream of the weir, ν = 10 -6 m 2 /s the kinematic viscosity of water, ρ = 1000kg/m 3 the reference density and C r = 5.2 is a constant. The mean (depth-averaged) velocity U is given by:

                 u(z) = u * 1 κ ln zu * ν + C r k(z) = u * 2 C µ 1 - z d (z) = u * 3 κz 1 - z d (38)
U = u * 1 κ ln u * d ν -1 + C r (39)
Knowing U and d, the shear velocity is iteratively computed. The imposed water depth, mean velocity and shear velocity are summarised in Table 2. At the outlet, we want to let the flow go out. To be correct, a hydrostatic pressure should be imposed. However, as the water depth is not known in advance, the pressure is imposed at zero.

Figure 8 shows the velocity field obtained at steady-state for the two configurations with ISPH-USAW, using δr = 0.0125H p . The free-surface position and pressure along the weir obtained with ISPH-USAW are compared to experimental results [START_REF] Usbr | Design of small dams[END_REF] for each configuration. The results are shown in Figure 9. The free-surface position obtained with ISPH-USAW shows quite good agreement with the experiment data, and the pressure along the weir shows reasonable agreement. The results obtained for the configuration with the lower head upstream match the experiments better than with the higher head. With the higher head, a pressure peak is observed in front of the weir in the numerical results. This is probably due to an issue with the convergence of the pressure The dots stand for experimental data [START_REF] Usbr | Design of small dams[END_REF] and the lines stand for ISPH-USAW. solver (in our case Bi-CGSTAB), but the cause of it has not been identified.

A similar behaviour was observed on a case of a schematic dam-break over an obstacle, both with an SPH formulation and with a two-phase Finite Volume formulation using the Volume of Fluid technique for the free-surface (see [START_REF] Leroy | Unified semi-analytical wall boundary conditions applied to 2-D incompressible SPH[END_REF]).

3-D cases

Laminar Poiseuille flow in a circular pipe with inflow/outflow boundaries

This case consists of a steady laminar flow through a 3-D pipe with a circular cross-section. Inflow and outflow boundaries are imposed at the extremities of the pipe. All lengths are made dimensionless by L, the radius of the crosssection. The dimensionless length of the pipe is equal to 4. The reference velocity of the flow U is the maximum velocity in the pipe, set to 1ms -1 by imposing the theoretical dimensionless velocity at the inlet:

v + = 1 -(y + -y + 0 ) 2 + (z + -z + 0 ) 2 e x (40)
where (y + 0 , z + 0 ) = (0, 0) are the dimensionless transverse and vertical coordinates of the cross-section centre. At the outlet, a zero-pressure is imposed. The Reynolds number is set to 10. The dimensionless time is defined by t + = tU L . The simulation is run until t + = 25, which corresponds to about 1.4 × 10 5

iterations with an initial dimensionless interparticular space δr + = 0.04. A steady-state is then achieved. Figure 10 shows the shape of the velocity field in the pipe at t + = 25. Figure 11 shows velocity and pressure profiles in the pipe at t + = 25. On the left, the horizontal dimensionless velocity v +

x obtained with ISPH-USAW is plotted as a function of z + along the vertical profile at the centre of the channel (x + = 2 , y + = 0). Note that the same results were obtained on a vertical profile in y + = 0 but at x + = 3.9 instead of x + = 2. Good agreement with the theoretical parabolic profile ( 40) is obtained. The dimensionless pressure is defined by:

p + = p ρU 2 /2 (41)
On the right of Figure 11, the dimensionless pressure p + obtained with ISPH-USAW is plotted as a function of x + along the horizontal profile at the centre of the channel. The agreement with the theoretical linear pressure distribution along the channel is good, the latter being given by:

p + = 8 Re 4 -x + (42)
where Re = U L ν . A small discrepancy close to the inflow boundary appears, where a homogeneous Neumann condition is imposed on the pressure. Nevertheless, the quality of the results shows that the 3-D ISPH-USAW model performs well with inflow/outflow conditions for this Reynolds number.

It should be noted that trying to increase the Reynolds number leads to numerical instability. With a Reynolds of 100 the results are as good as with a Reynolds of 10, but with a Reynolds of 1000 a blow-up of the simulation occurred after some time. This issue points out the fact that further work is necessary for high-Reynolds simulations in 3D. It is worth mentioning that increasing the Reynolds number in SPH can be challenging, as pointed out e.g. in [START_REF] Adami | A transport-velocity formulation for smoothed particle hydrodynamics[END_REF]. Finally, Figure 12 shows the evolution of the total mass of particles in the domain. The total mass tends to a constant which shows that a steady state is reached, however some mass peaks are observed with a periodic pattern. The time between two peaks is that spent by a particle to go through half of the pipe's length and the mass fluctuations in steady state are in the range of one layer of particles on the cross-section of the pipe. Figure 12: Laminar flow in a 3-D circular pipe with inflow/outflow conditions. Evolution of the total mass of particles in the computational domain using δr + = 0.04.

Preliminary application case: 3-D non-isothermal connected pipes

This case is a preliminary application case in 3-D that consists of two connected circular pipes at different temperatures. The geometry is described in the Figure 13. The reference length L is the diameter of the larger pipe. At the initial time, fluid with temperature T h is placed in a horizontal pipe with zero velocity and pressure and fluid with temperature T c (T c < T h ) is placed in a smaller inclined pipe connected to the first one with zero velocity and pressure.

As time goes by, fluid with temperature T h is injected through the left extremity of the horizontal pipe and fluid with temperature T c is injected at the highest extremity of the inclined pipe. The velocity is imposed at these inflow boundaries. In the horizontal pipe, it is imposed through a circular laminar Poiseuille profile:

v = U 1 - x -x 0 2 (L/2) 2 n 0 (43)
where x 0 is the position of the centre of the big pipe cross-section at the inlet, and n 0 = e x is the unit normal vector to that cross-section. U is the reference velocity of the flow and was set to 0.5ms -1 . The Reynolds number based on U and L was set to 10. On the other hand, in the inclined pipe the inlet velocity is imposed through

v = U 2 1 - x -x 1 2 (l/2) 2 n 1 (44)
where x 1 is the position of the centre of the small pipe cross-section at this inlet,

and

n 1 = -1 √ 3 , 0, -2 √ 3
is the unit normal vector to that cross-section. An was set to 0.162, and the Prandtl number ν K to 69. Figure 14 shows the shape of the temperature field at several instants during the simulation. a cut dam-break, a case of propagation of a solitary wave and a Creager weir, the latter involving the k -turbulence closure. The results show that the inflow/outflow conditions are able to let the fluid enter or leave the domain as prescribed, with few visible reflections in the domain. To go further, the outflow boundary conditions could be improved: the outflow boundary condition used for the wave case is quite a simple formulation for a radiative outlet and could be made more complex and more effective. On the other hand, a zero pressure was imposed at the outflow boundary on the cut dam-break, but a better condition could probably be designed. The conditions proposed in the literature for Eulerian methods could probably be applied in the framework of SPH, but this has not been tested extensively yet. To our knowledge, the question of free outlet conditions is quite problematic in the framework of Eulerian methods as well.

Finally, one case of validation in 3-D is also presented, as well as prelimi-nary results on a simple application case involving temperature effects in 3-D.

Further work in 3-D is necessary since instabilities were observed when simulating 3-D flows with high Reynolds numbers using the present formulation for inflow/outflow boundaries.

Figure 1 :

 1 Figure 1: Sketch of the different sets of entities involved in the representation of the USAW boundary conditions.

  are two main requirements for the imposition of open boundaries in an ISPH model: an algorithm to let particles enter and leave the domain, and the correct imposition of open boundary conditions on the fields, in particular on the pressure.5.1. Particles creation/destructionRegarding the algorithm for particles creation/destruction, the technique proposed by Kassiotis et al. in[START_REF] Kassiotis | Semi-analytical conditions for open boundaries in smoothed particle hydrodynamics[END_REF] is used. The idea is to let the masses of the inlet/outlet vertex particles v ∈ V i/o evolve over time as a function of the desired ingoing/outgoing mass flux through the inlet/outlet segments s ∈ S i/o directly connected to v. The vertex particles are then used to create/delete fluid particles. The mass evolution should not introduce any perturbations in the flow, so care must be taken that its evolution is smooth. The time-derivative of the mass, denoted by ṁn v , is determined by the Eulerian velocity u s imposed at the open boundaries:

Figure 2 :

 2 Figure 2: Sketch of the process of particles creation with vertices v and segments s at an inflow boundary [15]: a) the vertex masses grow due to the ingoing flux; b) their mass has reached the maximum threshold; c) new free particles are released and the vertex masses become negative (dashed line); d) the vertex masses start growing again.

Figure 3 :

 3 Figure 3: Open boundary technique based on the USAW boundary conditions: sketch of the notations for the computation of the fraction of segment area β a,v attributed to a vertex at the outlet [15]. (a) 2-D case; (b) 3-D case. The coefficient associated to the vertex v 0 is proportional to the red area so that the largest amount of mass goes to the closest vertex particle.

a

  boundary with prescribed velocity, or to enter the domain through a boundary with prescribed pressure. In other words, inlet and outlet conditions can be handled by any open boundary at the same time. This is necessary in case of a prescribed recirculation close to an inlet. Then, at boundaries where the velocity is prescribed a Dirichlet condition is also set on the temperature, and a homogeneous Neumann condition is set on the pressure. On the other hand, at boundaries where the pressure is prescribed a homogeneous Neumann condition is set on the velocity and temperature fields. Here the open boundary conditions on k and are not detailed: they are the same as for the temperature. The imposition of these open boundary conditions is done in a similar way as for wall boundaries in [12, 33]: Dirichlet conditions are imposed at the vertex particles whereas the Neumann conditions are imposed through the segments by setting the boundary terms of the Laplacian operators. Compatible conditions are then deduced and imposed on the complementary entities (i.e. segments in case of a Dirichlet condition, vertex in case of a Neumann condition).

  where S pp is the set of open boundary segments with prescribed pressure. The lines of the matrix corresponding to inlet/outlet vertex particles are removed.

Figure 4 :

 4 Figure 4: Dam-break with an outflow condition: velocity field shape obtained with the present ISPH-USAW model. Comparison with the free-surface shape of a non-cut dam-break (black dots).

Figure 5 :

 5 Figure 5: Sketch of the geometry of the case of propagation of a solitary wave on a slope.

Figure 7 :

 7 Figure 7: Creager weir geometry, from [37] after US Army Corps of Engineers (1990). Overall view (a), vertical profile in the indented area (b) and zooming in on the area surrounded by a dashed line (c).

Figure 8 :

 8 Figure 8: Creager weir case. Velocity field for a flow on a Creager weir for H/H d = 0.5 (top) and H/H d = 1.33 (bottom) with ISPH-USAW (flow from left to right).

Figure 9 :

 9 Figure 9: Creager weir case. Free surface (left) and pressure (right) on weir for H/H d = 0.5 (red) and H/H d = 1.33 (blue). Both are mean values over 3s once the flow has reached a steady state. The dimensionless pressure is p + = p/ρgd.The dots stand for experimental data[START_REF] Usbr | Design of small dams[END_REF] and the lines stand for ISPH-USAW.

Figure 10 :

 10 Figure 10: Steady laminar flow in a 3-D circular pipe with inflow/outflow conditions. Shape of the velocity profile obtained with ISPH-USAW using δr + = 0.04 at t + = 25 (flow from right to left).

Figure 11 :

 11 Figure 11: Steady laminar flow in a 3-D circular pipe with inflow/outflow conditions. Dimensionless velocity profile along the vertical line at the centre of the channel (left) and dimensionless pressure profile along the horizontal line at the centre of the channel (right). Comparison of the results obtained with ISPH-USAW using δr + = 0.04 with the theoretical velocity and pressure profiles.

Figure 13 :

 13 Figure 13: Laminar flow in two connected pipes at different temperatures. Sketch of the geometry.

  outflow boundary condition is imposed at the right extremity of the horizontal pipe: the pressure is imposed to zero and a homogeneous Neumann condition is imposed on the A homogeneous Neumann condition also imposed on the temperature at solid walls. The Grashoff number βg∆T L 3 ν 2

Figure 14 :

 14 Figure 14: Laminar flow in two connected pipes at different temperatures. Shape of the temperature field at several times.

  In this work a new formulation for inflow/outflow conditions in an ISPH model was proposed. The boundary conditions rely on the unified semi-analytical technique, which makes it possible to accurately prescribe Neumann or Dirichlet boundary conditions on the fields. The fluid enters and leaves the domain through the mass variation of boundary particles and their creation/destruction. Validation was performed on several 2-D cases presenting open boundaries:

Table 1 :

 1 Values of the k -model constants[START_REF] Launder | Mathematical models of turbulence[END_REF]

	κ	C µ	C 1	C 2	σ k	σ
	0.41 0.09 1.44 1.92 1.0 1.3

Table 2 :

 2 Parameters for the Creager weir numerical simulations

	H/H d		0.5	1.33
	H(m)		0.15	0.399
	q(m 3 /s/m)	0.115	0.560
	d(m)		0.547	0.78
	U (m/s) F r = U/ √	gd	0.21 0.09	0.72 0.26
	u * (m/s)	0.0089 0.0266
	Domain length after the weir (m)	1	2.6

Acknowledgements

This work was partly funded by the French Research Agency (CIFRE grant agreement # 2011-0264).


Agnès Leroy
email: agnes.leroy@edf.fr

Damien Violeau
email: damien.violeau@edf.fr

Martin Ferrand

Antoine Joly

Buoyancy modelling with incompressible SPH for laminar and turbulent flows

Keywords: SPH, incompressible, boundary conditions, turbulence, buoyancy, temperature

come

Introduction

This work aims at modelling buoyant, laminar or turbulent flows, using a 2D Incompressible Smoothed Particle Hydrodynamics (ISPH) model where wall boundary conditions are imposed through an appropriate and accurate technique. The word 'buoyant' herein refers to flows where the density varies in space and time due to an active scalar like temperature or salinity. For simplicity, we will always refer to 'temperature' and denote the scalar field by T .

Buoyancy plays an important part in many industrial and environmental flows.

Modelling its effects through a Lagrangian method presents the advantage of avoiding artificial diffusion, since the advection is inherently represented by the particle displacement.

The buoyancy effects are modelled here through the Boussinesq approximation coupled to a heat equation, so that density variations only act through an additional term in the gravity force of the momentum equation. As a consequence, the density is considered constant in all governing equations, which makes it possible to apply an incompressible algorithm to compute the pressure field from a Poisson equation. This approach has been widely used in CFD (Computational Fluid Dynamics) with mesh-based methods (see e.g. [START_REF] Reddy | The finite element method in heat transfer and fluid dynamics[END_REF] regarding the Finite Elements method or [START_REF] Oliveira | An improved piso algorithm for the computation of buoyancy-driven flows[END_REF] with the Finite Volume method), but is still rather new in the SPH publications. To our knowledge, the only similar works so far have been proposed by Szewc et al. [START_REF] Szewc | Modeling of natural convection with Smoothed Particle Hydrodynamics: Non-Boussinesq formulation[END_REF] and by Ghasemi et al. [START_REF] Ghasemi | 2D numerical simulation of density currents using the SPH projection method[END_REF], with classical but rather crude approaches for prescribing wall boundary conditions (ghost particles and dummy particles respectively). Based on our previous work [START_REF] Leroy | Unified semi-analytical wall boundary conditions applied to 2-D incompressible SPH[END_REF], we extend here the unified semi-analytical wall boundary conditions (hereinafter denoted USAW) to the present SPH model of buoyancy. This technique has proved its ability to work in the framework of ISPH. The main advantage of working with the USAW boundary conditions, besides their accuracy on complex wall geometries in arbitrary space dimension, is that they make it possible to exactly impose arbitrary wall boundary conditions on the temperature. Indeed, we propose a unified framework to prescribe isothermal (Dirichlet) or imposed heat flux (Neumann) wall boundary conditions in ISPH.

As in [START_REF] Leroy | Unified semi-analytical wall boundary conditions applied to 2-D incompressible SPH[END_REF], the proposed model is combined to a Reynolds-Averaged Navier-Stokes (RANS) approach to treat turbulent flows. The k -turbulence model is used, as it is today a standard approach for industrial and environmental flows. The second novelty proposed here is to apply the latter model to SPH with buoyancy. This is achieved by means of usual tools for grid-based methods, i.e. introducing a buoyant term in the k -equations (see e.g. [START_REF] Rodi | Turbulence models and their applications in hydraulics[END_REF]). The latter term is proportional to the temperature gradient, here computed with the SPH tools.

Again, the USAW technique is used to compute the relevant boundary terms.

To illustrate the ability of the present model to predict buoyant flows, a theoretical case is presented (laminar heated Poiseuille flow), where excellent agreement with the theoretical solution is obtained. Several benchmark cases are then proposed: a lock-exchange flow, two laminar and one turbulent flow in differentially heated cavities, and finally a turbulent heated Poiseuille flow. Comparisons are provided with a Finite-Volume (FV) approach using an open-source industrial code.

Governing equations and modelling choices

The system of equations to be solved is composed of the incompressible Reynolds-Averaged Navier-Stokes (RANS) equations coupled to a heat equation and to the k -turbulence closure. As explained in the Introduction, the formalism of heat transfer was chosen, but the reasoning applies to other active scalars like salinity. The Boussinesq approximation is used to account for density variations so that the system reads:

                                     ∇ • v = 0 dk dt = P + G -+ 1 ρ ∇ • (µ k ∇k) d dt = k (C 1 P + C 3 G -C 2 ) + 1 ρ ∇ • (µ ∇ ) dv dt = - 1 ρ ∇p + 1 ρ ∇ • (µ E ∇u) + g (1 -β(T -T 0 )) dr dt = v dT dt = K E ∇ 2 T (1)
In this system, r is the particle position, v is the Lagrangian velocity and u is the Eulerian velocity. Both velocities are equal in our SPH model except for wall boundary particles as we will see in Section 3.3. t is the time, ρ is the density (which remains constant due to the Boussinesq approximation), k is the turbulent kinetic energy field, is its dissipation rate, µ is the dynamic molecular viscosity (which may be a function of the temperature), µ T is the dynamic eddy viscosity, µ E = µ + µ T is the effective (i.e. total) dynamic viscosity. We also define the two variables

µ k = µ + µ T σ k and µ = µ + µ T σ . σ k , C 1 , C 2 and
σ are model constants described in Table .1. P is the production of turbulent kinetic energy and G is a buoyancy production/destruction term. C 3 is set to 1 if G ≤ 0 and 0 otherwise. P is calculated according to a mixed linear-quadratic model [START_REF] Guimet | A linearised turbulent production in the k-model for engineering applications[END_REF]:

P = min C µ kS, ν T S 2 (2)
where S = √ 2S : S is the scalar mean rate-of-strain. G is a buoyancy production/destruction term calculated as [START_REF] Rodi | Turbulence models and their applications in hydraulics[END_REF]:

G = βK T ∇T • g (3)
g is the gravity field (of magnitude g = 9.81m 2 s -1 ), p = p + 2 3 k with p the pressure, β is the coefficient of thermal expansion, T is the temperature field, T 0 is the mean temperature. K E = K +K T is the effective thermal diffusivity, with K the molecular thermal diffusivity and K T = ν T P r T , P r T being the turbulent Prandtl number, taken as 0.85 (see e.g. [START_REF] Rodi | Turbulence models and their applications in hydraulics[END_REF]). We also define ν = µ ρ and ν T = µ T ρ . ν T is calculated as a function of k and of as usual [START_REF] Launder | Mathematical models of turbulence[END_REF]:

ν T = C µ k 2 (4)
where C µ is a constant defined in Table .1. Note that the eddy viscosity is imposed equal to zero at the walls. In case of a laminar flow, ν T is set to zero and the k and equations are not solved.

Buoyancy modelling with incompressible SPH and the USAW boundary conditions

Space discretisation

In all this paper we work in 2D. We assume the reader is familiar with the standard SPH method (see [START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF] for more details). With the SPH method and the USAW boundary conditions the space discretisation is done through different sets of particles and boundary elements called segments. Fluid particles which do not belong to a boundary are called free particles a ∈ F , whereas particles belonging to a boundary are called vertex particles v ∈ V . The latter are linked together by segments s ∈ S, which compose a mesh of the boundary [START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF][START_REF] Leroy | Unified semi-analytical wall boundary conditions applied to 2-D incompressible SPH[END_REF] (see Figure 1). The set of all fluid particles is noted P = F ∪ V . In this framework, the discrete SPH interpolation of a field A at particle a with position r a reads:

[A] γ a = 1 γ a b∈P V b A b w ab (5)
where V b = m b ρ is the volume of particle b, m b being its mass, and w ab = w(r a -r b ), w being the SPH kernel. In general the value of a field A at a particle b is denoted by A b . γ a is the wall renormalisation factor, defined as in [START_REF] Kulasegaram | A variational formulation based contact algorithm for rigid boundaries in two-dimensional SPH applications[END_REF][START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF]:

γ a = Ω∩Ωa w(r a -r )d n r ( 6
)
where Ω is the fluid domain, Ω a is the compact support of the kernel at particle a and n is the space dimension. γ a is computed according to an analytical formula [START_REF] Leroy | Unified semi-analytical wall boundary conditions applied to 2-D incompressible SPH[END_REF]. We define two SPH gradient operators of the field A which will be used for different purposes:

v 1 s v 2
G γ,+ a {A b }= 1 γ a b∈P V b (A a +A b ) ∇w ab - 1 γ a s∈S (A a +A s ) ∇γ as (7) G γ,- a {A b } = - 1 γ a b∈P V b A ab ∇w ab + 1 γ a s∈S A as ∇γ as (8)
where A ab = A a -A b , A as = A a -A s and ∇γ as is the contribution of segment s to the gradient of γ a , defined as:

∇γ as = ∂Ωs∪Ωa w(r a -r )n s d n-1 r (9)
In this formula, ∂Ω s is the boundary area spanned by segment s and n s is the inward unit normal to the wall at segment s (Figure 1). ∇γ as is computed according to an analytical formula [START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF]. An SPH divergence of a field A is defined as:

D γ a {A b } = - 1 γ a b∈P V b A ab • ∇w ab + 1 γ a s∈S A as • ∇γ as (10)
The SPH Laplacian operator used in this work is the one proposed by Ferrand et al. [START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method[END_REF] which reads:

L γ a {B b , A b } = 1 γ a b∈P V b (B a + B b ) A ab r 2 ab r ab • ∇w ab - 1 γ a s∈S [B s (∇A) s + B a (∇A) a ] • ∇γ as (11)
where B is a (variable) diffusion coefficient for the field A and r ab = |r ab |. In case A is a vector, the Laplacian will be noted L γ a {B b , A b } and in case B = 1 the Laplacian will be noted L γ a {A b }. Later on, we will sometimes use the following notation for the boundary term of [START_REF] Kulasegaram | A variational formulation based contact algorithm for rigid boundaries in two-dimensional SPH applications[END_REF]:

L bound a {B b , A b } = - 1 γ a s∈S [B s (∇A) s + B a (∇A) a ] • ∇γ as (12)

Time discretisation

In the turbulent case, k a and a are calculated at the beginning of the iteration through a semi-implicit time-scheme making up the SPH form of the standard buoyant k -model (second and third lines of (1)):

k n+1 a -k n a δt = P n a + G a -n a k n+1 a k n a + 1 ρ L γ a {µ k,b , k n b } (13) n+1 a -n a δt = n a k n a C 1 P n a + C 3 G a -C 2 n+1 a + 1 ρ L γ a {µ ,b , n b } (14)
where the superscript n represents the iteration number and δt is the time step.

In [START_REF]Code_Saturne 3.0.0 Theory Guide[END_REF] and ( 14) the dissipative terms are treated implicitly in order to avoid negative values of k and (recall ρ is a constant, so we omit the particle subscript for this variable). P a is calculated according to [START_REF] Oliveira | An improved piso algorithm for the computation of buoyancy-driven flows[END_REF], where S a is computed with the following SPH approximation of the mean rate-of-strain tensor:

S a = 1 2 G γ,- a {u b } + G γ,- a {u b } T (15)
According to (3), G a is modelled by:

G a =        βC µ P r T k n a k n+1 a n a G γ,- a {T b } • g if G γ,- a {T b } • g ≤ 0 βC µ P r T (k n a ) 2 n a G γ,- a {T b } • g otherwise (16)
Again, G a is semi-implicited in case it is negative in order to avoid negative values of k and . After the computation of k, and ν T (using ( 4)), the spacetime discretisation of (1) follows an ISPH predictor-corrector scheme with a pressure Poisson equation, as in [START_REF] Leroy | Unified semi-analytical wall boundary conditions applied to 2-D incompressible SPH[END_REF], with an additional SPH diffusion equation for the temperature:

                                                     v * a -v n a δt = 1 ρ L γ a {µ E,b , u n b } -gβ(T n a -T 0 ) L γ a {p * n+1 b } = ρ δt D γ a {v * b } v n+1 a -v * a δt = - 1 ρ G γ,+ a {p * n+1 b } r n+1/2 a = r n a + δt 2 v n+1 a T n+1 a -T n a δt = L γ a {K E,b , T n b } rn+1/2 a = r n+1/2 a + δr a ṽn+1 a = v n+1 a + G γ,- a {v n+1 b } • δr a r n+1 a = rn+1/2 a + δt 2 ṽn+1 a ( 17
)
where v * is the predicted velocity field, p * = p + ρgz is the dynamic pressure and δr a is a particle shifting defined as in [START_REF] Leroy | Unified semi-analytical wall boundary conditions applied to 2-D incompressible SPH[END_REF]:

δr a = -0.35h 2 G γ,+ a {1} (18)
The particle shifting has been used to modify the velocity through a first-order

Taylor expansion (7th line of ( 17)). The fields u, k, and T are corrected the same way.

The last three steps of ( 17) correspond to a stabilising procedure (see [START_REF] Lind | Incompressible smoothed particle hydrodynamics for free-surface flows: A generalised diffusion-based algorithm for stability and validations for impulsive flows and propagating waves[END_REF][START_REF] Leroy | Unified semi-analytical wall boundary conditions applied to 2-D incompressible SPH[END_REF]).

Note that system [START_REF] Härtel | Analysis and direct numerical simulation of the flow at a gravity-current head. Part 1. Flow topology and front speed for slip and no-slip boundaries[END_REF] was written here for confined flows since all the validation cases considered in this paper are confined. This is why the resolution is done with the dynamic pressure instead of the total pressure, which yields more accurate results. In case of free-surface flows this can not be done because of the pressure Dirichlet condition on the free-surface, so that the term involving gravity in the first line of ( 17) becomes -g (β(T n a -T 0 ) -1) and p * is replaced by p in the second and third lines of (17).

Imposition of wall boundary conditions on u, v, k and

The present section summarises our wall boundary conditions. Our technique is based on an analogy with Finite Volume (FV) and was validated in [START_REF] Leroy | Unified semi-analytical wall boundary conditions applied to 2-D incompressible SPH[END_REF], especially for the k -turbulent model. The imposition of pressure wall boundary conditions is done as in [START_REF] Leroy | Unified semi-analytical wall boundary conditions applied to 2-D incompressible SPH[END_REF] and is not described in this paper. A Dirichlet boundary condition is imposed on v: the Lagrangian velocity of the walls is imposed so that the model includes the treatment of forced wall movement through the velocities of the vertex. The velocities of the segments are then defined by:

v s = 1 N s Vs v v (19)
where V s is the set of vertices linked to s and N s its size (in 2D, N s = 2).

Non-homogeneous Neumann boundary conditions are imposed on u, k and .

In such cases, the prescription of wall boundary conditions is done by imposing both the flux and the value of the field at the wall. Therefore, the compatibility of the fields values and fluxes at the wall must be ensured. This holds for the temperature field as well (see section 3.4). The Neumann wall boundary conditions are applied through the surface term of the Laplacian operator [START_REF] Lind | Incompressible smoothed particle hydrodynamics for free-surface flows: A generalised diffusion-based algorithm for stability and validations for impulsive flows and propagating waves[END_REF] like in mesh-based methods, whereas the Dirichlet boundary conditions are imposed at the vertex particles which are involved in the summations over b ∈ P in the Laplacian, gradient and divergence operators. Thus, in the aforementioned equations the particles a belong to F . In the subsequent subsections, the notation A ab = A a -A b will be widely applied as in Section 3.1.

Wall boundary condition on u

A non-homogeneous Neumann condition is applied to the Eulerian velocity field:

1 ρ L surf a {µ E,b , u b } = - 2 γ a ρ s∈S µ E,a ∂u ∂n a • ∇γ as ( 20
)
where a is a fictitious point placed at ras 2 . In the laminar case, the velocity distribution near the wall is almost linear, thus:

µ E,a ∂u ∂n a • n s = µ v as • t as δr as t as (21)
where

     t as = v as -(v as • n s ) |v as -(v as • n s ) | δr as = max(r as • n s , δr) (22)
with δr the initial interparticular space.

In the turbulent case, a two layers wall function is used for the velocity near the wall. The Neumann condition reads:

µ E,a ∂u ∂n a • n s = u 2 * ,a t as (23)
where u * ,a is the friction velocity at the wall seen by particle a, which is a solution of:

       v as • t as u * ,a = y + a if y + a ≤ y + lim v as • t as u * ,a = 1 κ ln δr as u * ,a ν + 5.2 if y + a > y + lim ( 24
)
where y + a = δrasu k,a ν , y + lim = 1 κ and κ is the von Kármán constant (see Table .1). The second line of (24) is solved through an iterative process.

On the other hand, the velocity at the vertex particles is left to evolve according to the viscous term:

u n+1 v = u n v + δt 1 ρ L γ v {µ E,b , u b } (25)
but its normal component is imposed to be equal to zero by projecting u n+1 v along the tangent to the wall.

Wall boundary condition on k

For the turbulent kinetic energy, a homogeneous Neumann condition is applied:

L bound a {µ k,b , k b } = 0 (26)
A compatible Dirichlet boundary condition on k is imposed at all vertex particles v:

k v = 1 N v s∈Sv k s , k s = 1 α s b∈F V b k b w sb (27)
where S v is the set of segments linked to v, N v is its size and α s is the Shepard filter:

α s = b∈P V b w sb (28)

Wall boundary condition on

Here again, a non-homogeneous Neumann condition is applied according to the theory of turbulent boundary layer:

L bound a {µ ,b , b } = 4C µ σ γ a s∈S k 2 a δr as |∇γ as | (29)
A compatible Dirichlet boundary condition is imposed on at all vertex particles v (see [START_REF] Leroy | Unified semi-analytical wall boundary conditions applied to 2-D incompressible SPH[END_REF] for more detail):

v = 1 N v s∈Sv s , s = 1 α s b∈F V b b + 4C 3/4 µ k 3/2 b κδr sb w sb (30)
L bound a {K E,b , T b } = - 2 γ a s∈S Q T s • ∇γ as (31)
where

Q T s = K ∂T ∂n s
is the heat flux at the wall (note that since ν T = 0 at the wall, this also holds for K T ). The values of the heat flux and of the temperature at the wall depend on the type of wall boundary conditions.

Neumann wall boundary condition

In this case the values of Q T s are imposed. A compatible Dirichlet condition is prescribed at the vertex particles, depending on the imposed flux:

T v = 1 N v s∈Sv T s , T s = 1 α s b∈F V b T b - Q T s • n s K δr sb w sb (32)
The previous choice for the temperature of the wall segments T s is inspired from the conditions described above for the k -model. It directly stems from a first-order Taylor expansion of the temperature near the wall.

Dirichlet wall boundary condition

In this case the value of the temperature is prescribed at the vertex particles.

The Q T s are imposed in a consistent way with the Dirichlet condition, so that the boundary term (31) reads:

L bound a {K E,b , T b } = - 2 γ a s∈S T * ,a u k,a |∇γ as | (33)
where u k,a is a friction velocity (which should not be confused with u * ,a defined above), and T * ,a =

Q T s •ns u k,a
by definition. In laminar mode, a linear temperature wall function is applied:

u k,a T * ,a = Q T s • n s = K T a -T s δr as (34)
with T s imposed through the Dirichlet condition, whereas in turbulent mode, u k,a is defined as:

u k,a = C 1/4 µ k 1/2 a ( 35
)
On the other hand, T * ,a is defined as:

T * ,a = T a -T s T + a ( 36
)
where T + a is computed through a three-layers wall function (see [START_REF]Code_Saturne 3.0.0 Theory Guide[END_REF]):

               T + a = P r y + a if y + a < y + 1 T + a = a 2 - P r T 2a 1 y + a 2 if y + 1 ≤ y + a < y + 2 T + a = P r T κ ln y + a + a 3 if y + a > y +
where y + a is defined as in (24), while

                                     y + 1 = a 4 P r 1/3 y + 2 = a 4 κ P r T a 1 = P r T a 4 a 2 = 15P r 2/3 a 3 = 15P r 2/3 - P r T 2κ 1 + ln a 4 κ P r T a 4 = 1000 (38)

Validation on laminar flows

The 5th order Wendland kernel [START_REF] Wendland | Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree[END_REF] was used for all the simulations with a smoothing length h = 2δr (recall δr is the initial interparticular spacing). In what follows the non-dimensional variables are denoted with a + superscript.

Hereafter, the present buoyant incompressible SPH model with our boundary conditions will be referred to as ISPH-USAW. The FV results were obtained

with the Code_Saturne open-source software [START_REF] Archambeau | Code_Saturne: a Finite Volume code for the computation of turbulent incompressible flows -industrial applications[END_REF].

Laminar Poiseuille flow

In order to check that the boundary conditions on the temperature are properly imposed by the method described in section 3.4, two configurations of a 2D laminar plane Poiseuille flow were tested. A schematic description of their geometries is provided in Figure 2. Our coordinates are denoted (x, z), the origin of the vertical axis being in the middle of the two walls, distant from 2L. In the first case (denoted TT), constant temperatures T 1 and T 2 were imposed on the lower and upper walls, respectively. In the second case (denoted QT), a constant heat flux Q T was imposed through the upper wall, while the lower wall remained isothermal at the temperature T 1 . The flow, with bulk velocity U , is driven by a constant volumic force.

In the case TT the flow is governed by three dimensionless numbers: while in the case QT the flow is governed by four dimensionless numbers:

             P r = ν K Re = U L ν Gr = βg∆T L 3 ν 2 (39) (
                     P r = ν K Re = U L ν Gr = βg∆T L 3 ν 2 N u = LQ T 0 K∆T (40)
where ∆T = max(T ) -min(T ). In the two cases the Prandlt number was set to 1 and the Reynolds number to 50. In the QT case, the Nusselt number was set to 0.5 and the Grashoff number was set to 196, while in the TT case the Grashoff number was set to 98.

The theoretical expressions of the dimensionless temperature and dynamic pressure p * = p+ρgz for the two cases are presented in the Table .2 as functions of Re, Gr, N u and the dimensionless coordinates x + = x L , z + = z L . In both cases the velocity field is that of the ordinary Poiseuille flow. The dimensionless fields were defined as p * + = p * ρβ∆T gL and T + = T -T1

∆T . The simulations were done with 902 particles ( δr L = 0.05). Figure 3 shows the vertical profiles of dimensionless temperature and dynamic pressure obtained with ISPH-USAW.

It can be observed that an excellent agreement with the theory is obtained in both cases. In the case QT, the vertical variation of the fields does not depend on the Grashoff number. Thus, different values of the heat flux should yield the same results, which was checked with ISPH-USAW and gave similar errors between the model and the theory in all cases (with a Grashoff number up to 1960). With the chosen discretisation, the maximum relative error was of the order of 1%, on the temperature and on the pressure. A convergence study was done for the case TT. The convergence plot is shown on the Figure 4, where the L 2 error on the temperature was computed through:

L 2 error = 1 V tot b∈P V b T -T theo T 2 2 (41)
where T theo is the theoretical temperature profile (linear profile) and V tot is the total volume of the particles (V tot = b∈P V b ). The order of convergence is close to 2 for lower discretisations, and decreases for finer simulations because the discretisation error threshold is being reached. These results show that the imposition of Dirichlet and Neumann boundary conditions on the temperature is properly done with our method.

Differentially heated square cavity

The second validation case consists of a laminar flow in a differentially heated square cavity of size L that was studied in [START_REF] Szewc | Modeling of natural convection with Smoothed Particle Hydrodynamics: Non-Boussinesq formulation[END_REF] with the SPH method. The left and right walls are isothermal, the cold wall being on the right, and the upper and lower walls are adiabatic. The molecular Prandlt number is P r = 0.71 and three values of the Rayleigh number Ra = P r × Gr were tested, namely 10 3 , 10 4 and 10 5 (the Grashoff number is defined as in the 3rd line of (40)). We define the dimensionless variables x + = x L , z + = z L , u + = uL ν , T + = T T h where T h is the temperature of the hot wall. A discretisation of 160 × 160 particles was used for the SPH simulations.

N u s = L ∆ T |G γ,- s {T b } • n s | (42)
Figure 7 shows the evolution of the local Nusselt number along the cold wall for the three values of Rayleigh number. It appears that the behaviour is globally well predicted. However, a discrepancy occurs on the top of the curve Ra = 10 5 , which corresponds to the top right corner of the flow in Figure 5, where temperature gradients are rather high. It is a consequence of a lack of accuracy of the G γ,-SPH gradient operator, used to compute the Nusselt number, since our temperature profiles are still in very good agreement with FV in this area.

Differentially heated lid-driven cavity

A differentially heated lid-driven cavity at Re = U L ν = 1000 was tested, L being the size of the cavity and U the velocity of the lid. The flow is driven by A discretisation of 160×160 particles was used for the SPH simulation. Figure 8 shows the shape of the temperature and velocity fields after convergence for a Rayleigh number of 10 5 . The results are compared to FV using a discretisation of 512 × 512 cells. We define the dimensionless variables

x + = x L , z + = z L , u + = u U , T + = T T h .
Figures. 9 and 10 show the velocity and temperature profiles along x + = 1/2 and z + = 1/2. Satisfactory agreement was obtained with FV for the velocity, while the temperature shows a very nice agreement.

The discrepancy in the velocity profile is due to the difference in the space discretisation, which is here significantly lower with the present SPH model.

Lock-exchange

The fourth validation case consists of a symmetric lock-exchange flow in a rectangular cavity of height 2L and width 30L. This case was studied in [START_REF] Ghasemi | 2D numerical simulation of density currents using the SPH projection method[END_REF] shown in Figure 13 were obtained with a symmetric operator for the pressure gradient [START_REF] Launder | Mathematical models of turbulence[END_REF], which reproduced the vortices the interface of the two fluids compared to the DNS results. Figure 14 shows the results obtained on this case with the antisymmetric SPH operator for the pressure gradient [START_REF] Guimet | A linearised turbulent production in the k-model for engineering applications[END_REF].

The numerical diffusion is higher than with the symmetric operator and the instabilities do not develop as well. For all other test-cases in the present work, an antisymmetric operator ( 7) was used since it conserves linear momentum, but no significant differences were observed when using a symmetric operator.

Infinite channel with a complex geometry

This case was designed in order to show that the model is able to simulate flows with complex wall geometries. It consists of an infinite horizontal channel

Validation on turbulent flows

The validation of the k -turbulence model with buoyancy was done on two classical validation cases: a turbulent plane Poiseuille flow with two isothermal walls and a rectangular differentially heated cavity.

Turbulent plane Poiseuille flow with two isothermal walls

This case consists of a turbulent flow between two parallel, infinite vertical isothermal walls. The temperature of the left wall is set to T h and that of the right wall to T c , with T h > T c . The flow is thus vertical, driven by a prescribed pressure gradient (the friction velocity is imposed) and by the temperature difference between the walls. The buoyancy force acts upward near the hot wall and downward near the cold wall. The friction Reynolds number, Re * = u * L ν is set to 150, where L is the half-width of the channel and u * is the friction velocity, which was set to 1ms -1 through the imposition of an upward volumic force:

f = u 2 * L -gβ b∈F V b (T b -T 0 ) b∈F V b (43)
where

T 0 = T h +Tc 2 .
The molecular Prandlt number is set to 0.71 and the Grashoff number to Gr = 9.6 × 10 5 . The following dimensionless variables are defined: and are visible close to the walls. They seem to be due to the differences in the imposition of the wall boundary conditions on those fields between FV and ISPH-USAW.

u + = u u * , x + = x L , T + = T T h , k + = k u 2 * and + = L

Differentially heated rectangular cavity in turbulent mode

The last validation case consists of a differentially heated rectangular cavity of aspect ratio 4. The geometry of the case is described in Figure 21. The SPH and FV is provided in Figure 21. Figure 22 shows the temperature profiles along x + = 0.125 and z + = 1. A satisfactory agreement is observed between ISPH-USAW and FV. Both methods present reasonable agreement with the DNS results, although some differences are observed, which were expected since a 2D RANS k -model is not meant to perfectly reproduce 3D DNS data.

Conclusions

In this work, a buoyancy model based on the Boussinesq approximation was introduced in an existing ISPH model with a RANS k -turbulence closure.

The production/destruction term due to buoyancy is taken into account in the k -equations. The wall boundary conditions are imposed through the unified semi-analytical technique, which makes it possible to impose arbitrary wall boundary conditions on the fields. In particular, the wall boundary condition imposed on the temperature can be of Dirichlet or Neumann type so as to model isothermal or adiabatic walls, or to impose an arbitrary heat flux through a wall.

The results obtained on several configurations of a 2D laminar plane Poiseuille flow show that the wall boundary conditions are properly imposed on the temperature. Several other laminar flows were modelled: a differentially-heated cavity, a differentially heated lid-driven cavity and a lock-exchange. Very good agreement was obtained with mesh-based methods in all cases. The only remaining issue is that the computation of the local Nusselt number is still inaccurate for high temperature gradients, due to a lack of accuracy of the SPH gradient operators.

The validation in turbulent mode was done on a 2D plane Poiseuille flow with two isothermal walls and a differentially-heated rectangular cavity. Comparisons with FV show that the boundary conditions in turbulent mode are satisfactorily imposed in the SPH model, in spite of some differences on k and compared to FV close to the wall.

Since buoyancy plays an important part in many industrial and environmental flows, the interest of its modelling is obvious. Modelling this process through a Lagrangian method presents the advantage of avoiding artificial diffusion, which makes SPH an appealing approach for this purpose. In view of this, further work should concern the extension of the present model to 3D, which will require GPU programming.

One should note that additional benefits could be obtained with LES (Large Eddy Simulation). The RANS model used here is relevant for performing simple simulations in the Industry, but fails to predict correlations between momentum and scalar concentration, which can be useful in some specific situations. Since SPH has already been applied to LES simulations with promising results (see e.g. [START_REF] Issa | Numerical assessment of the Smoothed Particle Hydrodynamics gridless method for incompressible flows and its extension to turbulent flows[END_REF][START_REF] Rogers | SPH modeling of breaking waves[END_REF][START_REF] Hu | Simulating 3D turbulence with SPH[END_REF][START_REF] Mayrhofer | Large eddy simulation with SPH: Mission impossible?[END_REF]), the present buoyancy model is expected to provide an interesting approach under the LES framework. This would of course avoid using the k -model, thus giving a better prediction of heat turbulent fluxes. Such an attempt is let to further investigation.

T + = 1 2 (1 + z + ) T + = -N u (1 + z + ) p * + = z +2 4 -2 Re Gr x + p * + = z + 1 2 -N u 1 + z + 2 - 2Re Gr x +

Figure 1 :

 1 Figure 1: Sketch of the different entities involved in the space discretisation with SPH and USAW wall boundary conditions.

3. 4 .

 4 Imposition of wall boundary conditions on the temperature We will now extend the ideas of Subsection 3.3 to the temperature. Again, what follows was inspired by the FV technique. Wall boundary conditions on the temperature may be of Neumann or Dirichlet type. In either case, their prescription is done by imposing both a heat flux and a temperature value at the wall. The compatibility of the temperature values and of the heat flux at the wall should be ensured, which is explained in the next two subsections. The boundary part of the diffusion term in the heat equation reads:

Figure 2 :

 2 Figure 2: Laminar Poiseuille flow: sketch of the two configurations.

Figure 3 :

 3 Figure 3: Laminar Poiseuille flow: profiles of dimensionless temperature (left) and dynamic pressure (right) on the vertical section of the channel obtained with ISPH-USAW, compared to the theoretical solutions of Table .2.

  Figure 3: Laminar Poiseuille flow: profiles of dimensionless temperature (left) and dynamic pressure (right) on the vertical section of the channel obtained with ISPH-USAW, compared to the theoretical solutions of Table .2.

Figure 4 :

 4 Figure 4: Laminar Poiseuille flow: convergence study concerning the error on the temperature field.

Figure 5 :

 5 Figure 5: Differentially heated square cavity at Ra = 10 5 . Shape of the temperature (right) and velocity (left) fields obtained with ISPH-USAW (top) and FV (bottom) after convergence.

Figure 6 :

 6 Figure6: Differentially heated square cavity at Ra = 10 5 . Profiles of velocity (left) and temperature (right) in x + = 1/2 and z + = 1/2 obtained with ISPH-USAW and FV after convergence. The horizontal profiles are also compared to the ones obtained by Wan et al.[START_REF] Wan | A new benchmark quality solution for the buoyancy-driven cavity by discrete singular convolution[END_REF] with the discrete singular convolution method.

Figure 7 :

 7 Figure 7: Differentially heated square cavity. Evolution of the Nusselt number along the cold wall of the cavity at Ra = 10 3 , 10 4 and 10 5 . Comparison of the results obtained with ISPH-USAW, FV and discrete singular convolution [16] after convergence.

Figure 11

 11 Figure 11 shows the repartition of local Nusselt number (computed according to (42)) along the upper and lower walls of the cavity. As observed in the previous test-case, the formula (42) used to compute the local Nusselt number lacks accuracy for high temperature gradients (left part of the dotted curve in Figure 11, i.e. top left angle of the flow in Figure 8).

Figure 8 :Figure 9 :Figure 10 :Figure 11 :

 891011 Figure 8: Differentially heated lid-driven cavity. Shape of the temperature (right) and velocity (left) fields obtained with ISPH-USAW (top) and FV (bottom) after convergence.

Figure 12 :

 12 Figure 12: Lock-exchange: shape of the temperature field obtained with ISPH-USAW at t + = 10.

Figure 13 :

 13 Figure 13: Lock-exchange: temperature contours at t + = 5, 10, 15, 20. Top: ISPH-USAW. Bottom: Härtel et al. [17].

Figure 14 :

 14 Figure14: Lock-exchange: temperature contours at t + = 5, 10, 15, 20 obtained with ISPH-USAW using an antisymmetric SPH gradient[START_REF] Guimet | A linearised turbulent production in the k-model for engineering applications[END_REF] for the pressure gradient.

Figure 15 :

 15 Figure15: Infinite channel with a complex geometry: sketch of the case. The upper wedges are isothermal at the temperature T c (in blue) and the lower wedges are isothermal at the temperature T h (in red).

Figure 16 :

 16 Figure 16: Infinite channel with a complex geometry: shape of the velocity (left) and temperature (right) fields obtained with ISPH-USAW (top) and FV (bottom) after time convergence.

Figure 17 :

 17 Figure17: Infinite channel with a complex geometry. Profiles of temperature along x + = 0.83, 1.66, 2.50 (left) and z + = 0.33, 1.0, 1.66 (right) obtained with ISPH-USAW (solid lines) and FV (points) after convergence.

Figure 18 :

 18 Figure18: Infinite channel with a complex geometry. Profiles of horizontal velocity along x + = 0.83, 1.66, 2.50 (left) and vertical velocity along z + = 0.33, 1.0, 1.66 (right) obtained with ISPH-USAW (solid lines) and FV (points) after convergence.

Figure 20 :

 20 Figure 20: Turbulent plane Poiseuille flow. Profiles of turbulent kinetic energy (top) and dissipation rate (bottom) after convergence. Comparison of ISPH-USAW and FV with DNS results provided by Kasagi & Iida [18].

Figure 21 :Figure 22 :

 2122 Figure 21: Differentially heated rectangular cavity. Sketch of the case (left) and shape of the temperature field after convergence with ISPH-USAW (middle) and FV (right).

Table . 1

 . : Values of the k -model constants[START_REF] Launder | Mathematical models of turbulence[END_REF]

	κ	C µ	C 1	C 2	σ k	σ
	0.41 0.09 1.44 1.92 1.0 1.3

Table . 2

 . : Laminar Poiseuille flow: theoretical solutions for T + and p * + .

	Case (a): TT	Case (b): QT

 (37)

Acknowledgements

This work was partly funded by the French Research Agency (CIFRE agreement # 2011-0264).