Martin Ferrand
email: martin.ferrand@edf.fr

Antoine Joly
email: antoine.joly@edf.fr

Christophe Kassiotis

Damien Violeau

Agnès Leroy

François-Xavier Morel

Benedict D Rogers

Unsteady open boundaries for SPH using semi-analytical conditions and Riemann solver in 2D

Keywords: Smoothed Particle Hydrodynamics, open boundaries, unsteady flows, confined flows, free-surface flows, Riemann invariants, inlet/outlet

Due to the Lagrangian nature of SPH, treating inlet/outlet boundaries (that are intrinsically Eulerian) is a challenging issue. An extension to the Unified Semi-Analytical boundary conditions is presented to deal with unsteady open boundaries in confined and free-surface flows. The presented method uses Riemann invariants to calculate flow properties near the open boundaries, thus allowing the possibility to treat complex shapes. Furthermore, details are presented for a parallel implementation of this method, including particle creation and deletion, updating properties of vertices and segments, and additional constraints on the time step. Simple validation cases are then displayed to illustrate the performance of the proposed method as well as the ability to deal with complex problems such as generation of water waves and free outlets.

Introduction

The intrinsic computational cost of the meshless numerical method Smoothed Particle Hydrodynamics (SPH) makes simulations where only part of a larger problem domain is treated a natural extension of the technique. More generally, open boundaries are required in many fluid mechanics simulations. For this kind of simulations, inlet/outlet conditions need to be developed. Efficient inlet/outlet boundary treatment is also required when solving coupled problems, especially when strong coupling algorithms are used, since the fields near the open boundaries are shared by the coupled models, and perturbations near the shared boundary can lead to instabilities [START_REF] Kassiotis | Coupling SPH with a 1-D Boussinesq-type wave model[END_REF].

The simplest way to treat inlet/outlet in SPH is to use a buffer layer, where the values of the fields at the boundary are imposed on several layers of particles ("buffer particles") that complete the kernel support of free particles close to the open boundary [START_REF] Vacondio | [END_REF][START_REF] Bouscasse | Multipurpose interfaces for coupling sph with other solvers[END_REF][START_REF] Braun | A framework for permeable boundary conditions in sph: Inlet, outlet, periodicity[END_REF]. In the inlet case, when entering the fluid domain a buffer particle is marked as fluid particle and is then free to move. This sudden modification can generate spurious shocks. Similarly, a fluid particle leaving the fluid domain through an outlet is first marked as a buffer particle, and some of its physical quantities are suddenly prescribed, generating shocks. In the framework of Weakly Compressible SPH (WCSPH), using Riemann solvers can partially solve this issue [START_REF] Mahmood | Absorbing inlet/outlet boundary conditions for SPH 2-D turbulent freesurface flows[END_REF][START_REF] Lastiwka | [END_REF], but modelling a complex boundary where the flow can enter and/or leave the domain remains a difficult issue.

A slightly more complicated approach would be to use mirror particles at open boundaries. This work is just starting and has been succesfully applied to Dirichlet pressure boundaries by Kunz et al. [7].

The Unified Semi-Analytical strategy proposed by Ferrand et al. [8] has shown promising results to treat both Dirichlet and Neumann conditions for wall boundary treatment in SPH (see also Leroy et al. [9]). In this approach, the boundaries are discretized using boundary elements hereafter referred to as 'segments', as well as 'vertex' particles at their intersections (see Figure 1). This treatment of the boundaries will be extended here to open boundaries. For this purpose, the main novelties are the following:

• Two additional terms in the SPH continuity equation will be introduced, naturally derived from the Unified Semi-Analytical approach;

• The masses of vertex particles evolve over time according to the desired ingoing/outgoing mass flux at the connected segments;

• The vertex particles are used to create or delete fluid particles which are free to move according to the SPH momentum equation;

• The fields at the open boundaries are specified from the Riemann invariants of the Euler equations.

All these changes allow the formulation to treat particles entering and exiting the domain without perturbations of the fields for both confined and free-surface flows.

The outline of the paper is as follows: first a short description of the Unified Semi-Analytical boundary conditions for walls as proposed in Ferrand et al. [8] is given in Section 2. These wall conditions will be extended to inlets and outlets (Section 3). It will be followed by a method prescribing the fields of the ingoing/outgoing particles from Riemann invariants. These developments will the be tested in Section 4, where the ability of the present method to treat wallbounded and free-surface flows with complex open boundaries will be tested.

The work presented herein will use only the WCSPH approach (for more details about this method, see Monaghan [10]). Developing similar open boundaries in an Incompressible SPH formulation is possible [START_REF] Leroy | A new incompressible SPH model : towards industrial applications. Modeling and Simulation[END_REF], however the equations in this case are elliptic, which will change the treatment of fields at the open boundaries. This article will focus on treating the hyperbolic equations of the WCSPH approach, and the treatment to impose compatible pressure and velocity fields on the open boundaries. Treating the mass fluxes is equivalent for both WCSPH and ISPH approaches.

Finally, this approach will be restricted to 2-D one-fluid flows.

Space discretization

As illustrated in Figure 1, the weakly compressible fluid domain Ω is discretised by a set of SPH fluid particles F denoted by the subscripts (.) a and (.) b , while the boundaries ∂Ω are discretized by a set of vertex particles V denoted by (.) v and connected to boundary segments S denoted by (.) s . Let P denote the particle set that is the union of the sets F and V. Each element (particle or segment) stores information such as mass m a , position center r a , particle velocity v a (i.e. the Lagrangian derivative of the position), fluid velocity u a (i.e. the velocity field at position r a ), density ρ a , dynamic viscosity µ a ≡ νρ a , pressure p a and volume V a = m a /ρ a .

It is important to underline that two sets of velocities are needed here, since the vertex particles and boundary segments at open boundaries are fixed in space but carry an information on the fluid velocity, namely u v and u s . In other words, we have u a = v a except for vertex particles and segments located onto the open boundaries. On the other hand, this equality holds for wall segments and vertex particles, thanks to the no-slip condition. Furthermore, in our model all particle masses m a are constant in time, with the exception of the vertex particles of open boundaries, as detailed later.

Let ρ 0 be the reference density of the considered fluid. In the Weakly Compressible SPH approach used herein, the pressure field is deduced from the density field using Tait's equation of state [START_REF] Tait | Report on the Voyage of Challenger[END_REF], given as:

p a = ρ 0 c 2 0 ξ ρ a ρ 0 ξ -1 , (1)
where c 0 is the numerical speed of sound and ξ is taken to 7 for water.

The SPH interpolation is based on a weighted interpolation from a kernel function w, generally compactly supported. Here we use the Wendland kernel of order 5 [13]. In our notation, Ω a refers to the support of the kernel function centred on r a . The subscripts (.) ab generally denotes the difference of a quantity between the particles a and b (unless stated otherwise). For instance u ab ≡ u au b and r ab ≡ r ar b . Some exceptions are made, including the following notations: w ab ≡ w (r ab ) and ∇w ab ≡ ∇ a w (r ab ). The symbol ∇ a denotes the gradient with respect to the position r a .

Boundary renormalization

Following the work of several authors, in particular Kulasegaram et al. [14] and Ferrand et al. [8], an additional field denoted by γ a is used to include boundary terms. This field measures the part of the kernel support which is inside the computational domain and is defined as:

γ a ≡ Ωa∩Ω w (r a -r) dr. (2)
The field γ a is computed from a dynamic governing equation (see Section 2.5). The discrete SPH operators presented in the next section also require the computation of ∇γ a , which is performed by a decomposition onto the boundary segments. Each segment contribution ∇γ as is defined as:

∇γ as ≡ s w (r) dS n s , (3)
where n s is the inward unit normal to the boundary segment s. Now ∇γ a is written as:

∇γ a = s ∈ S ∇γ as . (4)
This gradient can be calculated analytically, see Ferrand et al. [8].

Space differential corrected operators

