Simulating Flows with SPH: Recent Developments and Applications

Fluid motion is often numerically reproduced by means of grid-based methods such as Finite Difference Methods (FDMs) and Finite Elements Methods (FEMs). However, these techniques exhibit difficulties, mainly related to the presence of time-dependent boundaries, large domain deformations or mesh generation. This chapter describes a relatively recent meshfree and pure Lagrangian technique, the Smoothed Particle Hydrodynamics (SPH) method, which overcomes the above mentioned limitations. Its original frame has been developed in 1977 by (Gingold and Monaghan, 1977) and independently by (Lucy, 1977) for astrophysical applications. Since then, a number of modifications to ensure completeness and accuracy have been yielded, in order to solve the main drawbacks of the primitive form of the method. Because a calculation is based on short-range particle interactions, it is essential to limit the computational costs related to the neighbourhood definition. Available searching algorithms are then presented and discussed. Finally, some practical applications are presented, primarily concerning free surface flows. The capability to easily handle large deformations is shown.

. Dirac delta function centered at the point x (for one dimensional problems).
where  represents the domain of definition of f and x,y . Replacing  with a smoothing function W(xy, h), eq. (1) can be approximated as in which W is the so called smoothing kernel function or simply kernel and h, acting as spatial scale, is the smoothing length defining the influence area where W is not zero. While eq. (1) yields an exact formulation for the function f(x ), eq. (2) is only an approximation. The definition of W is a key point in the SPH method since it establishes the accuracy of the approximating function f(x ) as well as the efficiency of the calculation. Note that the kernel approximation operator is marked by the index I. The kernel function W has to satisfy several properties (Monaghan, 1988;Vila, 1999). The following condition W xy , h d y is known as partition of unity (or the zero-order consistency) as the integration of the smoothing function must yield the unity. Since W mimics the delta function, eq. (3) can be rewritten as a limit condition in which the smoothing length tends to zero Still, W has to be defined even, positive and radial symmetric on the compact support ( Fig. 2) where  is a positive quantity. A large number of kernel functions are examined in literature, e.g. quadratic to quintic polynomials, Gaussian etc. Among the others (Liu & Liu, 2003), a smoothing function satisfying the above condition is the cubic spline based kernel (Monaghan & Lattanzio, 1985) defined as: www.intechopen.com W = A n d 2 3 -q 2 + q 3 3 0 ≤ q < 1 W = A n d 2-q 3 1 ≤ q < 2 W = 0 otherwise (6) where  = 2, A(n d ), depending on the number of dimensions n d , denotes a scaling factor that ensures the consistency of eq. (3) whereas q denotes the dimensionless distance xy h ⁄ .

Fig. 2. Typical shape of the smoothing function W
Knowing the function f carried by a collection of moving particles, the integral representation given by eq. (2) can be converted into a discretized summation over all particle N within the compact support ( Fig. 2), yielding the particle approximation where the index k refers to particles within the compact support (see bold ones in Fig. 2), with mass m k and density  k being carried. Note that in this case the particle approximation is marked by the index a. The subscript will be avoided from now on. The presence of these two variables allows SPH to be easily and conveniently applied to hydrodynamics problems. The smooth estimate eq. (7) can be referred to a generic particle occupying the position x i , as follow . Particle approximation of spatial derivatives of a field function, such as divergence and gradient, is determined using the gradient of the kernel function rather than from the derivatives of the function itself where the nabla operator is referred to the location of particle I. The symbol "·" denotes the dot product. Eqs. (9)-(10) offer the great advantage of estimating their left hand side in term of the kernel gradient, i.e. allowing no special hypotheses on the particular field function (vector in the first case, scalar in the latter). A different formulation of the spatial derivative in eq. (9) can be achieved introducing the following identity (Liu and Liu, 2003)  inside the integral in eq. (2), yielding in this case Likewise the divergence, another particle approximation of the gradient can be derived, taking into account the following equation Since the field function on the right hand side of Eqs. (12) and (14) appears in the form of paired particles, such equations are conveniently employed in fluid dynamics, allowing the conservation of linear and angular momentum.

