Mathematical Physics - Volume II - Numerical Methods

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

176

Using the NCSPH approximations the conservation equations assume the following form: D ρ I D t = ρ I nnbr ∑ J = 1 m J ρ J ( v J − v I ) · ˜ ∇ ˜ W IJ , (5.56) D v I D t = − nnbr ∑ J = 1 m J σ I ρ 2 I + σ J ρ 2 J · ˜ ∇ ˜ W IJ , (5.57) 5.9 Tensile Instability A Von Neumann stability analysis of the SPH method was conducted Swegle et al. [87] and Balsara [7] separately. This has revealed that the SPH method suffers from a tensile instability. This instability manifests itself as a clustering of the particles, which resembles fracture and fragmentation, but is in fact a numerical artefact, see Figure 5.3 Swegle concluded that the instability doesn’t result from the numerical time integration algorithm, but rather from an effective stress resulting from a non-physical negative modulus being produced by the interaction between the constitutive relation and the kernel interpolation. In other words the kernel interpolation used in spatial discretisation changes the nature of original partial differential equations. These changes in the effective stress amplify, rather than reduce, perturbations in the strain. D e D t = − nnbr ∑ J = 1 m J ( v J − v I ) · ˜ ∇ ˜ W IJ . (5.58)

Figure 5.3: Typical colocated Eulerian SPH behaviour under tension. Although the linear elastic model was used for this simulation (i.e. no fracture is included in the constitutive model), unphysical fracture of the 2 − D specimen occurs as a consequence of numerical instability in areas of high tensile stresses. From Swegle’s stability analysis it emerged that the criterion for stability was that: W ′′ σ > 0 (5.59) where W ′′ is the second derivative of W with respect to its argument and σ is the stress, negative in compression and positive in tension. A stability analysis leading to the stability condition (5.60) is given at the end of this section. This criterion states that instability can also occur in compression, not only in tension. Indeed, if the slope of the derivative of the kernel function is positive, the method is