In the present section, space discretized differential operators defined by Ferrand et al. [8] will be used. These include G γ, + a {A b } which is a boundarycorrected gradient of the discrete scalar field {A b }, D γ,- a {A b } which is a boundarycorrected divergence of the discrete vector field {A b }, and L γ a {B b , A b } which is a boundary-corrected Laplacian of the discrete scalar (or vector) field {A b } with discrete diffusion field {B b }. These operators are defined by the following equations:

G γ,+ a {A b } ≡ ρ a γ a b ∈ P m b A a ρ 2 a + A b ρ 2 b ∇w ab - ρ a γ a s ∈ S A a ρ 2 a + A s ρ 2 s ρ s ∇γ as , (5)
D γ,- a {u b } ≡ - 1 γ a ρ a b∈P m b u ab • ∇w ab + 1 γ a s∈S u as • ∇γ as , (6)
L γ a {B b , A b } ≡ 1 γ a b ∈ P V b 2B ab A ab r 2 ab r ab • ∇w ab - 1 γ a s ∈ S (B s ∇A s + B a ∇A a ) • ∇γ as . (7)
The factor B ab is defined using a combination of B a and B b . In most of the SPH literature, it is taken as their arithmetic mean. However, to ensure continuity of shear flux, i.e. B a ∇A a • r ab = B b ∇A b • r ab , here it is their harmonic mean:

B ab = 2B a B b B a + B b . (8)
Furthermore, as a first order approximation the two boundary terms of Equation 7 can be equalised, so that

(B s ∇A s + B a ∇A a ) • ∇γ as ≈ 2 |∇γ as | B s ∇A s • n s .
This is the flux at the boundary, therefore if a Neumann boundary condition is imposed, then it will replace this term and if a Dirichlet condition is imposed then it can be replaced by -2 |∇γ as | B as A as /(r as • n s ), as explained by Ferrand et al. [8] and Leroy et al. [9].

Space discretized equations

The equations to be solved are the momentum equation, the equation of motion, the governing equation for γ a and the continuity equation:

ρ a du a dt = -G γ,+ a {p b } + L γ a {µ b , u b } + ρ a g, dr a dt = v a , dγ a dt = ∇γ as • v as , dρ a dt = -ρ a D γ,- a {u b }. (9)
It is important to underline that the particle velocity v a is used to move the particles and update γ a , while the fluid velocity u a is used in viscous forces and density evolution.

When time is considered as continuous and in the absence of open-boundaries, the last two lines of Equations ( 9) are fully equivalent to [8]:

d (γ a ρ a ) = d b∈P m b w ab , (10)
This formulation allows exact integration with time, and thus will be preferred in the following. However, in presence of inlets or outlets boundaries, the equivalence between Equations (9) and Equation (10) needs to be modified. The necessary developments are detailed in Section 3.1.

Time stepping without open boundaries

The time-stepping scheme used in Ferrand et al. [8] to solve the system of equations ( 9) is recalled hereafter:

f n a = L γ a {µ n b , u n b } -G γ,+ a {p n b } + ρ n a g, u n+1 a = u n a + ∆t ρ n a f n a , v n+1 a = u n+1 a , r n+1 a = r n a + ∆tv n+1 a , γ n+1 a = γ n a + ∆t s∈S 1 2 ∇γ n+1 as + ∇γ n as • v n+1 as , ρ n+1 a = 1 γ n+1 a γ n a ρ n a + b∈P m b w n+1 ab -w n ab , (11)
where p n = p(ρ n ) using Equation (1). Ferrand et al. [8] considered v a = u a since they did not deal with open boundaries, contrary to the present work.

In Section 3, this scheme will be adapted to account for ingoing or outgoing mass fluxes at open boundaries.

Volume diffusion correction

In the SPH operators used in Equations 9, the pressure and velocities are collocated, i.e. stored at the same points (the particle positions), and therefore require stabilization. Several authors have developed stabilization terms for SPH in order to solve this issue. We will refer the reader to the work by Fatehi and Manzari [15] or to the work by Ferrari et al. [16] whose work was adapted to Semi-Analytical SPH boundaries by Mayrhofer et al. [START_REF] Mayrhofer | Study of differential operators in the context of the semi-analytical wall boundary conditions[END_REF].

In the present work the stabilization factor is derived from the pressurevelocity coupling available in the literature for finite elements mesh-based method, namely Brezzi and Pitkäranta [START_REF] Brezzi | On the stabilization of finite element approximations of the stokes equations[END_REF], which was adapted to Semi-Analytical SPH boundaries by Ghaitanellis et al. [START_REF] Ghaitanellis | Application of the unified semi-analytical wall boundary conditions to multi-phase sph[END_REF]. The main principle behind this correction is to develop the term in u n+1 b in the continuity equation. In the absence of viscosity, this gives:

ρ n+1 a -ρ n a ∆t = -ρ n a D γ,- a u n b - ∆t ρ n a G γ,+ a {p n b } + ∆t ρ n a ρ n a g ≈ -ρ n a D γ,- a {u n b } -ρ n a D γ,- a - ∆t ρ n a G γ,+ a {p n b } + ∆tG γ,+ a {g • r n b } (12)
From this reasoning, a diffusion term ∆ can be defined as:

∆ n = -ρ n a -L γ a ∆t ρ n b , p n b + L γ a {∆t, g • r n b } (13)
Note that ∆ is not exactly equal to the second term of Equation ( 12) since with SPH operators the Laplacian is not tantamount to the divergence of a gradient. Following Brezzi and Pitkäranta [START_REF] Brezzi | On the stabilization of finite element approximations of the stokes equations[END_REF] this diffusion term will be weighted by an arbitrary factor Λ ∈ [0; 1], and it will be added to the continuity equation. The last of Equations [START_REF] Leroy | A new incompressible SPH model : towards industrial applications. Modeling and Simulation[END_REF] will therefore be rewritten as:

ρa = 1 γ n+1 a γ n a ρ n a + b∈F ∪V m b w n+1 ab -w n ab , ρ n+1 a -ρa ∆t = Λρ a L γ a ∆t ρb , pb -L γ a {∆t, g • r b } , (14)
where pb = p(ρ b ), and Λ is usually set to 0.1. In addition, the Laplacian operator in Equation ( 14) excludes the boundary elements. This correction is similar to the one developed by Mayrhofer et al. [START_REF] Mayrhofer | Study of differential operators in the context of the semi-analytical wall boundary conditions[END_REF], but it is still valid when the density presents large variations.

Boundary conditions

Discretising the boundary in terms of segments and vertices gives a method of dealing with the boundaries on the domain, however the method for imposing boundary conditions is not trivial. The approach chosen by Ferrand et al. [8], and later extended in Leroy et al. [START_REF] Leroy | [END_REF], is that Dirichlet boundary conditions will be imposed on the vertices and Neumann boundary conditions will be imposed on the segments. An averaging is then used to find the corresponding values on the segments or vertices.

For example, the shear stress on the walls (calculated either through the no-slip condition or through a turbulence wall function) is set by imposing the fluid velocity on the wall vertices, and the shear stress of the segments is then found by averaging the velocity of the connected vertices.

Formulation for unified semi-analytical open boundary condition

Time integration of the continuity equation with Unified Semi-Analytical

open-boundary conditions In this section, the continuity equation will be carefully integrated between two successive iterations. The aim is to integrate exactly all the terms which depend only on particle positions, in order to stop the imposed mass flux from introducing errors on the density at each time step. Otherwise, these errors would lead to density discontinuities for fluid particles near open boundaries, which would in turn result in spurious shockwaves when particles are created or deleted.

While integrating the continuity equation, new boundary terms will naturally appear due to the presence of open boundaries. Their time integration will be performed consistently with the time stepping scheme proposed in Section 2.5, so that no artificial density perturbation is induced by the open boundaries.

As stated above, the time-dependent SPH continuity equation is the fourth line of the system (9) involving the fluid velocity u b . Let us now rewrite this equation as a total derivative so that the time integration will be exact. To do so, the continuity equation is first rewritten using the divergence-operator definition of Equation ( 6):