SPH-form of governing equations
The mostly used governing laws ruling fluid motion are the Navier-Stokes equations, which specify that mass and linear momentum (also expressed in Newton's second law which state that the rate of change of the momentum of a particle is proportional to the resultant force acting on the particle) are preserved. Conservation laws in Lagrangian form are in the following provided www.intechopen.com in which  and v are respectively the density and velocity field, p is the isotropic pressure,  the laminar kinematic viscosity and f the external force. Different approaches (Monaghan, 1994;Liu and Liu, 2003;Oger et al., 2007) are available in order to derive the density particle approximation of the continuity equation (15.a) and momentum equation (15.b). For instance, referring to eq. (12), the density rate at particle i can be approximated as follows The SPH formulation of the velocity variation can be deduced from eq. (14) An artificial viscosity term (Monaghan & Gingold, 1983) can be added in eq. (14), allowing shock waves to be properly simulated (Liu et al., 2003a), damping spurious oscillation near the fronts The dissipative term above introduced is the most general viscosity used in SPH computations, since it provides good results when modelling shock fronts and prevents interparticle penetration ( term in eq. (19)). It is defined as The notation a = a i +a k 2 ⁄ , b =b i -b k is introduced above. Term c refers to the speed of sound which magnitude has conveniently to be at least ten times greater than the maximum estimate of the scalar velocity field (Monaghan, 1994),  = 0.1h is employed to prevent numerical divergences when two particle are approaching.  typically ranges from 0 to 0.01, depending on how viscous is the fluid. When other approaches for reproducing shear stresses are taking in, e.g. (Cleary, 1998), then it is not necessary to employ the first term inside eq. (19). Problem closure is achieved by combining conservation equations in discrete form (16), (18) with an equation of state, when the weakly compressible scheme is adopted (Lee et al., 2010). A relationship between pressure and density is given by (Dymond & Malhotra, 1988) where c 0 is the reference speed of sound, large enough to guarantee Mach numbers lower than 0.1 0.01, γ = 7, ρ 0 = 1000kg/m 3 , if the liquid is water.

Completeness and accuracy of the SPH method
In SPH the consistency order of a function (or more properly completeness or reproducing conditions) means the order of a polynomial required to evaluate the function exactly, when using values carried by particles. Next section 3.1 deals with how restoring consistency to jth order, while section 3.2 provides some aspects related to the spatial scale h, which has a strong influence on the accuracy.

J-th order consistency
Zero-order consistency for particle approximation is attained when the smoothing kernel satisfies eq. (3) in discrete form Eqs. (22), (23) are satisfied when a corrected smoothing function, e.g. the Shepard function (Bonet et al., 2004) is adopted. (23) are not satisfied in the case of the evaluation near the boundaries, where kernel truncation yields a lack of interpolation or when particles are exceedingly disordered. There are several ways to restore higher order consistency conditions. j-th order consistency (j ≥ 1) is achieved with updated versions of the original SPH method, e.g. Reproducing Kernel Particle Method (Chen et al., 2000;Liu et. al., 1995), Element Free Galerkin Mehod (Belytschko et. al., 1994;Krongauz & Belytschko 1997), Moving Least Square Particle Hydrodynamics (Dilts, 1999). An interesting method to ensure j-th order consistency has been proposed by (Liu et al., 2003b). Using the Taylor series expansion for the kernel and introducing the above equation in the following set, representing the reproducing condition of a function to the j-th accuracy in integral form (Liu and Liu, 2003)   a discretized consistency condition to j-th order is ensured by the resulting set eqs. (27) can be rewritten in a matrix fashion, such as which solution yields the j+1 coefficients G x , h . While the restoring condition to j-th order is ensured by eqs. (27), still some open issues arise. First of all, the corrected kernel may not be symmetric, in contrast with the requirements stated in eq. (5.a). Second, it might be negative somewhere, which may produce unphysical field results (e.g. non positive density). Finally, it might show secondary peaks aside the central (see Fig. 2), hence affecting the internal force between two approaching particles. Consequently, special care has to be taken when simulating flows with corrected version of SPH.
Concerning the stability of SPH, a pioneering work is given by (Swegle et al., 1994), whose results have then been confirmed by (Randles et al., 1999) and independently by (Belytschko, & Xiao, 2002).

Variable smoothing length
The introduced smoothing length h is crucial for the accuracy of the particle approximation. When h is too large, carried particle properties are excessively smoothed out in space, affecting the accuracy. Vice versa, too small values determine a little number of computing particles inside the compact support, determining a lack of interpolation. In SPH computations, h may be kept fixed or variable in time and space. Since the carried particle's mass is usually constant, the number of neighbouring particles inside the support domain should not vary, according to the following where h 0 and  0 are respectively the reference smoothing length and density. When h is chosen to be variable, special concern must be taken with reference to the third law of Newton (action -reaction principle). Indeed, if computing particles a and b have different h a and h b , then their compact supports are not the same, causing different exerting forces between such pairs. Symmetry of particle interactions must then be preserved. This condition can be achieved, taking into account a symmetrised smoothing length, as proposed by (Benz, 1989) Others approaches exist to ensure symmetry (e.g. the geometric average instead of the arithmetic average).

Neighbouring search methods
A limitation of SPH is that a single computation becomes very demanding in terms of running time as the number of particle increases. For each particle, since the integration of the governing equations is carried out on a limited number of adjacent particles located inside a cut-off distance r c = h, a large part of the computational burden depends on the actual searching procedure of the neighbouring particles. It is therefore crucial that efficient methods are adopted for such a search. The cut-off radius is indeed much lower than the typical domain's spatial dimension, hence the number of neighbouring particles N is a little fraction of the total number N tot . Straightforward determination of which particles are inside the interaction range would require the calculation of all pair-wise distances, a procedure whose computational burden would be of the order O(N tot 2 ), and therefore unpractical or totally impossible for large problems.
Two main approaches have been developed in the past in order to reduce the unnecessary computation of distances: the first based on dynamically storing each particle's neighborhood list, namely the Verlet list (Verlet L, 1967), and the second based on a framework of fixed cells (Allen and Tildesley, 2000). The first dataset gathers the N V,i , i = 1,…,N, neighbours contained in a range r v slightly greater than the cut off distance r c (Fig. 3). Fig. 3. Verlet list definition for the fluid particle "i"; particles marked with a cross between r c and r v , are included in the list but they do not give any contribution to the hydrodynamics properties of the targered particle.
The condition N V,i << N tot , will still hold, as long as r v is small enough, even thoughobviously -the number of neighbouring particles stored in the list is greater than those strictly required. The difference r v -r c = V skin represents a sort of safety "skin" around the cut-off distance r c and contains some elements which are initially unnecessary for the definition of the i th particle's properties, as shown by particles marked with a cross in Fig. 3. This takes into account the possibility that some or all of them will cross the interaction sphere and thus become part of the neighbourhood in the following time steps. With this approach the list is kept unchanged for some time steps, till a "refresh condition" is established. A number of criteria have evolved over the years as a condition to activate the list refresh operation, mostly based on the maximum possible distance traveled by the particle (Chialvo & Debenedetti, 1983;Blink and Hoover, 1985). The latter approach sets a partition of the physical space into fixed cells so that, for each particle, the neighbouring search only has to be performed within the surrounding particles, i.e. those inside the cell where the current particle is located and in the 8 (for 2D problems) or 26 (3D) adjoining ones. In order to define the neighbourhood, a data structure must be created so that each moving particle is connected to the grid cell it is located in; such a structure (linked-cell list) must be renewed at each time step, a relatively straightforward operation of O(N tot ) complexity. For this approach, as for the Verlet list method, the overall efficiency depends very much on the choice of parameters, balancing the size of the cells against their total number. Large cells imply a longer neighborhood search, while a greater number requires a longer data structure reconstruction. (Viccione et al., 2008) carried a numerical sensitivity analysis on the efficiency of the two procedures as function of the key spatial parameters, that is the Verlet list skin and the cell size.
=r v /r c Relative partial and total computational time Domain partitioned in cells whose dimensions are: Dx,cell=0,25m; Dy,cell=0,25m www.intechopen.com It appears that while both the linked cell list and the Verlet list do relieve the computational time, the comparative advantage of the linked cell increases to the point where it is practically of no use whatsoever. Later on, (Domínguez et. al., 2010) proposed an innovative searching algorithm based on a dynamic updating of the Verlet list yielding more satisfying results in term of computational time and memory requirements.

Applications of the SPH method
Smoothed Particle Hydrodynamics has been applied to a number of cases involving free surfaces flows.

Slamming loads on a vertical structure
The case of a sudden fluid impact on a vertical wall (Peregrine, 2003) has been examinated on a geometrically simple set up. (Viccione et al., 2009) shown how such kind of phenomenon is strongly affected by fluid compressibility, especially during the first stages. A fluid mass, 0.50m high and 4.00m long, moving with an initial velocity v 0 = 10m/s is discretized into a collection of 20.000 particles whith an interparticle distance d 0 = 0.01m. The resulting mass is at a close distance to the vertical wall, so the impact process takes place after few timesteps (Fig.  6). Timestep is automatically adjusted to satisfy the Courant limit of stability. Fig. 6. Initial conditions with fluid particles (blue dots) approaching the wall (green dots) The following Fig. 7 shows the results in terms of pressure at different times.   The rising and the following evolution of high pressure values is clearly evident. The order of magnitude is about 10 6 Pa, as it would be expected according to the Jokowski formula p = ρ C 0 Δv, with v = v 0 = 10m/s. After about 1/100 seconds most of the Jokowsky like pressure peak, generated by the sudden impact with the surface, disappeared, following that, the pressure starts building up again at a slower rate.

Simulating triggering and evolution of debris-flows with SPH
The capability into simulating debris-flow initiation and following movement with the Smoothed Particle Hydrodynamics is here investigated. The available domain taken from an existing slope, has been discretized with a reference distance being d 0 =2.5m and particles forming triangles as equilateral as possible. A single layer of moving particles has been laid on the upper part of the slope (blue region in Fig. 8).
Triggering is here settled randomly, releasing a particle located in the upper part of a slope, while all the remaining ones are initially frozen. Motion is then related to the achievement of a pressure threshold p lim (Fig. 9). The resulting process is like a domino effect or a cascading failure. While some particles are moving, they may approach others initially still, to the point for which the relative distance yields a pressure greater than the threshold value. Once reached such point, those neighbouring particles, previously fixed, are then set free to move. Runout velocity is instead controlled by handling the shear stress  bed with the fixed bed. Fig. 8. Spatial discretization. Red circles represent the area where local triggering is imposed. Fig. 9. Neighbour particle destabilization. a) Particle "i" is approaching the neighbour particle "j". b) Despite the relative distance "|r ij |" is decreased, particle "j" is still fixed because p ij < p lim . c) Particle "j" is set free to move because the pressure "p ij " has reached the threshold value "p lim ". As can been seen from the above Figures 10 to 12, by varying the location of the triggering area and the limit pressure p lim , the condition of motion are quite different. More specifically, the mobilized area increases when the isotropic pressure p lim decreases.

