A FIRST ASSESSMENT OF SMOOTHED PARTICLE HYDRODYNAMICS
IN FLUID MECHANICS

D. Violeau and M. Benoit

Laboratoire National. d’Hydraulique et Environnement (EDF),

6 quai Watier 78400 CHATOU, France

Tel. +33 (0)1 30 87 78 31 - Fax. +33 (0)1 30 87 80 86 - damien. E-mail: violeau@edf.fr

Abstract: The Smoothed Particle Hydrodynamics approach, a gridless lagrangian numerical method, has been used to develop a model capable to reproduce with good accuracy a large variety of flows. This model has been tested on various test cases, including a dam break and wind wave propagation, showing results in good agreement with experiments and theory. Some improvements have now to be explored in order to take into account the effects of turbulence.

Keywords: fluid motion, numerical models, lagrangien methods, smoothed particle hydronamics.

1    INTRODUCTION

Among the flows in which industry takes interest, some are quite difficult to reproduce using classical numerical methods. It is particularly visible when considering a flow including free surfaces and head losses (flow inside an outlet or a spillway, for example), or complicated free surface motions, such as wave breaking. In the last decades, however, new numerical approaches have been developed in this purpose. Smoothed Particle Hydrodynamics (SPH) is one of them. The SPH method is a fully lagrangian numerical technique which was first invented to simulate astrophysical phenomena, by Lucy (1977) and Gingold & Monaghan (1977). This method has proven its ability to reproduce fluid motion (Monaghan 1994). The Laboratoire National d’Hydraulique et Environnement (LNHE) of Electricité de France (EDF) has developed a model based on SPH to simulate free surface flows, which has proved to be capable to reproduce many flows with good accuracy (Violeau et al. 1999).

2    THE SPH GOVERNING EQUATIONS

SPH is a particle method, which does not require a grid to calculate spatial derivatives. Instead, we use an analytical differentiation of interpolation formulae based on a kernel (interpolating) function. The equations of motion and continuity become ordinary differential equations where the body forces are expressed as forces between particles (Monaghan 1992). To build the SPH formalism, we define fluid particles, considered as small bodies. Each particle a, located at r_a, carries a constant mass m_a, and variable parameters: the density r_a, the pressure p_a and the velocity u_a (and other physical properties if necessary). By using a kernel, the equation of motion for one fluid particle a can be written as:

                         (1)

in which the summation runs over all the particles b inside the domain. F_a is a body force (here gravity), and Ñ_a denotes the gradient with respect to the co-ordinates of the particle a. P_ab is a viscous term, depending on the relative locations and velocities of the particles :

                                  (2)

where x depends on the viscosity of the fluid. The kernel is a smooth, regular function of compact support h, called the smoothing length. As a consequence, when a particle b is far enough from a, its contribution to the motion of a vanishes. The very interesting thing is that the equation of motion (1) is not only a numerical approximation of the Navier-Stokes equations, but can also be derived from an action principle (Monaghan 1995, Violeau 2000). As a consequence, the total linear momentum of the fluid is conserved, as it is clear from (1), looking at the symmetry of the right hand side, representing the force exerted on particle a. It can be shown, under the assumption of a spherical kernel, that total angular momentum is also conserved (Violeau 2000). A suitable SPH form for the continuity equation is given by :

                        (3)

Equation (1) is used to calculate the velocities and positions of the particles at each time step. By the same way, equation (3) gives the evolution of density at the location of each particle. To calculate the pressure, we use an equation of state, such as following, valid for water :

                               (4)

where r₀ is a reference density and c₀ the speed of sound. However, for numerical consistency, it is necessary to chose a time step inversely proportional to c₀ (Monaghan 1992). When using the actual speed of sound, the time step is two small for practical applications. The usual method is to prescribe an artificial speed of sound, ten times higher than the maximum fluid velocity. Since the relative density fluctuations of a compressible fluid are of the order of magnitude of M² (M is the Mach number), here it does not exceed 1 %. Equation (4) is thus suitable for a weakly compressible flow (Monaghan 1994, Cummins et al. 1997).