dρ a = -ρ a D γ,- a {u b } dt = 1 γ a b∈P m b ∇w ab • u ab dt - ρ a γ a s∈S ∇γ as • u as dt. (15)
To transform the sums in Equation ( 15) into exact total derivatives, the fluid velocity u b should be replaced by the particle velocity v b . This is possible for all particles except for the sets of vertices and segments which belong to the open boundaries (denoted by V i/o and S i/o , respectively), and thus the density variation reads:

dρ a = 1 γ a b∈P m b ∇w ab • v ab dt - 1 γ a v∈V i/o m v ∇w av • (u v -v v ) dt - ρ a γ a s∈S ∇γ as • v as dt + ρ a γ a s∈S i/o ∇γ as • (u s -v s ) dt. (16)
In order to compute the first and third terms of Equation ( 16), we use the fact that in a Lagrangian frame dr = vdt, so that:

dw ab = ∇w ab • v ab dt, dγ a = s∈S ∇γ as • v as dt. ( 17
)
Therefore Equation ( 16) can be rewritten as:

d (γ a ρ a ) = γ a dρ a + ρ a dγ a = d b∈P m b w ab -γ a δρ i/o a + ρ a δγ i/o a , (18)
where

δρ i/o a and δγ i/o
a are defined by:

δρ i/o a ≡ 1 γ a v∈V i/o m v ∇w av • (u v -v v ) dt, ( 19
)
δγ i/o a ≡ s∈S i/o ∇γ as • (u s -v s ) dt. ( 20
)
It should be noted that Equation ( 18) is a generalized version of Equation (10). Furthermore the terms in Equations ( 19) and ( 20) are non-zero because particle and fluid velocities do not coincide at open boundaries, as already mentioned. In addition, for Equation [START_REF] Brezzi | On the stabilization of finite element approximations of the stokes equations[END_REF] to be valid, the particle masses should be constant over the time integration of γ a ρ a . Therefore, the particles masses m b are kept constant on the time interval t n ; t n+1 equal to m n b . Since the vertex masses will vary in time (see Section 3.3), this implies that these variations will occur after the time integration of the continuity equation.

Let us now integrate the continuity Equation ( 18) from time t n to t n+1 = t n + ∆t, to have a consistent time marching scheme:

(γ a ρ a ) n+1 -(γ a ρ a ) n = b∈P n m n b w n+1 ab -m n b w n ab - t n+1 t n γ a δρ i/o a + t n+1 t n ρ a δγ i/o a . (21)
Note that the summation is over the particles existing at time n, denoted by P n and the operator, dt, in the integrals is included in Equations ( 19) and [START_REF] Leroy | [END_REF]. The particle creation/deletion algorithm is described in Section 3.3.

The two time-integrated terms

t n+1 t n γ a δρ i/o a and t n+1 t n ρ a δγ i/o
a in Equation ( 21) will be called virtual variation terms, and are the only modification of the time-independent continuity equation proposed by Ferrand et al. [8] (last line of Equation 11).

The key point in the time integration of the virtual variations is that they are expected to maintain exactly a uniform flow field with a constant density. The term γ a δρ i/o a represents the variation experienced by v m v w av , with v ∈ V i/o , if the vertex particles v were moved with the velocity (u vv v ) while the fluid particle a is fixed. Similarly the term δγ i/o a represents the variation experienced by γ a if the segments s ∈ S i/o were moved with the velocity (u sv s ) while the fluid particle a is fixed.

Consequently the virtual displacement

δr i/o a ≡ ∆t (u n a -v n a
) is used to compute the virtual variation terms. Thus the time-integration of δρ i/o a and δγ i/o a can be computed consistently with the last two lines of Equation ( 11) as:

t n+1 t n γ a δρ i/o a = v∈V i/o m n v w r n av + δr i/o v -w (r n av ) , (22)
and:

t n+1 t n ρ a δγ i/o a = ρ n a 2 s∈S i/o ∇γ as r n as + δr i/o s + ∇γ as (r n as ) • δr i/o s , (23)
where

δr i/o s = ∆t(u n s -v n s )
, the velocity of the segments are defined by u s = (u v1 + u v2 )/2 and δr

i/o v = ∆t(u n v -v n v ).
In the above notation, v 1 and v 2 are the vertices connected to s.

Equation ( 23) is consistent with the way γ is calculated. If an analytical formulation of γ is used (such as in Leroy et al. [9]) it should be computed as:

t n+1 t n ρ a δγ i/o a = ρ n a s∈S i/o γ as r n as + δr i/o s -γ as (r n as ) . (24)
The proposed continuity Equation ( 21) allows to take account of ingoing/outgoing mass flux through the last two terms, while other terms depend only on particle positions; Equation (21) combined with Equations ( 22) and (23) gives:

ρ n+1 a = 1 γ n+1 a γ n a ρ n a + b∈P n m n b w n+1 ab -w n ab + v∈V i/o m n v w n av -w r n av + δr i/o v + ρ n a 2 s∈S i/o ∇γ as r n as + δr i/o s + ∇γ as (r n as ) • δr i/o s    . ( 25
)
For simplicity the sumations in the right hand side of Equation ( 25) will be noted as d i/o a , which gives:

ρ n+1 a = 1 γ n+1 a γ n a ρ n a + d i/o a . (26)

Time stepping with open boundaries

The time stepping scheme used in Equation ( 11) is therefore adapted to take into account the open boundaries with a volume diffusion correction (section 2.6):

f n a = L γ a {µ n b , u n b } -G γ,+ a {p n b } + ρ n a g, u n+1 a = u n a + ∆t ρ n a f n a , v n+1 a = u n+1 a for a ∈ F/V, r n+1 a = r n a + ∆tv n+1 a , γ n+1 a = γ n a + ∆t s∈S 1 2 ∇γ n+1 as + ∇γ n as • v n+1 as , (27)
ρa = 1 γ n+1 a γ n a ρ n a + d i/o a , ρ n+1 a -ρa ∆t = Λρ a L γ a ∆t ρb , pb -L γ a {∆t, g • r b } . As a reminder d i/o
a is defined through Equations ( 25) and (26). Furthermore, the constraints on the time step ∆t are described in Section Appendix B.4.

Mass update of vertex particles

The method used to create or remove particle needs to be carefully developed in order to keep a correct particle distribution near the open boundaries (this process can be understood as re-meshing). The choice made is to make the vertex particles grow (respectively decrease) over time on an inlet (respectively outlet) boundary. Thus their masses and volumes are not constant. However, the mass evolution should not create any perturbation neither on the density nor on the momentum. Therefore, as stated earlier a constant mass will be set between t n and t n+1 in Equations [START_REF] Jeffrey | Quasilinear Hyperbolic Systems and Waves[END_REF]. This is ensured by Equations ( 21), [START_REF] Archambeau | Industrial Applications[END_REF] and (23). The masses will therefore be updated outside of the marching scheme.

Nonetheless, the total mass of the system should change smoothly in time. To do so the masses of the vertex particles in v ∈ V i/o take into account the desired mass flux through the open boundary ṁv . This mass flux is computed as a weighted average of the mass fluxes through the segments s ∈ N sv , where N sv is the set of all the segments connected to v (usually two segments in 2-D):

ṁv = 1 2 s∈N i/o sv ρ s S s (u s -v s ) • n s , ∀v ∈ V i/o , (28)
where S s is the length of a segment. The flux ṁv is positive for an inlet and negative for an outlet.

. .

⃗ u s . s . v (a) t0.
.

.

. The flux ṁv is used to update the vertex masses. Instinctively one would assume that m n+1 v = m n v + ∆t ṁn v . However, in the case of an inlet, new fluid particles must be created near the open boundary when the fluid particles move away from it (see Figure 2). For this purpose, each time the vertex particle mass m v reaches a certain threshold θ v m ref , a fluid particle is created at the same position with the reference mass m ref and this mass is subtracted from the vertex particle mass to prevent from any density shock. For the sake of simplicity θ v is the open angle of the vertex divided by 2π (i.e. θ v = 0.5 for a plane boundary) and m ref is the mass of a fluid particle.

