D Violeau
email: damien.violeau@edf.fr

Investigation on conservative and dissipative forces for Lagrangian models in CFD and application to SPH

An investigation of conservative and dissipative forces for Lagrangian Computational Fluid Dynamics is conducted from Hamiltonian considerations including energy dissipation for macroscopic systems. It is shown than discrete forces must ful…ll particular rules to be in agreement with the fundamentals of Physics. Those rules are speci…ed in the case of the Smoothed Particle Hydrodynamics (SPH) numerical approach, leading to various numerical schemes satisfying variational consistency. A clear treatment of friction forces in SPH is proposed and studied in this framework. In particular, it is proved that the kernel function, which is at the heart of interpolation in SPH, must satisfy some constraints in order to be consistent with the dissipative properties of a real ‡uid. A numerical example is given to illustrate the abovementioned considerations.

Introduction

Lagrangian methods in CFD (Computational Fluid Dynamics) are becoming more and more used in the …elds of con…ned and free-surface ‡ows. All these methods, such as SPH (Smoothed Particle Hydrodynamics, see [START_REF] Monaghan | Smoothed Particle Hydrodynamics[END_REF], DPD (Dissipative Particle Dynamics, see e.g. [START_REF] Espanol | Smoothed Dissipative Particle Dynamics[END_REF] or MPS (Moving Particle Semi-implicit method, see e.g. [START_REF] Gotoh | Lagrangian particle method for simulation of wave overtopping on a vertical seawall[END_REF] are based on the idea that a ‡ow can be seen as a collection a macroscopic particles exchanging momentum, energy, etc. In this context, the question of the discrete form of forces is crucial, and can be investigated from Lagrangian and Hamiltonian viewpoints. A key element is how to treat macroscopic systems such as ‡uid particles, for which the laws of thermodynamics predict a completely di¤erent behaviour than for individual molecules, leading to the concept of energy dissipation. Hence, the correct treatment of conservative and dissipative forces in Lagrangian method, in connection to the principles of Hamiltonian and statistical mechanics yields some rules that developers must obey when writting the discrete equations of a Lagrangian numerical method consistently with variational ideas.

In section 2, we will brie ‡y remind the backgrounds of the Hamiltonian theory, before focusing on dissipative (or friction) forces, …rst in a general context, then applied to a system of macroscopic particles. We will show that, regardless of the considered numerical approach, the discrete particle friction forces should ful…ll some requirements deduced from the conservation laws, but also from considerations of volume decrease in the phase space. In section 3, we will apply these considerations to …nd the constraints that the conservative and dissipative forces must satisfy in the context of the SPH method. This work will be based on variational considerations applied to SPH, as well as a rigorous treatment of discrete viscous forces in the context of this method. The simulation of a steady ‡ow in a periodic hill channel will illustrate our conclusions.

2 Conservative and dissipative systems

General considerations on mechanical forces

In the present section we remind a few statements regarding Lagrangian and Hamiltonian mechanics for dissipative systems. Since the pioneering work from J.-L. Lagrange (see e.g. [START_REF] Landau | [END_REF] it is generally considered that a classical mechanical system may be fully described by a Lagrangian L, playing the role of a state function depending on its N generalized co-ordinates q i and their time derivative, or generalized velocities _ q i : L = L (fq i g ; f _ q i g ; t)

(1) = E K (fq i g ; f _ q i g) U (fq i g ; t)

where E K and U are the kinetic and potential energy, respectively (here and in the following, we use the notation fq i g to denote the collection of all q i 's). Then, the Lagrange equations give the behaviour of the system fq i g as 8i; _ p i + dp i dt = @L @q i (2) = @U @q i = F cons i where the generalized momenta are de…ned by

p i + @L @ _ q i = @E K @ _ q i (3)
In (2) also appears the notion of conservative forces F cons i . The name and notation suggest that they conserve energy, as it is well known, since they can be written as the derivatives of some state function (potential energy) with respect to the co-ordinates q i , that is to say as spatial gradients of a potential. The above equations are perfectly valid, for example, for a gas of interacting (colliding) molecules, and total energy is indeed conserved at the microscopic scale. However, the behaviour of such a system involving a very large number of degrees of freedom is better understood through a statistical viewpoint. This can be achieved through the H-theorem resulting from the Boltzmann equation, which states that entropy (disorder) increases with time while the shocks between molecules spread their velocities on a larger range (see e.g. [START_REF] Reif | Fundamentals of Statistical and Thermal Physics[END_REF]. This process tranfers energy from the macroscopic to the microscopic (thermal) scale, and thus yields energy dissipation at macroscopic scales. It also tends to make the macroscopic (Boltzmann-averaged) velocities more and more uniform while time goes on. At the macroscopic scale viewpoint, where the velocities are smooth, one must then admit that some forces F diss i exist which can be quali…ed as dissipative, for they no longer conserve (macroscopic) ernergy. As a consequence, those "friction" forces cannot be written as spatial gradients, but must depend on the macroscopic velocities f _ q j g. A Taylor expansion up to order 1 gives

F diss i = X j ij (fq k g) _ q j (4)
where the ij 's are N N unknown friction (or "kinetic") coe¢ cients depending on the co-ordinates fq k g only. Next, Onsager's principle states that these coe¢ cients satisfy the following symmetry condition

ji = ij (5)
(see e.g. [START_REF] Landau | Statistical physics part 1[END_REF]. This statement allows us to rewrite (4) in the following form

F diss i = @F @ _ q i ( 6
)
where we have de…ned the following quadratic form:

F = 1 2 X i;j ij (fq k g) _ q i _ q j (7)
The Lagrange momentum equations (2) can now be reformulated, for macroscopic systems, as 8i; _

p i = @U @q i @F @ _ q i (8) = F cons i + F diss i
The system is totally de…ned at each time by the data collection (fq i g ; f _ q i g). However, (3) shows that the choice (fp i g ; fq i g) is also possible; the energy (or Hamiltonian) H is then de…ned by

H (fp i g ; fq i g ; t) = X i p i _ q i L (9)
From considerations related to Galileo's relativility principle, it is also known (see e.g. [START_REF] Goldstein | Classical Mechanics[END_REF] that the kinetic energy is a quadratic function of velocities, with coe¢ cients depending on the generalized co-ordinates as

E K = 1 2 X i;j ij (fq k g) _ q i _ q j (10)
The positivity of E K shows that the f ij g are N N coe¢ cients of a positive de…nite quadratic form, and by de…nition they are symmetric; thus the matrix of coe¢ cients f ij g has an inverse, hereafter denoted by = 1 . From ( 3) and ( 10), we then obtain linear relations between the generalized momenta and velocities as

p i = X j ij (fq k g) _ q j (11)
Combining ( 9) and ( 11) now gives the Hamiltonian as the sum or kinetic and potential energy, the former being a quadratic function of generalized momenta:

H (fp i g ; fq i g ; t) = E K (fq i g ; f _ q i g) + U (fq i g ; t) (12) = 1 2 X i;j ij (fq k g) p i p j + U (fq i g ; t)
where the f ij g are the coe¢ cients of the matrix . Using (3) and ( 8), the di¤erential of H can be written as

dH = X i ( _ q i dp i + p i d _ q i ) @L @t dt X i @L @q i dq i + @L @ _ q i d _ q i (13) = @L @t dt + X i _ q i dp i _ p i + @F @ _ q i dq i
from which it immediately follows the set of Hamiltonian (macroscopic) dissipative equations:

@H @t = @L @t (14) @H @p i = _ q i @H @q i = _ p i @F @ _ q i
From (13) also follows the time behaviour of energy :

dH dt = @L @t X i @F @ _ q i _ q i (15)
When the system is isolated from any external in ‡uence, the Lagrangian is independant of time (@L=@t = 0). Moreover, the quadratic form (7) of the function F …nally yields the following relation:

dH dt = 2F (16)
For this reason, F is usually called the dissipation function. The fact that the macroscopic energy decreases with time immediately yields the positiveness of F , meaning that the quadratic form de…ned by the f ij g is positive de…nite. The dissipation of energy has an important geometrical consequence. Consider a small ensemble ~ of initial conditions (fp i g ; fq i g) for the system in the socalled phase space, and let us call Ṽ the volume of this ensemble. Let us look at the evolution of the system along time for all these possible initial values. In the phase space, the subsequent values of (fp i g ; fq i g) will draw an ensemble of curves following the "velocitiy …eld" (f _ p i g ; f _ q i g) according to the set of equations ( 14). We can now a¤ect to each point (fp i g ; fq i g) a probability

, considered as a "density" in the phase space. The conservation of total probability (which equals unity) is then formally equivalent to the conservation of mass of a virtual ‡uid in the phase space. For a real ‡uid of local volume V , density and velocity u in the physical space of spatial co-ordinates x i , mass conservation reads:

1 V dV dt = 1 d dt = div u = X i @u i @x i (17)
Transposing ( 17) with the notations of the phase space (i:e. V; ; r; u ! Ṽ ; ; (fp i g ; fq i g) ; (f _ p i g ; f _ q i g)) and using ( 4) and ( 14), we get

1 d dt = X i @ _ p i @p i + @ _ q i @q i (18) = X i @ 2 H @p i @q i @ 2 F @p i @ _ q i + @ 2 H @q i @p i = X i;j @ ij _ q j @p i
We mentioned earlier that the ij 's are functions of the q i 's only. Thus, with ( 14) and ( 12), we obtain

1 Ṽ d Ṽ dt = X i;j ij @ 2 H @p i @p j (19) = X i;j ij ij = tr ( )
where is the N N matrix of coe¢ cients ij . Equation ( 19) gives the time evolution of the volume e V of ~ in the phase space. Now, reminding that the trace of the product of two positive de…nite matrices is positive, one comes to the conclusion that the volume e V in the (macroscopic) phase space decreases under the e¤ect of dissipative forces. In the absence of dissipative forces, equation ( 19) turns back to Liouville's theorem, stating that the volumes are preserved in the (microscopic) phase space. On the other hand, at the macroscopic scale, the larger the friction forces, the more the volumes are decreased in the phase space. This means that the paths (fp i g ; fq i g) of the mechanical system get closer and closer, tending to an attractor, which is a known feature of dissipative systems. This deserves an interpretation: while disorder increases at the microscopic scale, the resulting friction at large scales, making the macroscopic velocities more and more homogeneous, increases macroscopic order, so that the qualitative behaviour of the system in the phase space depends on the scale at which we investigate it. In other terms, the measure of the volume in the phase space depends on the spatial scale of observation, as it is for a fractal object.

It is also known that the total linear and angular momenta of a closed system are both conserved, with or without dissipative forces. These conservation laws, according to Noether's theorem, are directely connected to the Hamiltonian form of the equations of motion [START_REF] Landau | [END_REF].

Systems of particles

Extending the above considerations to the case of a system of N p macroscopic particles labeled by the letters a, b, etc. is straightforward, the number of degrees of freedom being now N = nN p , where n is the dimension of the physical space (generally equal to 2 or 3, but sometimes equal to 1). Knowing the particle masses fm a g, their respective positions fr a g and velocities fu a = _ r a g (now playing the role of the fq i g and f _ q i g, respectively) and internal energies fE int;a g, the lagrangian (1) reads

L = X a 1 2 m a u 2 a X a E int;a U ext (fr a g) (20)
where the kinetic and potential energies are de…ned by

E K = X a 1 2 m a u 2 a = 1 2 X a;b u T a ab u b (21) U = X a E int;a + U ext (fr a g)
The kinetic energy is in the form (10) with n(a 1)+k;n(b 1)+l = m a kl ab (22) where k and l here denote spatial directions (from 1 to n), so that n (a 1)+k goes from 1 to N = nN p , as required. We may also, as pointed out by (21), write the kinetic energy of a particle using N 2 p matrices ab of order n attached to each pair of particles (one must insist on the fact that labels a and b do not represent here matrix indices):

ab = m a ab ( 23
)
ab being a Kronecker matrix of order n, i.e. ab = I n (the n th order identity matrix) if a = b, otherwise ab = 0. Note that here, the coe¢ cients ij 's do not depend on fq k g (in other words, the matrices ab are independant on the particle positions fr c g). The equations of motion without friction forces (2) may be written in a vector form for each particle as 8a;

dp a dt = @L @r a (24)

= @E int;a @r a @U ext @r a = F int a + F ext a
where the momenta are de…ned, according to (3) and ( 20), by

p a = @L @u a (25) = m a u a
The conservative force is thus now made of two parts, namely internal and external forces:

F cons a = F int a + F ext a ( 26
)
For a system of particles representing a ‡uid, the internal forces F int a correspond to pressure forces, as we will see later (section 3.1), while the external force F ext a is generally restricted to gravity:

U ext (fr a g) = X a m a g r a (27) = X a m a gz a
where g is the gravity acceleration (with norm g) and fz a g the particle elevations above an arbitrary reference level. The above two forces are conservative

(F int a + F ext a = F cons a
), since they are written in a gradient form. On the other hand, particle dissipative forces may now be taken into account, similarly to (4) and [START_REF] Bonet | Variational and momentum preservation aspects of Smoothed Particle Hydrodynamics formulations[END_REF], by the following model:

F diss a = X b ab (fr c g) u ab (28)
where

u ab = u a u b = u ba (29)
and the ab 's are N 2 p matrices of order n attached to each pair of particles (like with the matrices ab de…ned above, a and b do not represent matrix indices put particle labels) depending on the collection of particle positions fr c g and satisfying ba = ab (30) according to Onsager's principle (5) regarding kinetic coe¢ cients. The main di¤erence with (4) is that now the dissipative forces depend on the velocity di¤ erences. The reason is that a set of particles moving with a uniform velocity does not expericence any friction (in a continuous media viewpoint, the latter proposition states that shear forces depend on velocity gradients, see section 3.2). The friction force (28) can be separated into individual forces due to each particle b, namely F diss b!a = ab (fr c g) u ab . Since such a force can only involve the positions of the particles a and b, only r a and r b must appear in ab among the list fr c g. Besides, the friction force vanishes when a = b, so that ab must be a function of r ab , de…ned by

r ab = r a r b = r ba (31)
Lastly, for isotropy reasons the shear force F diss b!a do not depend on the orientation r ab . Hence, we must have

ab = ab (r ab ) (32)
where r ab = jr ab j, which is consistent with (30). The latter symmetry, as in [START_REF] Bonet | Variational and momentum preservation aspects of Smoothed Particle Hydrodynamics formulations[END_REF], allows to write the friction force experienced by particle a as

F diss a = @F @u a (33)
where

F = 1 2 X a;b u T ab ab u ab (34)
(we will sometimes drop the explicit dependency of ab in r ab , for simplicity). With these premisses, the dissipative particle momentum equation appears in a form analogous to (8):

8a; dp a dt = @E int;a @r a @F @u a @U ext @r a (35) = F int a + F diss a + F ext a
The conservation of total linear momentum of a closed ensemble of interacting particles directly stems from the fact that the sum of all forces is equal to zero. In particular, from (29) and (30) it immediately follows that the actionreaction law is ful…lled by particle friction forces:

F diss b!a = ab u ab (36) = F diss a!b
Similarly to (3), and using (20), ( 21) and ( 23), one de…nes the particles momentum vectors by

p a = @L @u a = X b ab u b = mu a (37)
The Hamiltonian now becomes

H (r a ; p a ) = 1 2 X a;b p T a ab p b + X a E int;a + U (fr a g) (38)
where the ab 's are N 2 p matrices of order n de…ned by

ab = 1 ab = 1 m a ab (39)
Equation ( 38) is the particle analogy to ( 12). The law of energy dissipation (15), for a closed set of particles, now reads

dH dt = X a @F @u a u a (40)
However, rearranging (40) to give energy losses as a function of F is now slightly di¤erent from ( 16). Indeed, using (28) and swaping the dummy labels a and b, we get

dH dt = X a;b ( ab u ab ) u a (41) = X a;b ( ab u ab ) u b
(we used the symmetry laws ( 29) and ( 30)). Then, taking the averaged of the former two expressions, one comes to

dH dt = 1 2 X a;b u T ab ab u ab (42) = F
Comparing to ( 16), we see that the dissipation rate is now assumed by F (instead of 2F ); this is due to the fact that the friction forces are linear function of the velocity di¤ erences, as mentioned above. The energy decrease then requires each matrix ab to be positive de…nite.

We will now investigate the appropriate form for the matrices ab . We already mentioned that shear forces satisfy the conservation of linear momentum (equation ( 36)). This is a direct consequence of the fact that these forces take their origin in additive microscopic phenomena which obey this conservation law. Similarly, the angular momentum having the same additive property at the microscopic scale, one must immediately deduce that shear forces should ‡ul…ll the condition of total angular momentum conservation, which is achieved by imposing that the angular momentum of individual shear forces are antisymmetric with respect to particle labels:

r a ( ab u ab ) = r b ( ba u ba ) (43)
or, using the symmetry properties:

r ab ( ab u ab ) = 0 (44)
The condition ( 44) is ful…lled only if the vector ab u ab is collinear to r ab ; thus there must exist a scalar quantity ab symmetric with respect to a et b, such that ab u ab = ab r ab . Besides, the latter relation shows that ab must be linearly depending on u ab , hence of the form ab = ab u ab , where ab is an antisymmetric vector. This gives

ab = r ab ab (45) ba = ab
The above condition is not su¢ cient to properly model shear forces in a particle formalism. In addition to the conservation of total angular momentum, these forces should vanish for a set of particles moving according to a rigid body rotation, thus for a velocity distribution de…ned by u a = u + r a ! where u and ! are arbitrary constant vectors, i.e. u ab = r ab !. Coming back to the general form (28), the latter condition is satis…ed if

8!; X b ab (r ab !) = 0 (46)
By similar arguments, it can easily be proved that this requires the existence of vectors 0 ab = 0 ba such that

ab = 0 ab r ab (47)
Conditions ( 45) and ( 47) are simultaneously ful…lled if and only if there exist symmetric scalar coe¢ cients ab so that ab = ab r ab = 0 ab . Finally, the matrices ab can be written as

ab = ab r ab r ab (48) ba = ab
With this formula, the particle friction forces (28) take the form

F diss a = X b ab (r ab ) (u ab r ab ) r ab (49)
(we remind that the ab 's, and thus the ab 's, are functions of the particle distance r ab ). The dissipation function (34) now reads

F = 1 2 X a;b ab (r ab ) (u ab r ab ) 2 (50)
As required, both the shear forces (49) and the energy dissipation (50) vanish for a rigid body motion (u ab r ab = 0). One could be surprised that the friction force exerted on a by b takes the form F diss b!a = ab (u ab r ab ) r ab , thus is on the axis r ab , contrary to what the theory of continuous viscous ‡uids suggests, i.e. a shear force collinear with u ab (see section 3.2). The latter case would involve matrices ab proportional to the identity I n , but would violate the conservation of total angular momentum and the non-dissipative character of rigid motions.

In order to ensure the de…nite positiveness of the ab 's, the ab 's must be positive non zero quantities. It is apparently surprinsing to …nd that, for b = a, r aa is simply the null vector 0, thus the aa 's seem to be null matrices, contrary to the what latter condition states. However, one should keep in mind that the ab 's are functions of the distances r ab . From ( 48), it appears that the nonnullity of the matrices aa 's can be obtained provided ab r ab r ab tends to a non-zero constant while r ab tends to zero, i.e. if ab varies according to r 2 ab when r ab tends to zero. Thus, there must exist a function g ab (r ab ) so that

ab = g ab (r ab ) r 2 ab g ab (r ab ) ! r ab !0 g a > 0 (51)
where the g a 's are positive quantities, constant whith respect to the particle distance (however, g a can depend on local thermodynamical properties of the particle a, like its temperature). Formulae ( 49) and (50) now read

F diss a = X b g ab (r ab ) (u ab e ab ) e ab (52) F = 1 2 X a;b6 =a g ab (r ab ) r ab (u ab e ab ) 2
where we have used unit vectors

e ab = r ab r ab = e ba (53)
To …nish this section, let us examine how the particle friction forces in ‡uence the volume in the macroscopic phase space. According to ( 19) and ( 39), one can write

1 Ṽ d Ṽ dt = X a;b tr ( ab ab ) (54) = X a tr aa m a
where the matrices are of order n (i.e. the space dimension), contrary to (19) where those were of order N (i.e. the number of mechanical parameters). With (54), it appears even more clearly than with (19) that the dissipative property of friction forces (or the positiveness of all the matrices ab ) is responsible for a decrease of volumes in the phase space. If the matrices aa were null, the decrease of volumes in the macroscopic phase space would not be ful…lled by the model. With (51), we …nally get

1 Ṽ d Ṽ dt = X a 1 m a lim b !a ab (r ab ) r 2 ab (55) = X a g a m a < 0
(keeping in mind the positiveness of masses).

3 Application to the SPH numerical method

Compatibility of discrete SPH operators

At the limit where the size of the particles tends to zero while their number tends to in…nity, the above equations lead to the Lagragian theory of continuous ‡uids (see e.g. [START_REF] Kambe | Variational formulation of ideal ‡uid ‡ows according to gauge principle[END_REF]. However, a Lagrangian computational approach can be directly constructed from the discrete viewpoint developed in the previous section. We will now investigate how the above principles can serve to write the discrete equations of the Smoothed Particle Hydrodynamics (SPH) numerical method in relevant forms. SPH is today widely used in CFD (see [START_REF] Monaghan | Smoothed Particle Hydrodynamics[END_REF] for a detailed review). It is not within the scope of this paper to repeat or summarize the numerous works already made on the accuracy of SPH interpolants and operators; we will …rst give a rapid overview of these principles. SPH is based on the idea that a ‡ow may be modelled using a collection of macroscopic particles fag following the ideas and notations of section 2. However, writing the equations of motion of a ‡uid requires further assumptions. Discrete approximations of continuous di¤erential operators are performed through the following process. For any scalar …eld A, the ensemble of its values at the points occupied by the particles is denoted by fA a g. First, being the ‡uid domain, A a is approximated by

A a = A (r a ) = Z A (r) (r r a ) d n r (56) Z A (r) w h (r r a ) d n r X b V b A (r b ) w h (r ab )
The …rst line of equation ( 56) is exact. The Dirac distribution is then approximated by a kernel regular function w h , then the integral is approximated by a discrete (Riemann) summation over the particles b, where fV a g are the volumes of the particles. The kernel is generally compactly supported, so that the discrete sum involves a limited number of neighbour particles. Besides, it is proved (see e.g. [START_REF] Oger | An improved SPH method: Toward higher order convergence[END_REF] that the approximation ( 56) is better if w h depends only on the particle distance, i.e. w h (r) = w h (r). We thus introduce a positive dimensionless kernel f such that

w h (r) = w;n h n f (q) ( 57
)
where q = r=h, h being the so-called smoothing length and w;n a dimensionless constant depending on the choice of the kernel and the space dimension n (usually equal to 2 or 3). The kernel is also usually required to be normalized, i.e. to satisfy

Z w h (r) d n r = 1 (58)
Under this assumption, considerations of spatial symmetry easily show that all momenta of odd orders of _ w h vanish, i.e. for all integers m:

Z _ w h (r) r r ::: r | {z } 2m+1 d n r = 0 (59)
(in our notations, for any set of m vectors A, B, ... Z, the quantity A B ::: Z is a tensor of order m whose components are A i1 B i2 :::Z im ). The constant w;n is chosen such that the normalizing condition (58) is ful…lled, which gives

w;n = 1 S n R R f 0 f (q) q n 1 dq (60)
where R f is the radius of the support of f and S n is the surface of the unit sphere in dimension n. We consider here the following example of kernel, namely the spline of order 5 given by

f 5 (q) = 8 > > < > > : (3 q) 5 6 (2 q) 5 + 15 (1 q) 5 if 0 q 1 (3 q) 5 6 (2 q) 5 if 1 q 2 (3 q) 5 if 2 q 3 0 if 3 q (61)
with support of radius R f = 3, while (60) gives

w;2 = 7 478 (62)
w;3

= 1 120

Let us now comme back to the SPH interpolation. Denoting w ab + w h (r ab ), (56) reads

A a X b V b A b w ab (63)
As an example, the density f a g of the particles are de…ned by

a = m a V a (64)
and can be interpolated using (63) to give

a X b m b w ab (65)
The assumption regarding the kernel dependency on r ab also yield the following rule:

@w h (r ab ) @r a = _ w h (r ab ) e ab ( 66
)
where

_ w h = w;n h n+1 _ f (q) (67)
is the …rst derivative of w h and e ab is de…ned by ( 53). Introducing the notation _ w ab + _ w h (r ab ), we can write

w ab = w ba (68) _ w ab = _ w ba
Then, taking the gradient of (63) with respect to r a and using (66) lead to a discrete approximation of the gradient of A:

(grad A) a X b V b A b _ w ab e ab (69)
It is known that such considerations can be applied to any kind of tensor. We will thus introduce discrete gradient (G) and divergence (D) operators of scalar (fA b g) and vector (fA b g) …elds at the location of the particle a as

G a fA b g + X b V b A b _ w ab e ab (grad A) a (70) D a fA b g + X b V b A b _ w ab e ab (div A) a
The formulae (70) are not the only possible ones allowing gradient and divergence discrete approximations. We may …rst observe that for every real number k, the following exact rules stand:

grad A = k grad A k + A k grad k (71) div A = k div A k + A k grad k
Writing the above two expressions at the location of the particle a, and using G and D to approximate the continuous operators, we get

(grad A) a k a G a A k b + A a k a G a k b (72) (div A) a k a D a A k b + A a k a G a k b
Then, using the de…nitions (70), we …nd a family of discrete gradient (G k ) and divergence (D k ) operators:

G k a fA b g + X b V b 2k b A a + 2k a A b ( a b ) k _ w ab e ab ( 73
)
D k a fA b g + X b V b 2k b A a + 2k a A b ( a b ) k _ w ab e ab
(see e.g. [START_REF] Monaghan | Smoothed Particle Hydrodynamics[END_REF][START_REF] Oger | An improved SPH method: Toward higher order convergence[END_REF]. One could also use the following rules:

grad A = 1 k grad k A A k grad k (74) div A = 1 k div k A A k grad k
to get another family of discrete operators as:

Gk a fA b g + 1 2k a X b V b ( a b ) k A ab _ w ab e ab ( 75
)
Dk a fA b g + 1 2k a X b V b ( a b ) k A ab _ w ab e ab
Under the assumptions of isotropy and normalization of the kernel (in particular, from the important consequence (59)), the various operators introduced above are accurate up to order 2 regarding h (see e.g. [START_REF] Oger | An improved SPH method: Toward higher order convergence[END_REF]. They have di¤erent properties which make them more or less relevant to model the equations of motion in a discrete form. Let us start with the continuity equation for a compressible ‡uid (17). Using the operator Dk (equation ( 75)) to approximate div u, we immediately …nd the following discrete form giving the evolution of particle densities along time:

8a; d a dt = 1 2k 1 a X b V b ( a b ) k u ab _ w ab e ab (76)
Equation ( 76) is not the only way to compute particle densities. Coming back to the initial SPH interpolation (63), one can write

Both ( 76) and ( 78) can be used to estimate the densities f a g. Although these approaches seem di¤erent, they are consistent to each other, which can be seen by writing the Lagrangian time derivative of (78), using (64):

8a; d a dt = a X b V b b a k dw ab dt (79)
Next, one can write the time derivative of w ab remembering that it only depends on r ab and using ( 66) and ( 68):

dw h (r ab ) dt = @w ab @r a dr a dt + @w ab @r b dr b dt (80) = _ w h (r ab ) e ab u ab
Equation ( 79), together with (80), gives back (76). In the following, we will respectively refer to these methods as the discrete continuity equation ( 76) and the density interpolation equation ( 78). The continuity equation has been proved to provide better prediction in the vicinity of a free surface [START_REF] Monaghan | Simulating free surface ‡ows with SPH[END_REF]. In addition, it has the advantage that a constant discrete velocity …eld fu a g (i.e. u ab = 0) keeps the density constant with time, which is not ensured by the density interpolation method.

We will now focus on the inviscid (i.e. conservative) momentum equation, written in its continuous form as

du dt = 1 grad p + g (81)
Approximating the pressure gradient with the operator G k (equation ( 73)), one obtains a discrete equation of motion in the form 8a;

du a dt = X b m b 2k b p a + 2k a p b ( a b ) k+1 _ w ab e ab + g (82)
The SPH equation of motion thus takes the form of a discrete particle momentum equation involving individual forces, just as in section 2. The choice of G k instead of Gk presents the advantage that the pressure (i.e. internal) forces are now written as

F int b!a = m a m b 2k b p a + 2k a p b ( a b ) k+1 _ w ab e ab = F int a!b (83)
The last equality stems from the symmetry rules ( 53) and (68). Thus, the internal forces satisfy the conservation of total linear momentum of a closed set of particles, as required by the theory. It is important to note that this property would not be true with Gk . Moreover, forces (83) are collinear with e ab and hence satisfy the conservation of total angular momentum (see e.g. [START_REF] Bonet | Variational and momentum preservation aspects of Smoothed Particle Hydrodynamics formulations[END_REF].

As a matter of fact, the above discrete conservation law suggests that the SPH equation of motion (82) may be derived from an action principle, which has been known for a long time for some particular forms of this equation (see e.g. [START_REF] Monaghan | Smoothed Particle Hydrodynamics[END_REF]. We will now show that the general form (82) of the momentum equation is consistent with the proposed formulations for the density estimation, i.e. equations ( 76) or ( 78), from a discrete variational principle. We start from the Lagrange equations ( 24), with the Lagrangian (20) and the form (27) for the potential of external forces. Introducing the internal energy per unit mass e int;b = E int;b =m b , we can write

@L @r a = X b m b @e int;b @r a X b m b g @z b @r a (84) = X b m b @e int @ b @ b @r a X b m b g ab e z = X b m b p b 2 b @ b @r a + m a g
where we used the well-known thermodynamical relation

@e int @ = p 2 (85)
which conducts from microscopic to macroscopic scales. It appears that internal forces depend on the derivative of the density of b with respect to the co-ordinates of a, i.e. @ b =@r a . The calculation of these terms depends on whether (76) or ( 78) are used to calculate the density; with the density interpolation, b explicitely depends on all the positions fr c g via the quantities w bc , while with the continuity equation b depends on the positions through the velocities u bc . We start our demonstration from the density interpolation (78). Using (66), we get Coming back to (87) and using ( 25), the Lagrange equation ( 24) …nally takes the following form:

@ b @r a = X c m c c b k 1 @w bc @r a (
8a; du a dt = X b m b 2k b p a + 2k a p b ( a b ) k+1 _ w ab e ab + g (89)
which is exactly (82). As a consequence, the operator G k used for pressure gradient is consistent with the density interpolation (78). We now consider the discrete continuity equation ( 76) for the density estimation; multiplying this equation by dt and noting that dr a = u a dt, one can write the di¤erential d b as

d b = 1 2k 2 b X c m c ( b c ) k 1 (dr b dr c ) _ w bc e bc (90) = X a @ b @r a dr a
(the …rst line of (90) can also be obtained multiplying (79) by dt and using (80)). The last line of ( 90) is nothing else than the de…nition of the derivatives @ b =@r a ; identifying the terms of both expressions, one gets

@ b @r a = 1 2k 2 b X c m c ( b c ) k 1 ( ab ac ) _ w bc e bc (91)
We may now insert (91) in (84) and using (68), one obtains:

@L @r a = X b;c m b m c p b 2k b ( b c ) k 1 ( ab ac ) _ w bc e bc (92) = X c m a m c p a 2k a ( a c ) k 1 _ w ac e ac + X b m b m a p b 2k b ( b a ) k 1 _ w ba e ba = X b m a m b 2k b p a + 2k a p b ( a b ) k+1 _
w ab e ab which goes back again to (82) (we intentionally dropped m a g here for simplicity). This shows the relevance of the choice made in ( 76) and ( 82) for the discrete form of the divergence ( Dk ) and gradient (G k ) operators, respectively. We may say that the discrete operators G k et Dk are compatible in a variational acceptation. Similarly, it can be shown that Gk et D k are compatible in the same way; however, Gk obviously does not ful…ll the antisymmetry principle (83), and thus violate the standard conservation laws, as already pointed out by [START_REF] Bonet | Variational and momentum preservation aspects of Smoothed Particle Hydrodynamics formulations[END_REF] for the particular case k = 0, although it was shown by numerical experience that the operator G0 provides satisfactory results [START_REF] Oger | An improved SPH method: Toward higher order convergence[END_REF]. As a …rst conclusion, we may recommend the choice of the operator Dk for the continuity equation (or, alternatively, equation ( 78) as a density interpolant) together with G k for the pressure gradient. With the choice k = 0 and k = 1, respectively, they give the following two sets of equations: Note that the interpolation form allowing the estimation of a in (93) may be understood as a numerical scheme: the density of particle a at time n, denoted by Both momentum equations ( 93) and ( 94) are known by SPH modellers, and it was already pointed out by [START_REF] Bonet | Variational and momentum preservation aspects of Smoothed Particle Hydrodynamics formulations[END_REF] that the choice of one of them should be made together with the correct equation for density estimation. However, these authors' developments lead to the (erroneous) conclusion that the momentum equation ( 93) is only compatible with the continuity equation while the momentum form (94) should be used together with the density interpolation. On slightly di¤erent backgrounds, [START_REF] Vila | On particle methods and Smooth Particle Hydrodynamics[END_REF] came to the same conclusions regarding the use of the continuity equation with k = 0 and k = 1, while [START_REF] Oger | An improved SPH method: Toward higher order convergence[END_REF] extended Bonet and Lok's approach to arbitrary k's, but they do not draw any conclusion regarding the density interpolation method. Finally, [START_REF] Monaghan | Smoothed Particle Hydrodynamics[END_REF] presents a variational formulation with the continuity equation and k = 1 extended to variable smoothing lengths, while [START_REF] Bonet | Hamiltonian formulation of the variable-h SPH equations[END_REF] add interactions with boundaries. Besides these developments, our calculations show that both the continuity equation and the density interpolation (with any value of k) can be used together with both forms of the momentum equation, provided k is properly chosen. Any mixed form (e.g. k = 0 for the continuity equation and k = 1 for momentum) is less relevant in terms of variational mechanics, and thus seems to be avoided in SPH applications, but only the numerical experience will tell us whether this viewpoint is correct (see below).

d a dt = a X b V b u ab _ w ab e ab or a = a X b V b w ab (93) du a dt = X b m b p a +
One should emphasize the fact that the compatibility of operators Dk and G k directly stems from the Hamiltonian form of discrete (particle) mechanics as presented in section 2.2. However, it can be understood in another way. Recent developments on truly incompressible SPH (see e.g. [START_REF] Lee | Comparisons of weakly compressible and truly incompressible SPH algorithms for 2D ‡ows[END_REF] have shown that important properties of the continuous di¤erential operators must be satis…ed in their discrete form. One of the most important of the gradient and divergence properties is their adjointness. Indeed, in the space of continuous pressure (p) and velocity (u) …elds, one can de…ne scalar products as

(p; q) = Z pq d (95) hu; vi = Z u v d
We then obtain

(p; div u) + hgrad p; ui = Z p div u d + Z (grad p) u d (96) = Z div (pu) d = I @ pu n d = 0
(the last integral vanishes due to boundary conditions u n = 0 on solid walls and p = 0 on free surfaces). In the discrete space of pressure (fp b g) and velocity (fu b g) …elds, (96) can be rewritten as

(fA b g ; fB b g) = X a V a A a B a (97) hfA b g ; fB b gi = X a V a A a B a
Using the symmetry rules ( 53) and ( 68), we can then calculate

G k a fA b g ; fB b g = X a;b m b m a 2k b A a + 2k a A b ( a b ) k+1 B a _ w ab e ab (98) = X a;b m a m b 2k b A a + 2k a A b ( a b ) k+1 B b _ w ab e ab = 1 2 X a;b m a m b 2k b A a + 2k a A b ( a b ) k+1 B ab _
w ab e ab (the dummy indices a and b have been swapped in line 2; line 3 is the average of the …rst two lines). Similary:

fA b g ; Dk a fB b g = X a;b m b m a 2k b A a ( a b ) k+1 B ab _ w ab e ab (99) = 1 2 X a;b m b m a 2k b A a + 2k a A b ( a b ) k+1 B ab _
w ab e ab which is obviously the opposite to ( 98). This gives a discrete and exact form of the relation ( 96), meaning that Dk and G k are skew adjoint. This is consistant with the variational form of continuous ‡uid dynamics, in which the relation (96) plays an evident role (see e.g. [START_REF] Kambe | Variational formulation of ideal ‡uid ‡ows according to gauge principle[END_REF], for both compressible and incompressible ‡ows. From these considerations, it would be possible to derive the general momentum equation ( 82) from a Lagrangian function L including the continuity equation in its form (76) as a constraint, using the technique of Lagrange multipliers (see [START_REF] Monaghan | Smoothed Particle Hydrodynamics[END_REF] for the case k = 0).

The conclusion of the above analysis is that one should carefully choose SPH operators; however, numerical tests leads to slightly di¤erent conclusions, as we will now see. To illustrate the ability of the presented SPH interpolants to predict velocity …elds, we now use the test case of the 2-D periodic steady laminar hill ‡ow based on the geometry presented on …gure 1 (the detailed geometry is available in [START_REF] Uribe | 10 th Ercoftac (SIG15) / IAHR / QNET-CFD Workshop on Re…ned Turbulence Modelling[END_REF]. The distances are H = 84:98 mm, L = 252 mm, l = 54 mm, h = 28 mm (obviously, the latter should not be confused with the smoothing length) and h I = 56.98 mm. About 20,000 ‡uid particles are driven by a horizontal propelling force updated in time to prescribe a constant mean velocity U = 1.785 10 3 m/s. The Reynolds number, based on the hill height h and the mean velocity U , is equal to 50, and there is no free surface, the upper boundary being a solid wall (hence, gravity is not considered here). In order to compute the pressure, a state equation is used (see [START_REF] Monaghan | Simulating free surface ‡ows with SPH[END_REF], based on a numerical speed of sound c 0 = 0:06 m/s (so that the ‡ow keeps almost incompressible). The results consist here of four vertical velocity pro…les at x=h = 0.5, 3.0, 5.0 and 8.0, respectively; they will be compared with the …nite volume simulations of [START_REF] Uribe | 10 th Ercoftac (SIG15) / IAHR / QNET-CFD Workshop on Re…ned Turbulence Modelling[END_REF]. Figure 2 shows that SPH presents a very good agreement to …nite volumes. It may be emphasized that the Lagrangian nature of SPH yields small ‡uctuations that have been removed here using time-averaging. We carried out various simulations with di¤erent options among the models proposed above; our numerical experience yields the following conclusions:

-The velocity pro…les are not sensitive to the choice of the exponent k; -Although our analysis suggests that the exponent k should be the same in the continuity equation and in the momentum equation (system (93) or system (94)), mixing the exponents (e.g. using the continuity equation of ( 93) and the momentum equation of ( 94)) does not a¤ects the results;

-While all the computations based on the use of a continuity equation provide good results, the use of a density interpolant like (78) induces instabilities and creates voids (see …gure 3). These voids remain even when increasing the speed of sound. Besides, even before the voids appear, the velocity pro…les are not as good as with a continuity equation (…gure 4). Thus, it seems that some of the conclusions from [START_REF] Bonet | Variational and momentum preservation aspects of Smoothed Particle Hydrodynamics formulations[END_REF] should be revised: the reason of the low quality of one of their computations is likely to be due to the use of a density interpolant, the continuity equations being more stable.

These conclusions may be modi…ed under di¤erent conditions, particularly for unsteady ‡ows or gravity ‡ow. However, this example shows that any conclusion based on theoretical considerations should always be subject to a numerical veri…cation. In other words, we have proved that various SPH interpolants can be deduced from variational considerations when satisfying some rules of consistency, but the numerical behaviour of the SPH method seems to be less demanding and allows more choices that the present theory suggests.

SPH dissipative forces

We now come to SPH approximations of dissipative forces, using the considerations developed in section 2.2. We …rst remind the general continuous form of the momentum equation for viscous ‡uids. The extension of Euler's momentum equation ( 81) to friction forces, namely the Navier-Stokes momentum equation, can be written as

du dt = 1 grad p + 1 div + g (100)
with

= s D + (div u) I (101) div = div ( s) + grad 2 3 div u s = 1 2 h grad u + (grad u) T i
where and are the ‡uid viscosities, s the rate-of-strain tensor and the exponent D denotes its deviatoric part (see e.g. [START_REF] Landau | Fluid Mechanics[END_REF]. We must then extend our SPH considerations to discrete operators of order 2.

It is known that using second order derivatives of the kernel leads to unphysical di¤usion properties [START_REF] Monaghan | Smoothed Particle Hydrodynamics[END_REF] and to numerical instabilities. It is thus recommended to …nd SPH interpolands of …rst order, through the following scheme. We …rst give an idea of the typical SPH form of such operators by focusing on the vector [div (K A grad A)] a , where K A is a scalar di¤usion coe¢ cient relative to the vector …eld A. Using a vectorial form of the operator D k (equation ( 73)) to estimate grad A gives the following approximation:

D k a f(K A grad A) b g = X b V b 2k b K A;a (grad A) a + 2k a K A;b (grad A) b ( a b ) k _ w ab e ab (102)
The gradients appearing in (102) could be estimated with operators G k or Gk (see e.g. Bonet and Rodriguez-Pas, 2005), but this method would lead to double discrete summations, hence to a too demanding numerical algorithm. We will then approximate the gradients by Taylor-expanding A b to the …rst order around r a :

A b = A a (grad A) a r ab + O r 2 ab (103)
(the "minus" sign stems from the de…nition r ab = r a r b ). We thus obtain

(grad A) a e ab A ab r ab (104)
Swaping indices in (104) also gives

(grad A) b e ab = (grad A) b e ba (105) A ba r ba = A ab r ab
Introducing ( 104) and ( 105) into (102) gives a new family of "Laplacian" operators:

L k a [fK A;b g ; fA b g] + X b V b 2k b K A;a + 2k a K A;b ( a b ) k A ab r ab _ w ab [div (K A grad A)] a (106) For k = 0, we …nd L 0 a [fK A;b g ; fA b g] + 2 X b V b K A;ab A ab r ab _ w ab (107)
with, for any arbitrary …eld B:

B ab = B a + B b 2 (108)
However, this method can hardly be extended to arbitrary second derivatives. Following [START_REF] Monaghan | Smoothed Particle Hydrodynamics[END_REF], we will use the following general discrete approximation of the ij component of the tensor [grad (K A grad A)] a :

@ @x i K A @A @x j a X b V b K A;ab A ab r ab [(n + 2) e ab;i e ab;j ij ] _ w ab (109)
where e ab;i is the i th component of the vector e ab (n is still the dimension of the problem). This form was …rst suggested by [START_REF] Espanol | Smoothed Dissipative Particle Dynamics[END_REF] for a similar Lagrangian numerical method (Smoothed Dissipative Particle Dynamics) in the case n = 3 and constant K A , and extended to the cases n = 2 and n = 3 for SPH and arbitrary K A by [START_REF] Monaghan | Smoothed Particle Hydrodynamics[END_REF]. However, the latter paper does not provide any accurate proof of this result for non constant di¤usion coe¢ cients. We give an accurate demonstration of (109) for arbitrary dimensions and di¤usion coe¢ cients in the Appendix of the present paper. Note that the case of a non-constant viscosity is important for modelling multiphase ‡ows; it can also be required when considering turbulent ‡ows modelled through Boussinesq's eddy viscosity assumption as a closure for the Reynolds stress tensor (see e.g. Pope, 2000 or Violeau and[START_REF] Violeau | Numerical modelling of complex turbulent free surface ‡ows with the SPH Lagrangian method: An overview[END_REF] for its application to SPH).

All second order operators can be deduced by applying ( 109) to an arbitrary component A l of a vector, then multiplying by the unit vector corresponding to the m th axis:

@ @x i K A @A l @x j a e m X b V b K A;ab A ab;l r ab [(n + 2) e ab;i e ab;j ij ] _ w ab e m (110) Contracting indices l = m and i = j gives [div (K A grad A)] a 2 X b V b K A;ab A ab r ab _ w ab (111)
which is identical to (107). Next, coming back to (110) and setting l = i and m = j yields

h div K A (grad A) T i a X b V b K A;ab r ab [(n + 2) (A ab e ab ) e ab A ab ] _ w ab (112) Then, setting l = j and m = i in (110) leads to [grad (K A div A)] a X b V b K A;ab r ab [(n + 2) (A ab e ab ) e ab A ab ] _ w ab (113)
(note that the discrete forms ( 112) and ( 113) are identical, although the corresponding continuous operators are equal only if K A is a constant). The last three approximations allow to write an SPH form of shear forces (101) as

(div ) a X b V b r ab (n + 2) ab 3 + ab (u ab e ab ) e ab + 5 ab 3 ab u ab _ w ab (114)
For n = 2 and constant viscosities, the latter equation gives the viscous force proposed by [START_REF] Espanol | Smoothed Dissipative Particle Dynamics[END_REF]. The case of an incompressible ‡ow (div u = 0, or = 2 3 according to (101)) simpli…es to:

(div ) a X b V b ab r ab [(n + 2) (u ab e ab ) e ab + u ab ] _ w ab (115)
If in addition to the incompressibility assuption, we state that the viscosity is constant in space, we observe that div h (grad u)

T i = grad div u = 0 (116)
Thus, according to (111) we get

(div ) a = [div (grad u)] a (117) 2 X b V b u ab r ab _ w ab
This form of SPH viscous forces was proposed …rst by [START_REF] Morris | Modelling low Reynolds number incompressible ‡ows using SPH[END_REF] with a non-constant viscosity, similarly to (107). Finally, under the same asumptions we may invoke (116) to write

(div ) a = h div grad u + 2 (grad u) T i a (118) 2 (n + 2) X b V b
u ab e ab r ab _ w ab e ab which agrees with the standard SPH arti…cial viscous forces proposed by [START_REF] Monaghan | Smoothed Particle Hydrodynamics[END_REF], as well as the molecular of turbulent viscous forces suggested on identical backgrounds by [START_REF] Violeau | Numerical modelling of complex turbulent free surface ‡ows with the SPH Lagrangian method: An overview[END_REF], although both papers uncorrectly assumed this model to be valid for variable viscosities. We thus got four formulae ((114), ( 115), ( 117) and ( 118)) to estimate shear forces, the last three being applicable to incompressible ‡ows and the last two to constant viscosities only. Refering to section 2.2, we notice that all of them consist of individual friction forces F diss b!a between pairs of particles, according to the general form ( 28) and (32). Equation ( 114

Equation ( 117) gives

ab = 2 V a V b r ab _ w ab I n (121)
and ( 118) yields

ab = 2 (n + 2) V a V b r ab _ w ab e ab e ab (122)
(note that, with a constant , (120) is the average of ( 121) and ( 122)). In formulae ( 119) to (122) the dependency of the ab 's on particle distance r ab stands through the quantities _ w ab = _ w h (r ab ). The symmetry condition (30) of the matrices ab (i.e. the action-reaction law (36)) is satis…ed by these four models, but not the de…nite positiveness. Indeed, the general form (119) yields the following dissipation function from (34), after a few algebra:

F = 1 2 X a;b V a V b r ab 2 4 ab 3 h (n + 2) (e ab u ab ) 2 + 5u 2 ab i + ab h (n + 2) (e ab u ab ) 2 u 2 ab i 3 5 _ w ab ( 123
)
whose sign is not …xed, due to the factor u 2 ab . In contrast, note that the original Navier-Stokes equation ( 100) gives the energy lost by dissipation (for a closed ‡uid, ideally without contact forces at its boundary) as

dH dt = Z : s d (124) = Z h s D : s D + (div u) 2 i d
Equation ( 124) gives the energy dissipation from a positive form (provided and are positive, see [START_REF] Landau | Fluid Mechanics[END_REF], which can be seen as an extension of (42) to continuous ‡uids. Equation ( 123) is thus a discrete approximation of (124), and its non-positiveness points out a weakness of the proposed model ( 109) for the SPH second-order derivatives. However, for incompressible ‡ows, the speci…c forms (120) to (122) respectively read

F = 1 2 X a;b V a V b ab r ab h (n + 2) (e ab u ab ) 2 + u 2 ab i _ w ab (125)
for (120); then

F = X a;b V a V b r ab u 2 ab _ w ab (126)
for ( 121), and …nally

F = (n + 2) X a;b V a V b r ab (e ab u ab ) 2 _ w ab (127)
for (122). Their positiveness is ensured provided _ w ab < 0 for b 6 = a, i.e. if the following condition is satis…ed:

r > 0 : _ w h (r) < 0 (128)
Thus the kernel derivative must be negative for all non-zero particle distances, a condition satis…ed by the kernel (61) introduced above (the behaviour of _ w h (r) in the vicinity of the origin will be stated below). On similar (thermodynamical) ideas, it has been noticed by [START_REF] Monaghan | Smoothed Particle Hydrodynamics[END_REF] that this condition is necessary to model scalar (e.g. temperature) ‡uxes consistent with their gradients. The de…nite positiveness of energy dissipation through dissipative matrices (120) to (122) are then ensured; however, they possess di¤erent properties in view of the considerations made in section 2.2. The proposition (117) presents the advantage to give shear forces collinear to the velocity di¤erence u ab , in accordance to what the theory of continuous incompressible ‡uids suggests with = s. Nevertheless, our considerations regarding discrete dissipative forces show that the forces given by this model do not conserve angular momentum nor vanish for a rigid body motion. The latter two conditions require ab = ab r ab r ab (equation ( 48)), which is only satis…ed by (122) (i.e. ( 118)) with

ab = 2 (n + 2) V a V b r 3 ab _ w ab > 0 (129)
which thus seems to be the best of the three proposed formulae, but again, only experimental tests will drive us to …rm conclusions about that point. In this case, the dissipation function ( 127) can be recovered from (50). Moreover, equation ( 51) tells us that ab (r ab ) r 2 ab must tend to a positive, non-zero quantity g a when r ab tends to zero. With (129), we can conclude that _ w h (r) must be equivalent to r when r tends to zero, in other words, using the dimensionless kernel (57), one must have _

f (q) q !0 C w q (130)
where C w is a positive, non zero constant depending on the kernel. As a consequence, at the origin (r = 0) the kernel …rst derivative must be null and its second derivative non-zero, negative:

_ f (0) = 0 (131) • f (0) < 0
The kernel (61) satis…es these conditions with C w = 120. Dimensionally speaking, with (67), (130) reads

_ w (r) r !0 C w w;n h n+2 r (132)
A kernel which would not satisfy the conditions (131) would yield a wrong estimation of volumes in the phase space, and thus an incorrect behaviour of the particle motion, as we will see later.

With model (118) for discrete viscous forces, the momentum equation (82) now takes the form 8a;

du a dt = 1 a X b V b " 2k b p a + 2k a p b ( a b ) k + 2 (n + 2) u ab e ab r ab # _ w ab e ab + g (133)
The three viscous models ( 115), ( 117) and ( 118) were tested on the case of the hill ‡ow presented above (section 3.1). It appears that, due to the incompressibility of the ‡ow and the constantness of the viscosity , they all give identical velocity pro…les (not plotted here, since the curves are almost superimposed). Hence, despite the fact that the …rst two models do not conserve angular momentum, the velocity distributions are not a¤ected here. So far, we used in all our computations the kernel (61), which satis…es the conditions (131). In a second step, we tested two kernels violating those rules. First, we considered the following kernel:

f (q) = 1 2 q 2 + 2q 2 (134) _ f (q) = q + 2
having a non-zero …rst derivative at the origin. Its radius is R f = 2, while (60) gives w;2 = 3 28 (135) in two dimensions. Figure 5 shows that the velocity pro…les, even timeaveraged, are very scattered and underestimated. The kernel used here does not ful…ll the rule (51), so that (55) shows that it preserves the volumes Ṽ in the phase space instead of reducing them, which is the reason of the observed instability. Next, we considered the following kernel:

f (q) = 1 4 q 4 4 3 q 3 + 3 2 q 2 + 9 4 (136)
_ f (q) = q (q 1) (q 3) which has a zero …rst-order derivative at the origin, but a positive secondorder derivative for small q's. Its radius is R f = 3, while (60) gives, in two dimensions:

w;2 = 20 243 (137)
This kernel gives better results than (134) in our test case when looking at time-averaged velocities (not plotted here), but …gure 6 reveals that the instantaneous velocity pro…les are again scattered, while with the kernel given by (61) they keep rather smooth and stable while time runs.

To …nish, note that from (129), one can …nd the quantities g a introduced at the end of section 2.2:

g a = lim r ab !0 ab r 2 ab (138) = 2 (n + 2) C w w;n h n+2 V 2 a
Thus, according to (55), the relative decrease of phase space volumes is given

1 Ṽ d Ṽ dt = 2 (n + 2) C w w;n h n+2 X a V a a (139)
As expected, it is proportional to the ‡uid viscosity. One should not be surprised of the fact that equation ( 139) depends on the smoothing length h (i.e. on the spatial discretization) and on the choice of the kernel (through C w and w;n ). The considerations at the end of section 2.2 show that the process of volume decrease in the phase space is strongly dependent on the scale at which the system is observed. As a matter of fact, dimensional considerations show that there cannot exist any formula giving the phase volume decrease in the continuous formalism of ‡uids, such as (124) for the energy dissipation.

Conclusions

We have seen that the variational viewpoint can easily be applied to Lagrangian numerical methods such as SPH. This approach leads to consistent equations, which should be carefully chosen to make the estimation of density and the momentum equation consistent to each other. Nevertheless, there is numerical evidence that in practice, the SPH method requires less precautions in writting the equations. One of the conclusions of our tests, however, is that the estimation of density bene…ts from using a continuity equation, while using a density interpolation is likely to yield instabilities, even when satisfying the abovementioned consistency conditions. Furthermore, we have given a solid demonstration of the general forms of second-order derivatives in SPH. Among those forms, our investigations on discrete particle dissipative forces reveal that only a limited number are consistent with the laws of Physics, in particular with conservation of angular momentum. However, again, SPH seems to accept non-conservative discrete viscous forces without any damage. On the other hand, a major requirement is that at the origin, the kernel w h (r) must have a …rst-order derivative which vanishes while its second-order derivative must be non-zero, negative in order to ensure stability. These constraints were found to be closely linked to the theoretical behaviour of a dissipative mechanical system in the phase space. In the past, other conditions have been found to be necessary to keep SPH simulations stable in various acceptations: for example, it is known that the second-order derivative of the kernel must be non-zero, positive for all non-zero r ab 's in order to avoid particle clustering (see e.g. [START_REF] Swegle | Smoothed Particle Hydrodynamics Stability Analysis[END_REF]. The latter condition, together with our conditions, is obviously not easily satis…ed. These considerations should be carefully considered by SPH developers, who often use arbitrary kernels. Besides, similar considerations could be extended, mutatis mutandis, to other particle methods like DPD and MPS.

Appendix

We give here a demonstration of formula (109): @ @x i K A @A @x j a where G is the following second order tensor:

G = grad [K A (grad A)] = @ @x i K A @A @x j e i e j (142)
(in which the e i 's are the unit vectors of a cartesian basis). We …rst write an integral approximation of (141), according to the SPH interpolation (56):

J = Z K A;ab A ab (r) r (n + 2) r r r 2 I n _ w h (r) d n r (143)
Next, we use a Taylor expansion of K A A in the vicinity of the particle a, similarly to (103) but up to the second order: On this basis, the integral (143) can be developed as

J = J 1 (K A grad A) a + 1 2 J 2 : G a + J 3 (147)
The integrals J 1 , J 2 and J 3 are tensors of respective orders 3, 4, and 2, de…ned by The property (59) shows that J 1 and J 3 are null; thus, (147) gives

J 1 = Z 1 r 2 (n + 2) 1 r (r r
J = 1 2 J 2 : G a + O h 4 (149)
To calculate J 2 , we write the volume of integration in dimension n as d n r = r n 1 drd!, where d! is the in…nitesimal angle of dimension n, which allows separating the variables:

J 2 = Z ht 0 _ w h (
Z R f 0 f (q) q n 1 dq = n S n
The last two integrals in (150) can be calculated from symmetry conditions. We begin with the second one: the components of the second-order tensor e e are e i e j , where e i = x i =r denotes the i th component of the unit vector e (not to be confused with the i th unit vector of the basis e i , written in bold). Obviously, changing e i into e i keeps the integral unchanged, since the sphere is invariant through plane re ‡ections. Thus, the only integrals of e i e j which do not vanish are those with i = j. Besides, they are all equal to each other for reasons of isotropy, while their sum is the integral of e i e i = 1. From this it immediately follows that Z n e e d! = S n n I n (152)

For the integral of e e e e in (150), similar considerations show that the only non-zero components are of the form e 4

i or e 2 i e 2 j with i 6 = j. Observing that a rotation of angle =4 around the origin in the (i; j) plane keeps the sphere unchanged, we get Z

  ac e ac + m a gIn the second summation, only the term involving b = a survives, to give X

V

  b w ab at time n + 1 to give (n+1) a .

  [(n + 2) e ab e ab + I n ]

  ab A ab r ab [(n + 2) e ab;i e ab;j ij ]

K

  A;b A b = K A;a A a K A;a (grad A) a r ab A a (grad K A ) a r ab (K A ) a r ab ] [(grad A) a r ab ] + O r 3 abwhere the double contraction of two tensors A and B of respective orders p and q is de…ned by A : B = A i1i2:::ip 1ip B ipip 1i3:::iq 1iq , whith Einstein's summation convention regarding the dummy indices i p 1 and i p . We then use …rst order developments of K A and A analogous to (104): (grad A) a r ab = A ab + O r 2 ab (145) (grad K A ) a r ab = K A;ab + O r 2 ab Combining (144) and (145) gives K A;ab A ab (K A grad A) a r ab 1 2 G a : (r ab r ab ) + O r 3 ab (146)

  r) I n r _ w h (r) d n r n _ w h (r) d n r

  where e is the local unit vector r=r, running over the n sphere n of radius unity and surface S n . The …rst integral in (150) is easily calculated from the kernel normalization condition (60):Z ht 0 _ w h (r) r n dr = [w h (r) r n ]

  and Revenga(2003) and[START_REF] Monaghan | Smoothed Particle Hydrodynamics[END_REF]. The above results can be used to rewrite the expected integral in the very simple form Z n e e e e d! = S n n (n + 2) ( ij kl + ik jl + il jk ) e i e j e k e l (160) (the terms involving the integral of e 4 i are implicitely included into the other terms). Finally, (150), (152) and (160) lead to J 2 = ( ik jl + il jk ) e i e j e k e l (161)

Coming back to (149), we now get the expected result:

x n a e i e j (162) = @ @x i K A @A @x j a e i e j = [grad (K A grad A)] a