Boundary conditions are not difficult to prescribe, except for solid boundaries. A usual method is to introduce boundary particles, which exert repulsive forces on fluid particles (Monaghan 1994, Morris et al. 1997). These forces are included into the term F_a of equation (1). One advantage of this approach is that these boundaries can easily move during the computation.

3    COLLAPSE OF A WATER COLUMN IN A TANK

SPH is robust and quite easy to program. The method is able to reproduce flows which cannot be simulated with classical numerical methods, including large free surface motions, like plunging waves. LNHE model is a 2DV model based on the SPH formulation, using the method and equations previously described. Since only near neighbour particles interact with one given fluid particle, link lists are built at each time step in order to reduce computational time and required memory.

The model has been tested on several test cases, particularly a schematic dam break (collapse of a water column in a tank). The tank length is 2 m, particles being located on the left side at the beginning of the computation, in a square of size H = 1 m (see Fig. 1). This test case contains 10000 particles of water, falling into the tank and making a plunging wave. The time evolution of the dimensionless maximum x-position of water (x / H) and the water depth on the left side boundary (h / H) is plotted on Fig. 2, versus dimensionless time. They show a very good agreement with experimental data (Monaghan 1994).

Fig. 1    Collapse of a water column with SPH (10 000 particles). The units are in meters.

Fig. 2    Collapse of a water column with SPH. Maximum x-position and the water depth on the left boundary with experiments: comparison of the computed results with experiments.

4    GRAVITY WAVE PROPAGATION

The propagation of gravity waves in a flume has also been simulated. The number of fluid particles is about 4000, and the flume length is 70 m. On the left side, a moving wave paddle creates the wave signal with given amplitude and period. The model correctly reproduces the propagation of waves, including wave breaking (Fig. 3). With a given amplitude of paddle, the computed motion gives correct values for the wavelength L, according to Stokes first order wave theory, as shown in table 1.

Fig. 3    Wave propagation in a flume using SPH.

Fig. 4    Wave modelling: velocity profiles compared to Stokes’ theory.

An example of maximum horizontal velocity profiles on a vertical section (between bottom and free surface) is shown on Fig. 4, compared to theoretical (Stokes) profiles. We see that the computation is again in good agreement with theory. Nevertheless, some discrepancies appear. They can be explained by the fact that to date, the effects of turbulence are not taken into account. This does not mean that SPH is a laminar model, but large-scale diffusion is considered through a constant bulk viscosity included in the diffusion term (2). As a consequence, large turbulent eddies can develop in the flow, while small eddies cannot because particles are considered as rigid bodies (Cleary & Monaghan 1999b).

5    FURTHER WORK: TURBULENCE MODELLING

To avoid the problem previously mentioned, a turbulence model should be developed. Since SPH is based on particle motion, it seems convenient to ponder a method based on LES. This approach should be particularly suitable since convolution is at the origin of the SPH equations. The size of eddies that can be represented is related to the average particle separation. In classical LES, structures smaller than the mesh size are represented by a turbulent viscosity function of the cell size. In SPH, such a turbulent viscosity, can be introduced. In addition, the diffusion term given by equation (2) may lead to questions of efficiency (Takeda et al. 1994, Cleary & Monaghan 1999a). Other numerical approaches should be examined for the diffusion term, such as stochastic methods (Welton & Pope 1997), which consist in adding a random walk to the standard convection of SPH. Laurence (1995) conducted an LES of grid turbulence with such a stochastic method, with a random velocity constructed on the Smagorinsky model and obtained the same spectra as with the standard diffusion model. An advantage is that the method is much more general than the standard eddy-viscostity subgrid model. Indeed, Durbin & Speziale (1993) show that a stochastic process can easily reproduce the effect of a second moment model, thus introducing anisotropy effects.

6    CONCLUSION

A 2DV numerical model, based on the SPH method, has been developed and tested. It has proved to be of a great interest, because of its simplicity. In addition, SPH reproduces correctly the motion of fluids, in accordance with theory and measurements. However, no turbulence model has been implemented yet. Further investigation will focus on this point, considering recent improvements in this field.