In the case of an outlet, particles are removed when they cross an open boundary segment. In order to remain consistent the mass of an outgoing fluid particle is redistributed over the open boundary vertices connected to this segment. Each of these vertices receives the mass β a,v m a , where β a,v is a distribution (or weighing) factor.

The factor β a,v is constructed so that most of the mass is distributed to the closest vertex. For a segment s with vertices v 0 and v 1 , vector p i is defined as the vector between vertex v i and the projection of the particle a on the segment, i.e. p i = r avi -(r avi • n) n (see Figure 3). Assuming that the projection belongs to the segment, then for v 0 and v 1 connected to s, β a,vi is defined as:

r v 0 r v 1 r a p 0 n p 1
β a,v0 = p 1 • r v0v1 r v0v1 2 , β a,v1 = p 0 • r v1v0 r v1v0 2 = 1 -β a,v0 . (29)
Otherwise, β a,v is set to zero. Therefore, the masses of vertex particles are updated as follows:

m n+1 v = m n v + ∆t ṁn v + δm n v , ∀v ∈ V i/o , (30)
where δm n v quantifies the mass variation due to particle creation/destruction, i.e. -m ref each time a fluid particle is created at position v, +β a,v m a each time a particle a crosses a segment connected to v. As stated above, the step described by Equation ( 30) is added at the end of the time stepping scheme [START_REF] Jeffrey | Quasilinear Hyperbolic Systems and Waves[END_REF].

In Equation [START_REF] Hérault | Achieving the best accuracy in an sph implementation[END_REF] only the mass flux positiveness determines whether an open boundary repulses or attracts fluid particles. As illustrated later in the periodic wave test case, i.e. Figure 12, it is possible to prescribe time-dependent velocities moving from negative to positive values. Henceforth, the proposed strategy makes it straightforward to treat both inlet and outlet at the same open boundary, depending on space and time.

Free-surface flows can also be treated that way. For example, for the periodic wave shown in Figure 12 and later in section 4.3, the mass flux ṁv at the vertical open boundaries oscillates consistently from positive values to negative values as the velocity is prescribed according to fifth-order solution to Stokes wave theory.

Imposing pressure and velocity

The Riemann problem formulation described in Appendix A is used to define compatible pressure and velocity fields. The principle is to used a method common in finite volumes, where characteristic waves (λ -1 , λ 0 and λ +1 ) model discontinuities between the exterior state (the boundary conditions) and the interior state (the fluid) through Generalized Riemann Invariants (GRI).

Imposed velocity

When a velocity field is imposed on a boundary, for instance an incoming free-surface wave velocity profile at an inlet, it is imposed on the open boundary vertices (as it is a Dirichlet boundary condition). Since the imposed velocity field is needed on the open boundary segments, its value will be averaged from the neighbourghing vertices:

u s = 1 |N vs | v∈Nvs u v . (31)
This value will then be the prescribed velocity field on the exterior state, i.e. u ext = u s . In this case the linearised Riemann problem is used to compute the value of the density (or pressure) fields at the exterior state, and therefore at the boundary segment and vertices.

The GRI of wave λ +1 will therefore be used to calculate the external state pressure p ext from the internal state. However the type of discontinuity of the wave λ +1 needs to be defined. It is known that since λ +1,ext = λ +1,2 , then if λ +1,ext > λ +1,int the wave is a shock, however since p ext is unknown then c ext cannot be calculated. Therefore assumptions needs to be made regarding the relationship λ +1,ext > λ +1,int . It will be assumed that this condition is verified if u n,ext > u n,int , but this assumption will need to be verified once c ext is known. Therefore if:

1. If u n,ext > u n,int , the λ +1 wave is a shock wave. The Rankine-Hugoniot relationships have to be used through the following equation:

p ext = p int + ρ int u n,int (u n,int -u n,ext ) . ( 32
)
2. If u n,ext ≤ u n,int , the λ +1 wave is an expansion wave. The corresponding Riemann invariant R + must must be used (i.e. R +1,ext = R +1,int ):

u n,ext -ψ (ρ ext ) = u n,int -ψ (ρ int ) , (33)
where ψ is defined by Equation (A.8).

Once the pressures are known on the segments (since p s = p ext ), then the pressure on the vertices are found by averaging the connected segments:

p v = 1 |N sv | s∈Nsv p s . (34)
The densities are once again calculated using Equation (1). As mentioned earlier, if the imposed condition corresponds to an inlet then the tangential velocity components need to be defined, otherwise they will be defined from the interior state.

Imposed pressure (or density)

When a pressure field is imposed on a boundary, for instance hydrostatic pressure at an outlet, it is imposed on the vertices (as it is a Dirichlet boundary condition). Since the imposed pressure field is needed on the open boundary segments, to do so its value will be averaged from the neighbourghing vertices:

p s = 1 |N vs | v∈Nvs p v . (35)
This value will then be the prescribed velocity field on the exterior state, i.e. p ext = p s . In this case the linearised Riemann problem is used to compute the value of the velocity fields at the exterior state, and therefore at the boundary segment and vertices. The GRI of wave λ +1 will therefore be used to calculate the normal velocity component from the internal state. However the type of discontinuity of the wave λ +1 also needs to be defined. In this case it will be assumed that this condition is verified if c ext > c int , but this assumption will need to be verified once u n,ex t is known. Therefore if:

1. If c ext > c int , the λ +1 wave is a shock wave. The Rankine-Hugoniot relationships have to be used through the following equation:

u n,ext = p int -p ext ρ int u n,int + u n,int . (36)
Note that due to the presence of u n,int in the denominator of Equation (36), when the pressure is imposed on a boundary adjacent to still water this formulation could pose difficulties. This is why, when u n,int → 0 the limit of Equation (36) should be considered, i.e.:

lim un,int→0 u n,ext = u n,int . (37)
Practically, this condition will be imposed when u n,ext > c 0 /10. 2. If c ext ≤ c int , the λ +1 wave is an expansion wave. The corresponding Riemann invariant R + must be used (i.e. R +1,ext = R +1,int ):

u n,ext -ψ (ρ ext ) = u n,int -ψ (ρ int ) . ( 38
)
Once the velocities are known on the segments (since u s = u ext ), then the velocities on the vertices are found by averaging the connected segments:

u v = 1 |N sv | s∈Nsv u s . (39)
Again, if the imposed condition corresponds to an inlet then the tangential velocity components need to be defined, otherwise they will be defined from the interior state. Finally, the densities are still calculated using Equation (1).

Imposed open boundary conditions algorithm

The linearised Riemann problem is used to calculate the open boundary conditions before the time stepping is done. Therefore for an open boundary problem, the calculation steps are listed below:

1. Calculation of the fields at the boundaries, using the method described in Ferrand et al. [8] for wall boundaries, and the method described in Section 3. The parallel implementation for this algorithm is described in Appendix B.

Numerical tests

For all the simulations, the quintic Wendland kernel with smoothing length h = 2∆r will be used.

Non-orthogonal flux on inlet/outlet in a square

The first test case presented is that of a 2-D square fluid domain of size L with only inlet/outlet boundary conditions (see Figure 4). The flow is at an angle of π/4 to the normal to the boundaries; in other words, the velocity imposed at all boundaries is u 0 = |u 0 |/ √ 2(1, 1) T . It illustrates one of the advantages of present Semi-Analytical inlet/outlet boundary treatment: the absence of buffer layer makes it possible to prescribe a non-orthogonal flow quite easily, regardless of the open boundary shape.

The initial conditions are of constant velocity and density fields with the same values as imposed at the boundary. Particles are initially distributed onto a Cartesian grid, the inter-particle distance ∆r and segment size are nondimensionalized using the problem size, i.e. ∆r + = ∆r/L for all different tests performed (from 0.8k to 41.6k particles). Furthermore the speed of sound at rest is set to c 0 = 10|u 0 |. The fluid viscosity ν is set to 10 -2 m 2 .s -1 and the volumic diffusion factor Λ will be set to 0, so that the focus of this case is on the open boundary conditions.

The evolution of the L 2 errors for both density and velocity is computed as:

ε ρ = 1 N P a ∈ P (ρ a -ρ 0 ) 2 ρ 2 0 and ε u = 1 N P a ∈ P |u a -u 0 | 2 |u 0 | 2 , ( 40
)
where N P is the number of fluid and vertex particles. As a first test the Riemann invariants will not be used and both the velocity and density will be imposed as u 0 and ρ 0 respectively. The evolution of the errors is displayed in Figure 5 with respect to non-dimensional time t + = t|u 0 |/L. It is therefore shown that the errors are very small; they are mainly linked to . the round-off errors. Since the simulations presented are computed on single precision GPU (Graphics Processing Units), the values observed are acceptable.

It should be emphasised that the additional terms in the continuity equation (Equation 18) are essential, as the authors have observed that the simulations that ran without those terms crashed immediately. As a second test, the velocity are imposed on inlets and the pressure on the outlets together with Riemann invariants. As shown by Figure 6, the errors are slighlty larger than previously. This is because the previous case was overconstrained. Futhermore, when using Riemann invariants at the open boundaries the error is almost independent of the particle size. However the error on the velocity is of an order of magnitude higher than that of the density. This can easily be explained by linearizing the Euler equations for the numerical perturbations δρ et δu:

δu τ ∼ c 2 0 ρ 0 δρ λ , ( 41
)
δρ τ ∼ ρ 0 δu λ , (42)
where τ and λ are time and length scales corresponding to the perturbation propagation.

The order of magnitude of the errors are ε ρ ∼ δρ/ρ 0 and ε u ∼ δu/|u 0 |. Combining the last two equations in order to remove τ and λ yields ε u ∼ (c 0 /|u 0 |) ε ρ . The error on the velocity is thus amplified by the numerical Mach number, which is taken as 10 as explained above.

Rapidly expanding pipe

In this next test case, the flow from a small pipe (named pipe 1) will enter a large pipe (named pipe 2), creating a sudden expansion of the flow. See Figure 7 for details on the geometry. Far from the expansion, the flow should follow the Hagen-Poiseuille equations, which for a flow between two infinite plates is given by:

Imposed Velocity Imposed Pressure 2W 1 2W 2 L 1 L 2 x z
u x = U max 1 - z 2 W 2 , ( 43
)
dP dx = - 2ρ 0 νU max W 2 , ( 44
)
where the origin of z lies on the central pipe axis, W is the half width of the pipe and U max is the maximum velocity. To ensure that the flow stays laminar, this maximum velocity will be defined according to a Reynolds number Re:

U max = νRe 2W . (45)
Equation ( 43) will therefore be used to impose a Poiseuille flow velocity profile on the inlet, and a constant pressure will be imposed on the outlet. The geometrical parameters are chosen to correspond to one of the experiments presented by Hammad et al. [21]: Figure 8 shows the horizontal velocity and density profiles at the inlet and outlet. As can be seen the velocity profile at the outlet is very close to the Poiseuille flow. In addition, the density at the inlet is greater than at the outlet, which is a necessary condition as the pressure is linked to the density through the state Equation (1). Finally, integrating the profile in Figure 8 • Mass flux in: 0.523 kg.m -1 .s -1

W 1 = 0.13 m, L 1 = 3W 1 , W 2 = 2W
• Mass flux out: 0.525 kg.m -1 .s -1

• Analytical mass flux: 0.521 kg.m -1 .s -1

The mass flux differences appear because even though the fluid particles have fixed masses throughout the simulation, the open boundary vertices have a variable mass. The mass variations are consistent with the boundary mass flow, and the differences between the inlet and the outlet are a result of the compressible medium which cannot achieve a perfect steady state. The normalised error in the mass flux is therefore equal to 4 × 10 -3 , which is acceptable relative to the spatial discretisation, and it does not have a significant impact on the simulation. In Figure 9 the streamlines obtained from the stream function are compared with the experimental results presented in Hammad et al. [21] and a finite volume solution using Code_Saturne [START_REF] Archambeau | Industrial Applications[END_REF]. The calculation of the stream function Ψ from a velocity field can be performed using different methods and can provide insight into sources of disagreement with the experimental values. To analyse the simulation results the following equation was used:

1 1 1 1 2 2 2 2 -2 -1 0 1 2 3 4 -1 0 1 x/W 2 z/W 2 Poiseuille flow

Poiseuille flow Recirculation

Ψ(x, z) =          r(x,z) r(L2,W2) (u x dz -u z dx)
above the central axis

r(x,z) r(L2,-W2) (u x dz -u z dx) below the central axis (46)
Amongst the differences between the SPH simulation and the experimental results, the stream line at Ψ = 0 does not cross the wall at the same distance from the expansion, meaning that the recirculation zone is smaller in the simulations. This could be because in the state Equation (1) the value of ξ for diethylene glycol (which was used in the experiment) is not readily available, and therefore ξ = 7 was imposed, which is the value for water. This is also true for the simulations using Code_Saturne, although the incompressible model used in this case allows these results to be closer to the experimental values.

Nonetheless, at the outlet the streamlines of the simulations are close to the experimental results. As a further check of the validity of the simulation, the head loss will be modelled. Looking at the streamlines of the flow, plotted in the flow can be separated into 3 regions: a Poiseuille flow inside pipe 1, a recirculation zone created by the sudden expansion and a Poiseuille flow inside pipe 2 (see Figure 10). One can estimate the head loss in each section by (see Hammad et al. [23] for more details):

             1 2 U 2 1 -U 2 1 = 1 ρ 2 (P 1 -P 1 ) , 1 2 U 2 1 -U 2 2 + 1 ρ 2 (P 1 -P 2 ) = 1 2 U 1 -U 2 2 , 1 2 U 2 2 -U 2 2 = 1 ρ 2 (P 2 -P 2 ) , (47)
where the subscript represent the interface between different zones of the flow. Subscript 1 denotes the interface between the inlet and pipe one, subscript 1 the interface between the Poiseuille flow of pipe one and the recirculation area, subscript 2 the interface between the recirculation area and the Poiseuille flow of pipe two and subscript 2 the interface between pipe two and the outlet (see Figure 10). Furthermore

U 2 = 1/(2W ) W -W U 2 d z and P = 1/(2W )
W -W pd z. These three equations can be combined to give:

1 2 U 2 1 -U 2 2 + 1 ρ 2 (P 1 -P 2 ) = 1 2 U 1 -U 2 2 + 1 ρ 2 (P 1 -P 1 ) + 1 ρ 2 (P 2 -P 2 ) + O ν U 1 W 2 . ( 48
)
Defining the head of the flow by:

H = U 2 2g + P gρ 2 . ( 49
)
The theoretical head loss between the inlet and the outlet is thus equal to (neglecting the last term in Equation 48): The expected head loss is therefore equal to the sum of head loss of the Poiseuille flows, and Borda-Carnot head loss in the recirculation zone. Figure 11 shows the head along the x-axis. Extracting the simulated head loss and calculating the theoretical head loss gives:

H 1 -H 2 = 1 2g U 1 -U 2 2 + 1 gρ 2 (P 1 -P 1 ) + 1 gρ 2 (P 2 -P 2 ) (50)
-ν 2 R e / ( g H 1 W 3 1 ) -ν 2 Re /(g H 1 W 3 2 ) 1 2gH 1 U 1 -U 2 2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 3.5 4 0.4 0.6 0
• Head loss from simulations = 5.04 × 10 -7 m

• Theoretical head loss = 5.00 × 10 -7 m Therefore, the simulated head loss has less than a 1% difference with the theoretical head loss.

2-D periodic free-surface water wave

The next test case shows that these new open boundary conditions can be used for both inflows and outflows, and that open boundaries can alternate between the two without any difficulties, even with a free-surface. To do so periodic water waves propagating on a flat bed will be imposed. The waves imposed are calculated from the fifth-order solution to Stokes wave theory given by Fenton [24]:

kη(x, t) = 5 i=1 i i j=i B ij cos jk x -ct + θ k (51) U k g 1/2 = C 0 + 2 C 2 + 4 C 4 (52) c = U c + U (53) u x (x, z, t) = c -U + C 0 g k 3 1/2 5 i=1 i i j=1
A ij cosh (jkz) jk cos jk xct + θ k (54)

