Mathematical Physics - Volume II - Numerical Methods

Chapter 5. Review of Development of the Smooth Particle Hydrodynamics (SPH) Method

162

Benz [16] and Monaghan [58] cover the early development of SPH. Right from the early days of SPH the importance of the smoothing kernel function as the essential feature of the SPH scheme was recognised. The Gaussian and the cubic B spline kernel functions are the most widely used in SPH, see Monaghan and Lattanzio [60]. However, most practical work relies on the monotone splines which, when used with small supports, allow for more accurate numerical solutions and higher numerical efficiency according to Balsara [7]. Liu et al. [50] among a number of other researchers demonstrated that in general, regardless of the choice of kernel function, the SPH method is not even zero order consistent. This is a consequence of the fact that the accuracy of the kernel interpolation depends on the distribution of the interpolation points within the kernel support. This effect is especially pronounced in the vicinity of boundaries, where the kernel support extends beyond the domain of the problem considered and consequently becomes incomplete Liu et al. [52]. Libersky and Petchek [47] extended SPH to work with the full stress tensor in 2D. This addition allowed SPH to be used in problems where material strength is important. The development SPH with strength of materials continued with its extension to 3D by Libersky [48]. Applications of SPH to model solids, i.e. material with strength, further highlighted shortcomings in the basic method: consistency, tensile instability, zero energy modes, treatment of contact and artificial viscosity. These shortcomings were discussed in detail in the first comprehensive analyses of the SPH method by Swegle [87], Wen [105]. The problems of consistency and accuracy of the SPH method, identified by Belytschko [11], were addressed by Randles and Libersky [74], Vignjevic and Campbell [97] and a number of other researchers. This resulted in a normalised first order consistent version of the SPH method with improved accuracy. The attempts to ensure first order consistency in SPH resulted in emergence of a number of variants of the SPH method, such as Element Free Galerkin Method (EFGM) Belytschko [13], Kongauz [44], Reproducing Kernel Particle Method (RKPM) Liu [52][54], Moving Least Square Particle Hydrodynamics (MLSPH) Dilts [22], Meshless Local Petrov Galerkin Method (MLPG) Atluri and Zhu [1]. These methods allow the restoration of consistency of any order by means of a correction function. It has been shown by Dilts [22] and independently by Atluri et al. [5] that the approximations based on corrected kernels are identical to moving least square approximations. A comprehensive stability analysis of particle methods in general by Belytschko [14], Xiao and Belytschko [106], and independently by Randles [74] who worked specifically on the SPH method provided improved understanding of the methods analysed and confirmed the conclusions from Swegle’s initial study. Randles’ unique analysis, which included space and time discretisation, showed that SPH can be stabilised by precise choice of time step size and predictor corrector type of time integration. Rabezuk et al. [73] demonstrated that if used within a total Lagrangian framework SPH does not exhibit the tensile instability. Tensile instability in SPH has, as a consequence, non-physical motion of particles which form clusters. This was first observed in materials loaded in tension (negative stress), however the instability can develop under compressive loading, see Swegle [87]. In simulations of solids the instability may result in non-physical numerical fragmentation. Another unconventional solution to the SPH tensile instability and zero energy mode instability problems was first proposed by Dyka [23] in which the stresses are calculated at