Conclusion
Recent theoretical developments and practical applications of the Smoothed Particle Hydrodynamics (SPH) method have been discussed, with specific concern to liquids. The main advantage is the capability of simulating the computational domain with large deformations and high discontinuities, bearing no numerical diffusion because advection terms are directly evaluated. Recent achievements of SPH have been presented, concerning (1) numerical schemes for approximating Navier Stokes governing equations, (2) smoothing or kernel function properties needed to perform the function approximation to the Nth order, (3) restoring consistency of kernel and particle approximation, yielding the SPH approximation accuracy.
Also, computation aspects related to the neighbourhood definition have been discussed. Field variables, such as particle velocity or density, have been evaluated by smoothing interpolation of the corresponding values over the nearest neighbour particles located inside a cut-off radius "r c " The generation of a neighbour list at each time step takes a considerable portion of CPU time. Straightforward determination of which particles are inside the interaction range requires the computation of all pair-wise distances, a procedure whose computational time would be of the order O(N 2 ), and therefore unpractical for large domains.
Lastly, applications of SPH in fluid hydrodynamics concerning wave slamming and propagation of debris flows have been discussed. These phenomena -involving complex geometries and rapidly-varied free surfaces -are of great importance in science and technology.

Acknowledgment
The work has been equally shared among the authors. Special thanks to the C.U.G.Ri. (University Centre for the Prediction and Prevention of Great Hazards), center, for allowing all the computations here presented on the Opteron quad processor machine. The constant evolution of the calculation capacity of the modern computers implies in a permanent effort to adjust the existing numerical codes, or to create new codes following new points of view, aiming to adequately simulate fluid flows and the related transport of physical properties. Additionally, the continuous improving of laboratory devices and equipment, which allow to record and measure fluid flows with a higher degree of details, induces to elaborate specific experiments, in order to shed light in unsolved aspects of the phenomena related to these flows. This volume presents conclusions about different aspects of calculated and observed flows, discussing the tools used in the analyses. It contains eighteen chapters, organized in four sections: 1) Smoothed Spheres, 2) Models and Codes in Fluid Dynamics, 3) Complex Hydraulic Engineering Applications, 4) Hydrodynamics and Heat/Mass Transfer. The chapters present results directed to the optimization of the methods and tools of Hydrodynamics.