u z (x, z, t) = C 0 g k 3 1/2 5 i=1 i i j=1 A ij jk sinh (jkz) sin jk x -ct + θ k (55) Rk g = 1 2 C 2 0 + kD + 2 E 2 + 4 E 4 (56) p(x,z,t) ρ0 = R -gz -1 2 (u x -c) 2 + u 2 z ( 57
)
where η is the free-surface elevation (with the water depth h = D + η), A is the wave amplitude (i.e. half of the wave height), D is the mean water depth, g is the acceleration due to gravity, k is the wave number (defined as k = 2π/L, where L is the wave length), c is the wave velocity, U c is the mean current velocity, U is the mean horizontal velocity, θ is a phase constant and R is the Bernouilli constant. The wave period T can be calculated from the wave number and wave velocity, i.e. T = 2π/(kc). Furthermore the z-axis has its origin at the bed. The constants are given in table 1.

For the present test case a depth of 0.5 m, an amplitude of 0.05 m and a wave length of 2.5 m have been chosen. Furthermore there will be no mean Table 1: Constants for the fifth-order solution to Stokes wave theory given by Fenton [24] and used in equations ( 51) -(57).

A11 = 1/ sinh (kD) A22 = 3S 2 /[2(1 -S) 2 ] A31 = (-4 -20S + 10S 2 -13S 3 )/[8 sinh(kD)(1 -s) 3 ] A33 = (-2S 2 + 11S 3 )/[8 sinh(kD)(1 -s) 3 ] A42 = (12S -14S 2 -264S 3 -45S 4 -13S 5 )/[24(1 -S) 5 ] A44 = (10S 3 -174S 4 + 291S 5 + 278S 6 )/[48(3 + 2S)(1 -S) 5 ] A51 = (-1184 + 32S + 13232S 2 + 21712S 3 + 20940S 4 + 12554S 5 -500S 6 -3341S 7 -670S 8 )/[64 sinh(kD)(3 + 2S)(4 + S)(1 -S) 6 ] A53 = (4S + 105S 2 + 198S 3 -1376S 4 -1302S 5 -117S 6 + 58S 7 ) /[32 sinh(kD)(3 + 2S)(1 -S) 6 ] A55 = (-6S 3 + 272S 4 -1552S 5 + 852S 6 + 2029S 7 + 430S 8 ) /[64 sinh(kD)(3 + 2S)(4 + S)(1 -S) 6 ] B11 = 1 B22 = coth(kd)(1 + 2S)/[2(1 -S)] B31 = -3(1 + 3S + 3S 2 + 2S 3 )/[8(1 -S) 3 ] B33 = -B31 B42 = coth(kd)(6 -26S -182S 2 -204S 3 -25S 4 + 26S 5 /[6(3 + 2S)(1 -S) 4 ]) B44 = coth(kd)(24 + 92S + 122S 2 + 66S 3 + 67S 4 + 34S 5 )/[24(3 + 2S)(1 -S) 4 ] B51 = -(B53 + B55) B53 = 9(132 + 17S -2216S 2 -5897S 3 -6292S 4 -2687S 5 + 194S 6 + 467S 7 +82S 8 )/[128(3 + 2S)(4 + S)(1 -S) 6 ] B55 = 5(300 + 1579S + 3176S 2 + 2949S 3 + 1188S 4 + 675S 5 + 1326S 6 + 827S 7 +130S 8 )/[384(3 + 2S)(4 + S)(1 -S) 6 ] C0 = [tanh(kD)] 1/2 C2 = [tanh(kD)] 1/2 (2 + 7S 2 )/[4(1 -S) 2 ] C4 = [tanh(kD)] 1/2 (4 + 32S -116S 2 -400S 3 -71S 4 + 146S 5 )/[32(1 -S) 5 ] E2 = tanh(kD)(2 + 2S + 5S 2 )/[4(1 -S) 2 ] E4 = tanh(kD)(8 + 12S -152S 2 -308S 3 -42S 4 + 77S 5 )/[32(1 -S) 5 ]
Where S = sech (2kD).

current in the flow. The Ursell number for this case is 2AL 2 /d 3 = 5 32π 2 /3, which would mean that linear wave theory is applicable. However, a higher order solution is used, as phase difference can occur after several time periods. Finally, in the simulation the fluid viscosity ν is set to 10 -6 m 2 s -1 , the particle spacing ∆r is set to one tenth of the amplitude and the numerical parameters are c 0 = 20 m.s -1 and Λ = 0.1. An illustration of this test case can be found in Figure 12 where it is shown that the open boundaries alternate between inlets and outlets. A simulation will be run where the fluid velocities calculated from the analytical solution (Equations 51 -57) will be imposed along the normal of the open boundaries, and Riemann invariants will be used to calculate the pressure. This will be compared to the analytical solution and a simulation with periodic conditions. Periodic conditions mean that particles exiting the domain on one side will enter the domain on the other side, but more importantly particles on one side will have as neighbours particles close to the other side, acting as if the domain was truly infinite. This will therefore serve as a reference, because ideal open boundary conditions would give the same values as an extension of the domain filled with fluid and wall particles. Furthermore, it should be noted that if both the pressure and the velocities would be imposed, then the problem would be overconstrained and particles would exit the domain.

Plots of the free-surface are then presented in Figure 13. The first conclusion is that even after several wave periods the free-surface stays very close to the analytical solution for all type of boundaries. Secondly, because of its nature, the simulation that uses periodic conditions needs the domain to be a multiple of the wave length. This is not the case with open boundary conditions (as can be seen from the fact that the green free-surface particles extend over a longer range than the blue free-surface particles). Finally, the waves propagating through a periodic domain are slightly slower than the analytical solution, which is not true for simulations with open boundary conditions as the analytical solution is imposed on the velocities.

Furthermore, let us define the error for the free-surface location using the following equation:

ε η = 1 N F S a ∈ F S (η sim -η) 2 A 2 , (58)
where FS represents the free-surface particles, N F S is the number of free-surface particles, η sim is the free-surface elevation of the simulations and η is the analytical solution given in Equation (51). These errors are plotted over time in Figure 14, which shows that for these discretisation the error (the solid line) is of the order of the particle spacing (the dashed line). Unfortunatly, this is not the case for particle size smaller than shown in Figure 14 as an unstability develops, increasing the error. These instabilities appear as a checkerboarding effect close to the inlets, and work on the volume diffusion term close to open boundaries will be necessary in the future to reduce those errors for very refined simulations.

t/T = 0 -2 0 2 η A t/T = 2 -2 0 2 η A t/T = 4 -2 0 2 η A t/T = 6 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 -2 0 2 x/L η A

Analytical solution

Periodic boundary conditions Imposed velocities In all cases the open boundary conditions seem to model properly the propagating waves, and no discrepancies appear between the fluid and imposed freesurfaces even after seven wave periods.

2D solitary wave

The Riemann invariants become particularly useful when the boundary conditions are not well known, or not easy to predict. To illustrate this concept, a test case for the propagation of a solitary wave on a bed with a step, that will cause it to break will be presented. See Figure 15.

At the inlet the solitary wave will be imposed from the Korteveg-De Vries solitary wave equation [START_REF] Dean | Water wave mechanics for engineers and scientists[END_REF]:

η(t) = A sech 2 (k(ct -x 0 )) , (59) h
(t) = D + η(t), (60) p(z, t) = ρg(z -h(t)), (61)
u x (z, t) = c η(t) h(t) , (62)
u z (z, t) = z h(t) ∂η ∂t (t). ( 63
) η(x, t) = A sech 2 (k(x -ct -x 0 )) , (64)
where A is the amplitude and x 0 is wave shifting length. The wave number and celerity are computed as k = 3A/4D 3 , c = g(A + D).

At the outlet a hydrostatic pressure corresponding to the reference level D will be imposed and the velocities will be calculated using Riemann invariants.

The geometry chosen for the test case illustrated in Figure 15 is a the reference depth D = 0.6m, wave amplitude A = D/2 and characteristic length L = 2.5D. The wave shifting length is set to ensure that the waves are not in the domain at the initialisation, i.e. x 0 = -4k.

The fluid viscosity ν is set to 10 -6 m 2 s -1 , the particle spacing will be set to one thirtieth of the amplitude (∆r = A/30) and the numerical parameters are c 0 = 20 m.s -1 and Λ = 0.1.

The evolution of the solitary wave within the domain at different instances is shown in Figure 16. It shows that the step causes the wave to break, creating a highly disturbed flow to exit the domain. Nonetheless, imposing a hydrostatic pressure and using Riemann invariants is sufficient for the wave to exit the domain without any visible reflection in WCSPH. This shows the high flexibility of these new open boundary conditions.

Finally, it should be noted that the small peak that appears in the last frame of Figure 16 appears because when the wave goes over the step, a portion of its energy is reflected back to the inlet, which in turn reflects it back towards the step. This last reflection occurs because after the wave has entered the domain the velocity is imposed as zero, making it essentially a wall.

Conclusion

The Unified Semi-Analytical boundary condition have been extended to treat complex inlets and outlets without spurious shocks. This was achieved by allowing the mass of the vertices to vary according to the imposed flow rate, and then creating or deleting particles accordingly. Furthermore, it was necessary to take into account the imposed Eulerian velocity field in order to integrate properly in time the continuity equation, and thus avoid artificial density pertubations near the boundary. Using Riemann invariants it is possible to impose compatible density and velocity fields that can support both Neumann and Dirichlet boundary condition.

These new boundary conditions have been tested on several test case. Firstly a simple flow entering and leaving the domain in a non-orthogonal way, on periodic waves where the open boundaries were used as both inlets and outlets and on a solitary waves for which the outlet boundary conditions were unknown. In all cases this formulation of open boundaries in a Smoothed Particle Hydrodynamics framework has proved to be appropriate to complex confined and free-surface flows.

Finally, the extension to 3-D should be straightforward, but since the boundary elements will be more complex greater care will be needed to allow fluid particles to enter and leave the domain without creating any artifical pertubations. and the matrix B as:

B (Y ) ≡                 u n ρ 0 • • • 0 c 2 ρ u n 0 • • • 0 0 0 u n 0 • • • . . . . . . 0 . . . 0 0 0 • • • 0 u n                 . (A.3)
Furthermore the speed of sound c is defined from the state Equation (1) as:

c ≡ ∂p ∂ρ = c 0 ρ ρ 0 ξ-1 2
.

(A.4)

The eigenvalues of B are λ -1 ≡ u nc, λ 0 ≡ u n (with a multiplicity d -1) and λ +1 ≡ u n + c. Along the open boundary there is a discontinuity between the exterior state (upon which we want to impose our values) and the interior state (i.e. the fluid domain). The eigenvalues of the problem, i , represent the slope of the characteristic curves of the Riemann problem (see Figure A.17), which are discontinuities between the states on their right and left sides (for example the wave λ 0 is a discontinuity between the data state Y 1 and Y 2 ). The characteristic waves can therefore be used to link the exterior state and the interior state. To go from one state to another the Generalized Riemann Invariants (GRI) will be used, as these define relations that hold true, for certain waves, across the wave structure. These are calculated from the relations defined by Jeffrey [START_REF] Jeffrey | Quasilinear Hyperbolic Systems and Waves[END_REF]: condition ∂p/∂n = 0 will be used to calculate the pressure on the segments, and on the vertices it will be the average of the connected segments.

When dealing with an open boundary, either the velocity or the pressure will be imposed on the vertices and Riemann conditions will be used to calculate the corresponding pressure or velocity on the segments. In all cases the missing values on the segments and vertices will be calculated from the connected vertices or segments (see Sections 3.4.1 and 3.4.2). Furthermore the mass flux will only affect the masses of vertices, even though it uses the densities, velocities and lengths of the connected segments. However, as illustrated by Figure B.18, dealing with the vertices that are in between a wall and an open boundary segments is not trivial. For the vertices that are in-between a segment with an imposed velocity and a wall, it is recommended to impose the wall velocity, and to deduce the pressure from the connected segments. For the vertices that are in between a segment with imposed pressure and a wall, it is recommended to impose the wall velocity and the open boundary pressure. Furthermore, it is recommened to keep the masses of the vertices that are inbetween wall segments and open boundaries constant. Therefore, even though these particles will be used to calculate the mass repartion factor β their masses will not evolve. This implies that the open boundary vertices can be considered as an interpolation of the Eularian fields imposed.

Imposed Velocity

This means that in the implementation 4 loops would be necessary:

1. Calculate the values that are imposed 2. Update the vertices and segments according to the connected segments or vertices

2 .Figure 1 :

 21 Figure 1: Space discretization. Vertex particles are shown in green and segments in orange.

Figure 2 :

 2 Figure 2: Mass evolution and particle creation for an inlet: the mass of vertex particles v is growing due to a positive flux us • ns at the segment s (times t 0 to t 3 ). When the masses of the inlet vertices reach the defined threshold of m ref /2 (time t 2 ), new fluid particles are created and the mass of vertex particles is decreased.

Figure 3 :

 3 Figure 3: Mass distribution process for an outgoing fluid particle a crossing a segment s connecting vertices v 0 and v 1 . The dark gray line represents the factor βa,v 0 .

4 . 1

 41 and 3.4.2 for open boundaries. 2. Time stepping of the fluid particles in the domain using Equations (27). 3. Update the mass of the vertex particles using Equation (30), create particles if m v > θ v m ref and delete particles that have crossed an open boundary.

Figure 4 :

 4 Figure 4: Non-orthogonal inlet/outlet in a square.

Figure 5 :

 5 Figure 5: Non-orthogonal inlet/outlet in a square. Evolution of the L 2 error in density (left) and velocity (right) for different discretizations without using Riemann invariants.

Figure 6 :

 6 Figure 6: Non-orthogonal inlet/outlet in a square. Evolution of the L 2 error in density (left) and velocity (right) for different discretizations using Riemann invariants with imposed velocities on the inlets and imposed pressure on the outlets.

Figure 7 :

 7 Figure 7: Configuration of the rapidly expanding pipe. The origin is placed at the beginning of the second pipe, and on the central axis of the pipes.

1 and L 2 =

 2 4W 2 . The particle spacing ∆r = W 1 /26 and the physical parameters chosen are ρ 0 = 1190 kg.m -3 , ν = 3.19 × 10 -5 m 2 .s -1 and Re = 20.6. The numerical parameters are c 0 = 0.03 m.s -1 and Λ = 0.1. A background pressure is used to help the flow stabilize.

2 (

 2 to calculate the mass fluxes at the inlet and outlet we get the following time averaged values: a) Imposed velocity at inlet (b) Computed velocity at outlet (c) Computed density at inlet (d) Imposed density at outlet

Figure 8 :

 8 Figure 8: Horizontal velocity and density profiles at the inlet and the outlet of the rapidly expanding pipe. The dashed line shows the analytical solution for a Poiseuille flow, and the solid line the simulation results.

Figure 9 :

 9 Figure 9: Streamlines in the rapidly expanding pipe from the SPH simulated flow (above), experimental results by Hammad et al. [21] (middle) and simulation result using Code_Saturne [22]. The values displayed at the streamlines are Ψ/Ψ(0, 0).

Figure 10 :

 10 Figure 10: The three regions in the rapidly expanding pipe that can be defined to claculate the head loss from the simulated flow.

Figure 11 :

 11 Figure 11: Head profile plotted along the x-axis of the rapidly expanding pipe. The blue line is the computed head, the dashed lines show the theoretical head gradient of the Poisseuille flows (from Equations 44 and 45) and the arrows around the dotted line shows the head loss from the Borda-Carnot equation.

-

  Horizontal velocity at t = 0.73s (b) Horizontal velocity at t = 1.35s velocity u x /U :

Figure 12 :

 12 Figure 12: Propagation of a regular waves on a flat bottom with open boundaries given at 2 instances. A time varying velocity profile is prescribed on both vertical boundaries so they alternate between inlets and outlets.

Figure 13 :

 13 Figure 13: Free-surface particles for simulations using different boundary conditions compared to the analytical solution of equations (51) -(57).

4 Figure 14 :Figure 15 :

 41415 Figure 14: Errors on the free-surface for simulations with an imposed velocity profile on the open boundaries for two particle spacing ∆r. The dashed line show the two values ∆r/A.Imposed VelocityImposed Pressure

Figure A. 17 :

 17 Figure A.17: Sketches of the Riemann

FluidFigure B. 18 :

 18 Figure B.18: Illustration of the different boundary conditions that can be imposed on vertices and segments. The symbol D represents a Dirichlet boundary condition and the symbol N represents a Neumann boundary condition. The superscripts show if this conditions are imposed on the velocity u or on the pressure p.

Acknowledgements

Most of the GPU fork development of GPUSPH, as well as some of the computations presented here, were run on the NVIDIA device offered as Liberski price to Martin Ferrand by the SPHERIC Steering Committee in 2010. The authors are also grateful to Arno Mayrhofer for his help during this work.

Appendix A. The Riemann problem formulation

In the present approach, the mass flux through an open boundary is calculated at the position of the segments (although it is applied to the connected vertices), therefore the values of the pressure and velocity fields on those segments are required. Since the WSPH method is used, the boundary conditions can be derived using a linearised Riemann problem. This methodology is widely used in the literature for finite volumes, but it is not as common for SPH problems. The approach chosen here follows the work of Blondel et al. [26], with simplifying assumptions of barotropic fluid and subsonic flow. The latter assumption is justified by the fact that weakly compressible SPH assumes the speed of sound to be significantly larger than the flow velocity.

Let us consider a segment s on an open boundary and work in the local reference frame of space dimension d relative to s, denoted by (n, t 1 , • • • , t d-1 ), n being the local normal unit vector and t k the tangential unit vectors. In this section, for simplicity we will drop the label s relative to all fields. The linearised Riemann problem (LRP) can now be written as:

where the unknown vector is defined as:

where Y i is a components of vector Y and V λ r,i is a components of the right eigenvector of B associated to the eigenvalue λ.

The GRI corresponding to the celerities λ -1 and λ +1 are respectively denoted by R -1 and R +1 :

where: and λ 0 = u n < 0 (outgoing mass flux). One should bear in mind that λ -1 is supposed to be always negative and λ +1 always positive as the flow is to be subsonic.

As displayed in

Therefore to link the exterior state a first assumption is made, that the wave λ -1 is a "ghost" wave, as it is outside of the domain. This means that the data Y ext of the state is assumed equal to the data Y 1 of state 1.

Wave λ 0 is known as a contact discontinuity, i.e. λ 0,1 = λ 0,2 [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF]. This implies that u n,1 = u n,2 , which is consistent with the GRI relations. Through these relations it is also considered that dρ = 0 across λ 0 . Therefore ρ 1 = ρ 2 , but also p 1 = p 2 as the pressure is defined by Equation [START_REF] Kassiotis | Coupling SPH with a 1-D Boussinesq-type wave model[END_REF].

Therefore to link the fluid velocity along the normal of the segment and the pressure, one needs only to find a relation between state 2 and the interior state; i.e. across the characteristc wave λ +1 (the tangential velocities willdiscussed later). In terms of notation, a distinction will be made between this characteristic wave within state 2, λ +1,2 , and its value within the interior state, λ +1,int . As defined in Toro [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF] the characteristic wave λ +1 can belong to three possible type of discontinuity:

1. Contact wave In this case the characteristics are parallel in both states, this means that the following condition is true:

and the GRI relations hold across the characteristic wave:

Expansion (or rarefaction) wave

In this case the characteristics are diverging, and the two data states are connected through a smooth transition. This case is defined by the following condition:

and the GRI relations hold across the characteristic wave:

(A.12)

Shock wave

In this case the characteristics are converging towards the wave λ +1 , meaning that this case is defined by the following condition:

In this case the GRI relations do not hold across the characteristic wave, and therefore the Rankine-Hugoniot relationships have to be used:

Conservation of mass:

For the case of the tangential velocities the Riemann invariants relations make it possible to write du t1 = ... = du t d-1 = 0 across the wave λ +1 (we will assume these relations still hold in the case of a shock). Therefore the tangential velocities are equal between state 2 and the interior state. However, the wave λ 0 acts as a discontinuity.

Therefore if the problem is that of an inlet then the discontinuity prevents a relationship between state 2 and state 1 to be defined for these velocities, all that is known is that du t1 = ... = du t d-1 = Const. This means that these velocities need to be defined by the user for inlet open boundaries. However for the case of an outlet then the λ 0 will be assumed to be a "ghost" wave and the exterior tangential velocities will be assumed to be equal to the tangential velocities of state 1 and 2, and therefore equal to the velocities of the interior state.

In mesh-based methods, such as finite volumes, knowing the interior state data is simple as boundary segments are only connected to one element. However in the SPH formulation there is more than one fluid particle connected to a boundary segment. Therefore defining the interior state is non-trivial. The method chosen is to use the SPH interpolations as defined by Ferrand et al. [8]:

where, as a reminder, F contains only fluid particles and α s is the Shepard filter defined for a segment s by:

Similarly, the interior velocity is defined as:

Finally, the density will be calculated using Equation [START_REF] Kassiotis | Coupling SPH with a 1-D Boussinesq-type wave model[END_REF].

Appendix B. Handling varying number of particle in a parallel implementation

Handling a varying number of particles in a parallel implementation can be difficult. To allow the algorythm presented here to be easily implemented, several key numerical considerations will be described in the following sections.

Appendix B.1. Creating and deleting particles

As mentioned earlier this algorithm has been implemented in a parallel framework, and therefore particle creation/deletion can be non-trivial. Namely, adding an extra particle needs to be an atomic operation (i.e. an operation that will be done sequentially on all processors). The same is true when deleting particles and the mass of the concerned vertex particles.

The implementation tested is based on a branch of the GPUSPH software [START_REF] Hérault | [END_REF][START_REF] Hérault | Achieving the best accuracy in an sph implementation[END_REF], and therefore most of the memory is handled in the same manner. In addition to the aforementioned problem, all the arrays are of a fixed size at the initialization. As such, changing the size of any array when a particle is removed or added can be quite expensive, especially when it happens at each iteration. Thus it is suggested to allocate larger arrays. The user has to evaluate the total number of particles that will be added during the simulation time. Regardless of this number, the loops are only done on the current number of particles at each iteration.

Furthermore it is very important that the masses of the open boundary vertices is kept constant during the stepping algorithm, so that the condition m b = m n b during the open time interval ]t n , t n+1 [ is maintained (i.e. excluding both t n and t n+1 ). The same is true for the set of particles, i.e. during the continuity equation step, the particles that have crossed an open boundary are considered, but not the particles that will be created.

Appendix B.2. Updating vertices and segments

When dealing with the semi-analytical boundary conditions, it has been established that Dirichlet boundary conditions will be imposed on the vertices, and Neumann bounday conditions on the segments [START_REF] Leroy | [END_REF]. This implies that for walls the velocity will be imposed on the vertices, and that velocity on the segments will be the average of the connected vertices. For the pressure the

Appendix B.3. Mass of the vertices

For open boundaries segments it is strongly recommended that their masses stay between ±m ref . For the vertices that are slightly above the free-surface or near walls, the SPH interpolation can create mass fluxes that will break the condition -m 0 /2 < m v < m 0 /2. This is particularly critical at the outlets, and therefore it is recommended to clip the masses and for open boundaries with imposed pressures, the mass flux above the free-surface is only allowed to be negative. Finally, above the free-surface the mass is initialised at 0.

Appendix B.4. Constraints on the time-step

The time marching scheme uses an adaptive time-stepping, where the timestep is calculated using the following condition:

where the speed of sound c 0 is chosen such that the following condition is always satisfied: Additional details about the numerical stability of WCSPH are available in Violeau and Leroy [START_REF] Violeau | [END_REF].