The nonlocal, local and mixed forms of the SPH method

From its early days the SPH method has been criticised for its shortcomings namely tensile instability and consistency. Without thorough understanding of the method attempts were made to make the classical SPH method consistent and stable which resulted in the local and Total Lagrangian forms of SPH similar to the finite element method. In this paper we derived and analysed a consistent nonlocal SPH which has similarity with Bazant’s imbricate continuum. In addition, the paper provides comparison and discussion of different SPH forms including: Classical SPH, Nonlocal, Local and Mixed SPH. The partition of unity approach was used to define the following two mixed forms: Local–Nonlocal and Local–Classical SPH. These mixed forms were intended for modelling of physical processes characterised with local and nonlocal effects (local and nonlocal constitutive equations), e.g. progressive damage and failure. The stabilising effect of the Local form on the Classical SPH, which is inherently unstable (tensile instability), are also illustrated. The stability analysis, presented in appendices A and B, demonstrate stability of the continuous and discrete form of the nonlocal SPH based on Eulerian kernels for elastic continuum.


Introduction
The development of the Smoothed Particle Hydrodynamics (SPH) method follows a typical path of any novel idea in science where an author/group of authors publishes the basic ideas.Hesitatingly at first, little by little, other original contributions appear, until a certain body of knowledge is reached.At this point, overview articles are published, topical conferences are organised, and a first mention is made in textbooks, until specialised books and monographs are written.SPH has reached that stage now.
The SPH method was originally proposed by Lucy [1] and Gingold and Monaghan [2] to simulate a wide range of problems in astrophysics, e.g. the formation and evolution of proto-stars, galaxies and the three-dimensional motion of compressible fluids at different spatial scales.
The meshless nature of the SPH method led to its successful application to other fields.According to the comprehensive bibliographic survey by S. Li and W.Z. Liu [3], the SPH method has been employed to solve problems of both compressible flow [4], incompressible or weakly compressible flow [5][6][7], free-surface flows [8,9], multiphase flow [10][11][12][13][14][15], plasma-fluid motion [16,17], general relativistic hydrodynamics [18,19], heat conduction [20,21], nonlinear dynamics [22] and solid mechanics [23,24].In this paper our interest is primarily in solid mechanics.Here we refer only to few references selected to illustrate development and applications of the SPH method with no disrespect to the, now, large body of related literature.
It is important to recognise the role of SPH rEsearch and engineeRing International Community (SPHERIC, w ww.spheric-sph.org)the international organisation which has been bringing together the community of researchers and industrial users of SPH since 2005.SPHERIC activities include training, methods development through a range of grand challenges and applications to real world problems.
This paper revisits derivation of the Eulerian version of the SPH method and offers a number of new explanations of different forms of the discretised balance equations.These developments follow our initial observation of nonlocal properties of the SPH method in modelling damage induced material softening [1].Furthermore, a local-nonlocal version of the SPH (LNL-SPH) method is proposed based on partition of unity.In order to provide a basis for the LNL-SPH derivation we first review relevant theories on which the method is based.Parallels between Eringen's nonlocal continuum theory, Bazant's Imbricate continuum and Classical SPH are discussed in Section 2. Emphasis is placed on the SPH properties common with these theories.Derivation of LNL-SPH formulation is given in Section 3. Section 4 outlines stability analysis of the new formulation.Section 5 contains a number of numerical examples which illustrate performance of LNL-SPH in analysis of elastic solids.This includes tensile instability consideration and comparison for different SPH forms.The main conclusions are stated in Section 6. Appendix A provides stability assessment of the continuous nonlocal SPH approximation for elastic continuum resulting in the stability requirements for kernel functions.Appendix B -provides stability assessment of the discretised nonlocal SPH method for elastic continuum according to Swegle [25].

SPH parallels with nonlocal and imbricate continuum theories
In classical (local) solid mechanics material points are regarded as fundamental entities which form a threedimensional differentiable manifold representing a deformable body [26,27].These points have three translational degrees of freedom.The response of the material to displacement of its points (deformation) is characterised by development of internal forces which can be quantified by a symmetric stress tensor defined locally at each point.
However, the atomic theory of solids provides evidence of the existence of long-range, cohesive forces and their effect on the dispersion of elastic waves.The inclusion of long-range effects has motivated developments of nonlocal continuum mechanics.
The nonlocal theories are characterised by postulates for the whole body which include localisation residuals in the balance equations which are obtained from the global balance statements.Further, the related formulation of constitutive relations allows for dependence on the nonlocal material properties (nonlocal effects).The inclusion of nonlocal effects is typically achieved by combining local and global postulates: local laws of balance of mass, momentum and moment of momentum in conjunction with a global law for the balance of energy or a global Clausius-Duhem inequality, [28].Classical local continuum models comprise partial differential equations (PDEs) based on local information in an infinitesimal neighbourhood and are derived, from the assumption of smooth fields.
Nonlocal continuum models have the potential to provide alternatives or to be used in combination with the local PDE models in many applications including solid mechanics.While mathematical analysis and numerical solution of local PDEs are well established branches of mathematics, the development of rigorous theoretical and computational frameworks for nonlocal models requires further research [29].
Local continuum theories are characterised by local approximations, e.g. the strain ε and stress σ tensors at the macroscopic material scale are regarded as certain averages of the microscopic (actual) strains and stresses taken over a representative volume V (appropriate for the material; considered).In the case of smooth (differentiable) strain fields the constitutive model can be defined as a relation between ε (x) and σ (x) at the same location (material point) x.However, if the strain field is characterised by high gradients in V, statistical analysis shows [30,31] that the local approach is no longer adequate, and that the entire macroscopic (averaged) stress distribution should be related to the entire macroscopic (averaged) strain distribution.

Eringen's nonlocal continuum formulation
One of the nonlocal continuum theories used to solve elasticity problems was proposed by Eringen [32].According to this theory the stress-strain relation is where x is the position vector; ε and σ are the macroscopic stress and strain tensors; C is the fourth order material stiffness tensor, C is the material nonlocal stiffness tensor; α (s) is a given weighting function and Ω is the weighting function domain.In applications of this theory it is often assumed that C ( x, x ′ ) = C (x) α (s) , s = ⏐ ⏐ x ′ − x ⏐ ⏐ ; where α (s) is given symmetric function with the following properties: ∫ Ω α (s) ds = 1, α (s) = 0 ∀ s / ∈ Ω ;.The above equation can be restated and interpreted as σ (x) = C (x) : is the nonlocal strain tensor.Note that when the size of the kernel domain Ω tends to zero, the weighting function tends to the Dirac delta function, i.e. α and the nonlocal strain tends to the local strain ε (x) ⇒ ε (x).
For clarity reasons we limit the considerations in this paper to linear elasticity and small-strain formulations.Consequently, Eq. ( 2) can be restated in 1D as: where σ (x) is stress (macroscopic); E is elastic modulus; ε (x) is local strain ; ε (x) is weighted average (nonlocal) strain; H is averaging operator; and α (s) is a given weighting function, which is zero or negligibly small outside the region − l /2 ≤ x ≤ l /2, where l is a parameter which defines size of the representative volume.According to Eringen's nonlocal continuum theory the strain averaging defined by Eq. ( 3) is interpreted as the dependence of strain ε (x) at point x on the strains ε (x + s) at the locations x + s.This definition of nonlocal strain was questioned and an alternative definition offered by Bazant who stated that the above interpretation is inappropriate because the material at one point cannot "feel" the local strain at another point.The strain measure should be physically interpreted through its relation to the change of distance between points, e.g.x + s and x − s [33,34].The key aspects of the nonlocal approach proposed in those two papers is summarised and partially reproduced below, for 1D elastic continuum.This provides the background for application of a similar framework for interpretation and derivation of the SPH discrete forms of balance equations.

Bazant's imbricate continuum formulation for 1D
To illustrate effects of the strain averaging defined by Eq. ( 3) one can substitute ε (x + s) = du(x+s) ds (where u is displacement) into Eq.( 3) and integrate by parts.Using the property of the weighting function α (s) = 0 ∀ − l /2 ≥ s ≥ l /2, where l is a characteristic length, and noting that the symmetry property, α (s) = α (−s) requires anti-symmetry of the derivative dα(s) ds = − dα(−s) ds , nonlocal strain ε (x) becomes ε To simplify Eq. ( 4), the weights, α (s) are assumed to be uniform, i.e., α (s) = 1 /l.This assumption results in ds ds = 0 (to be discussed further in relation to SPH).Hence, the mean strain ε (x) defined by Eq. ( 4) becomes a finite difference expression Eq. ( 5) indicates that the nonlocal strain ε and stress σ depend on the relative displacement between the ends of segment l or, more precisely, that the stress resultant σ • A (force) over the cross section area A of the representative volume depends on the change of length of segment 1-2 in Fig. 1(a).Note, in the concept of imbricated elements points 1 and 2 are symmetrically located relative to cross section PQ (Is this a general requirement or specific to uniform α (s)).Similarly, a change of length of segment 3-4 in Fig. 1(a) determines the stress on cross section RS.
Following Bazant, we continue with consideration of a 1D continuum, i.e. a bar of a unit cross section.The bar is subdivided into N equal segments of length h by nodes k = 1, 2, 3, . . .,N + 1 as shown in Fig. 1(b).The bar is also divided into K overlapping segments of length l, see Fig. 1(c), subject to the condition that l = n • h in which n is an integer of 2 or greater and l is the given size of the representative volume.The mean strain ε, defined by the difference expression Eq. ( 5), is the strain obtained in one-dimensional elements of length l.These elements are shown as circles in Fig. 1(c).They span over n segments h, and are connected only at their ends, so that the strain in each of them is uniform.These elements overlap, i.e., are imbricated (see Fig. 1(c), the nodes located above each other have the same displacement).A cross section between two nodes, Fig. 1(c), is intersected by a total of n identical imbricate (regularly overlapping) elements acting in parallel.The combined cross-section area of all imbricate elements intersected by the same cross section of bar must be a constant, independent of n.Denoting this combined area as (1 − c) where 0 ≤ c ≤ 1, we see that the cross section of each imbricate element is equal (1−c) /n.Obviously, it decreases as the subdivision of l is refined and tends to 0 as n → ∞.
In their development [33] proceed by formulating the following hypothesis, for 1D continuum, in which l is a material property (a constant) called the material characteristic length.
Hypothesis I.The nonlocal stress σ at any point x (except in a boundary layer of thickness l/2): 1. Depends on the change of distance between points x + l/2 and x -l/2; but 2. does not depend on the change of distance between any other two points lying a finite distance apart.
It follows that the continuum must consist of all possible (continuously distributed) material elements connecting any two points lying at distance l from each other.Moreover, part 2 of Hypothesis I suggests that these elements should have a parallel and overlapping (imbricated) arrangement, in which the neighbouring elements are not mutually joined, see Fig. 1(c).A proof that such an imbricated arrangement unequivocally follows from Hypothesis I may be given by variational calculus.
The principle of virtual work was used to prove Hypothesis I.According to Hypothesis I the incremental virtual work, δℵ, depends on the local strain, ε, and on the nonlocal strain, ε.Therefore, for a bar of length L: where x is the length coordinate of the bar; u is displacement; δℵ b is the work done by stresses for segments of length l at each end; δℵ 0 , δℵ 1 is the works of nonlocal stresses σ and local stresses τ within the rest of the bar; p(x) is distributed load; δu(x), δε(x), δε(x) are kinematically admissible variations; and c is a coefficient with values between 0 and 1, defined above and illustrated in Fig. 1.
for ε, and ε = ∂u /∂x for ε, we have Introduction of new variables, ξ = x + l /2, η = x − l /2, and then renaming in separate integrals both ξ, η as x, yields After integrating by parts in δW 0 , and substituting p (x) = −ρ ∂ 2 u(x) ∂t 2 (ρ = mass density) according to d'Alembert principle, we obtain where F 1 and F 2 are certain functions which are zero for l /2 ≤ x ≤ L − l /2 and have effect only within end segments of length l /2.Since δℵ = 0 for any kinematically admissible δu(x), the expression in the braces { } of the first integral must be zero.This yields the one-dimensional continuum equation of motion, Eq. 19 of Bazant at.al. [33] (given below).
Eq. ( 19) in [33] ε Eq. ( 4) in [33] Note that Eq. ( 4), as well as Eq.19 of Bazant at.al. [33], applies only at points whose distances from the ends are at least l.The boundary conditions, which ensue from functions F 1 and F 2 , do not refer just to the end points, but are evidently spread out over boundary segments of length l (which can be generalised in two or three dimensions as boundary layers, i.e. volume constrained domains).Bazant proposes that the boundary conditions are applied to the finite element part of the discretisation [33].
The foregoing variational derivation shows that, if there is a difference operator (or gradient-averaging operator) in the strain-displacement relation, the same operator must be applied to stress in the continuum equation of motion (or equilibrium).The fact that this equation can be derived from an energy principle is important.It means that the operators are symmetric and that discretisation must lead to symmetric matrices.This is, of course, automatically implied if one begins with a finite element system, as was done in Bazant [34].The foregoing properties, which are not shared with the existing nonlocal continuum theory, represent the characteristic features of the imbricate nonlocal continuum.Similar, more generic observations were made in [35].

SPH discretisation
This section outlines the main aspects of the SPH method relevant to the considerations to follow.SPH uses kernel interpolation to approximate the field variables at any point in a domain.For instance, an estimate of the value of a function f (x) at the location x is given, in a continuous form, by an integral of the product of the function f (x) and a kernel (weighting where: the angle brackets ⟨⟩ denote a kernel approximation h is a parameter, known as the smoothing length, that defines size of the kernel support Ω x ′ is the additional independent variable.The kernel function W usually has the following properties: It is symmetrical, typically bell-shaped, function with compact support, i.e.W is zero everywhere except on a finite domain.In Classical SPH this domain is taken to be all points within twice the smoothing length, h, of the centre: The kernel gradient ( 13) is antisymmetric, i.e.
where is r = ⏐ ⏐ x − x ′ ⏐ ⏐ distance from the centre of the kernel support.These requirements, formulated by Lucy [1], ensure that the kernel function reduces to the Dirac delta function when h tends to zero, i.e. lim h→0 W ( ) and therefore, it follows that lim h→0 ⟨ f (x)⟩ = f (x).These requirements are consistent with the requirements for the nonlocal weighting function α (s) in Eringen's nonlocal continuum theory.
It is important to observe that the approximation (10), when applied to displacement and velocity fields in continuum, results in nonlocal continuum description.The classical, local continuum is obtained as a special case when the kernel support size h → 0, i.e. lim h→0 W ( ) and the kernel function tends to Dirac's delta function.
At this point, it is useful to recall Eringen's statement [36], In order for this type of approximation, to be applicable to an elastic continuum it must satisfy the following two requirements: Requirement I -If the stresses, σ(x), are everywhere zero, the strains in a stable material must be also zero, i.e., no unresisted deformation (zero -strain energy deformation mode) may be permitted by the theory.
Requirement II -In a stable material, the wave propagation velocity, v, must be real.
From these two requirements it follows that the Fourier transform of the weighting function, W ( ⏐ ⏐ x − x ′ ⏐ ⏐ , h), must be positive for all real k, i.e.
These requirements are referred to and used in the consideration of SPH stability in Appendices.
If the function f (x) is only known at N discrete points, the integral of Eq. ( 10) can be approximated by a summation: In the above equation, the subscript I and J denote particle number, m J and ρ J the mass and the density of particle J, N n the number of neighbours of particle I (number of particles that interact with particle I, i.e. within the kernel support), m J ρ J is the volume associated with the point/particle J and W (|x I − x J | , h).The discrete Eq. ( 15) constitutes the basis of the SPH method.The value of a variable at a particle, denoted by superscript I, is calculated by summing the contributions from a set of neighbouring particles (Fig. 2), denoted by superscript J and for which the kernel function is not zero: Note that in general ⟨ f (x , unlike in finite element approximations where shape functions satisfy the condition of being equal to 1 or 0 over the interpolation points, i.e.N I (x J ) = δ I J .Now we consider a kernel approximation of gradient ∇ f (x) of a scalar valued function, i.e.
Integrating by parts gives: Using Green's theorem, the first term on the right-hand side can be rewritten as: The surface integral is zero if the domain of integration is larger than or equal to the compact support of W or if the field variable assumes zero value on the boundary of the body (free surface).If none of these conditions are satisfied, modifications should be made to account for boundary conditions.When discretised over a set of points, contained in the kernel support, Eq. ( 19) becomes where the kernel gradient is evaluated for points J at distance r I J = ⏐ ⏐ x − x ′ ⏐ ⏐ from the centre I of the kernel domain However, in SPH codes, instead of the approximation (20) velocity gradient and stress divergence are typically approximated using alternative expressions: where v(x) in Eq. ( 21) is artificially added in order to ensure that the velocity gradient vanishes for a uniform velocity distribution [37].Similarly, σ(x) in Eq. ( 22) is added to ensure conservation of linear momentum.The standard derivation of these equations can be found in [38,39].Note that Eq. ( 21) represents gradient of relative velocity instead of velocity.Similarly, the difference form of Eq. ( 22) represents an approximation of divergence of stress resultant [σ(x At this point we refer again to Bazant at.al. [34,36] and the statement that strain measure should be physically interpreted through its relation to the change of distance between points x + s and x − s.In the case of the SPH approximation, given below for one dimension (1D) for clarity reasons, nonlocal strain ε is related to the change of distance between points x and x ′ (or x and x − s, s where D i represents differential operator as proposed by Bazant [34].The term [(u (x) − u (x − s)) W (s, h)] 2h −2h = 0 due to compact kernel support.
Note, Eq. ( 23) represents smoothed gradient of relative displacement and that in the SPH method the displacement field is not smoothed.Further, when the kernel support size h → 0, i.e. lim h→0 W ( ) and the kernel function tends to Dirac's delta function the nonlocal strain in Eq. ( 23) tends to zero ε (x) → 0, not to the classical local continuum strain as illustrated below.
This represents a very important difference from Eringen's nonlocal continuum where lim h→0 ε tends to local strain ε (x) as the kernel domain size tends to zero.The same observation holds for the velocity gradient approximation (21) and the difference form of the stress divergence approximation (22).
Finally, a discrete SPH form for ε can be stated as where I & J are particle labels not indexes.Bu differentiating Eq. ( 25) with respect to time the nonlocal strain rate can be defined in terms on particle velocities as In the derivation of the expressions for strain and strain rate, Eqs. ( 25) and ( 26) the kernel anti-symmetry, , was used.These expressions are well known SPH discrete forms originally derived using a heuristic approach to enforce Galilean invariance.An alternative explanation is offered in [40] where the SPH balance equations were derived in a moving reference frame.
When used in a rate form of a constitutive equation, the gradient of relative velocity between particles I and J given by Eq. (26) indicates that the nonlocal stress rate σ (x I ) at particle I depends on relative velocity between the particles.
Observe that the discrete nonlocal strain rate and stress rate (nonlocal effects) also tend to zero as the reduction of kernel domain size tends to zero, i.e. when the kernel function tends to Dirac delta function.This represents the main difference with Eringen's nonlocal continuum theory where this limit process results in the classical local continuum theory.
For the principle of virtual work in the SPH formalism we assume that the incremental work, δℵ, depends on the local strain, ε, and on the smoothed (nonlocal) strain, ε.Therefore, for a bar of length L virtual work is: where δℵ b is the work done by stresses for the boundary segments of length 2h at each bar end (where kernel support is incomplete); δℵ 0 , δℵ 1 are the works of nonlocal stresses σ and local stresses τ for the rest of the bar; δu(x), δε(x), δε(x) are any kinematically admissible variations; and c is the partition of unity coefficient with values 0 ≤ c ≤ 1, to be defined later.Note, in the expression for δℵ 1 x is the length coordinate of the bar; u is displacement.Substituting δℵ 0 and δℵ 1 into Eq.( 27) yields where F 1 and F 2 are certain functions which are zero for 2h ≤ x ≤ L − 2h and have effect only within end segments of length 2h.Since δℵ = 0 for any kinematically admissible δu(x), the expression in the braces { } of the first integral must be zero.This yields the one-dimensional equation of motion given in (30).
The bar of a unit cross section is subdivided into N equal segments of length ∆ p by particles K = 1, 2, 3, . . ., N + 1 as shown in Fig. 3.The bar is also divided into K overlapping segments of length 4h which defines the size of the kernel support, i.e. representative volume.The nonlocal strain ε, given by Eq. ( 25), obtained in one-dimensional overlapping domains of length 4h, is used to determine nonlocal stress σ (x) for each particle.These domains, shown as circles in Fig. 3, overlap at each particle.The number of overlapping segments (kernel domains) for a particle I is equal to number of particles within the domain of particle I, i.e. n = N n + 1 where ∆ p is number of neighbours of particle I (n includes particle I).In order to continue with the discussion of the SPH approximation it is necessary to recall the concept of area vectors within the SPH method originally outlined by Swegle [41] which is based on the fundamental definition of the stress tensor where a resultant force F acting on a surface A with the unit normal n due to stress σ is given by F = σ • nA.He showed that the SPH approximation for stress divergence (smoothed stress divergence) can be rewritten in terms of a particle interaction area as: where A I J = V J ∇W I J is the interaction area vector which has a direction normal to the surface and a magnitude equal to the area of the surface (shown in Fig. 4) and V is the particle volume.Note I and J indicate particles not vector or tensor components and that A I J is defined for each IJ pair.Thus, the gradient of the kernel function can be thought of as defining the area on which stress acts to produce a force between the particles.
When this concept is applied in 1D and more specifically to the bar illustrated in Fig. 3 where each particle has only two neighbours one on each side, A I J represents the IJ pair interaction area and at the same time the bar cross section.If the interparticle distance ∆ p is decreased so that the number of neighbours is increased to 4, the interaction areas A I J for the two particle pairs (particle I and the neighbouring particles) on each side of particle I can be determined.This can be interpreted as two overlapping bar elements acting in parallel on particle I.The combined cross-section area of all overlapping elements intersected by the same cross section of bar should be a constant and independent of the number of overlapping elements.

Variational derivation for difference operator (3D)
As before, we start by adopting Hypothesis I Bazant at.al. [33] and assume that the stress depends on the change of distance between two particles I and J, |x I − x J |, in this case in 3D.
Therefore, the stress depends on the kernel approximation of the relative displacement gradient, D ⊗ u (x) = D i u j (x), where D is a differentiation operator defined as: Note, D ⊗ u (x) represents smoothed gradient of relative displacements, i.e.
Consequently, nonlocal strain ε (x) defined using the differential operator D represents a nonlocal smoothed strain measure and can be expressed in a continuous and then discretised form as where i and j are indices, I and J are particle labels, x is coordinate vector; u is displacement vector; and 2h is the material characteristic length (a material property).
The work of the variations δ where σ i j is the nonlocal stress tensor (repeated indices imply summation).Based on symmetry of σ i j we can write, ε in Eq. ( 34) defines the nonlocal smoothed strain tensor.According to Hypothesis I [33], the variation of the total work in the body, δℵ, depends on the nonlocal stress σ i j and nonlocal strain ε i j as well as on the local stress τ i j and local strain ε i j .Therefore where: Ω ′ is the domain of the body without a boundary layer of thickness 2h, Fig. 5 Ω b is domain of boundary layer of thickness 2h, Ω = Ω ′ ∪ Ω b , see Fig. 5 S 0 is surface of the body as shown in Fig. 5, [33] dv = d x 1 d x 2 d x 3 is differential volume element δu i is any kinematically admissible virtual displacement f i is distributed volume force in Ω ′ cδℵ 0 is the virtual work of the local stresses τ i j within domain Ω ′ (1 − c) δℵ 1 is the virtual work of the nonlocal stress σ i j within domain Ω ′ ds is differential surface element The combined volumes of n imbricate elements overlapping at a point must be constant and independent of n.Denoting this combined volume as (1 − c) where 0 ≤ c ≤ 1, we see that the volume contribution of each imbricate element is (1−c) /n.Obviously, it decreases as the subdivision of V is refined and tends to 0 as n → ∞.
Note, c characterises the fraction of the material behaving in a local manner.For c = 0, we have the nonlocal continuum and for c = 1, we have the classical local continuum as special cases.There is a number of possible ways of determining c, two possibilities are discussed in Section 5.
Let us now try to separate δu i , in the integrand of δℵ 1 .Using Eq. ( 34) we can write where repeated subscripts (both i and j) imply summation over 1, 2, 3.
where is Ω is kernel domain (support).Note that σ is evaluated at x and x ′ and that, according to Gauss theorem, Eq. ( 38) results from the smoothed divergence of the stress resultant as shown in Fig. 6, i.e.
is not smoothed and for that reason the symbol σ is not used.However, if instead of the differential operator D we use conventional local differentiation the expression for δℵ 1 has the form given by Eq. (38) which is not consistent with the definition of the nonlocal strain.Note that this is the form most commonly used in the SPH balance of linear momentum [37][38][39].We proceed with virtual work δℵ 1 for Ω ′ and δℵ b for Ω b as The term δℵ 0 by Gauss' integral theorem, can be restated in which S is the boundary surface of Ω ; n is the unit normal of the surface; and ds is the surface element of S and τ (ε).
The description of continuum motion based on local stress τ and nonlocal stress σ allows for two different constitutive equations.The nonlocal constitutive equation for σ offers the possibility to include nonlocal aspects of material behaviour which cannot be incorporated into the local constitutive equation for τ, see for instance [42].
There are few possible ways to evaluate local stress τ and strain ε including the possibility to calculate τ by smoothing nonlocal stress σ or by using local strain calculated using normalised corrected kernel interpolation.The latter option is described below.
The local strain ε and the local stress divergence ∇ • τ are calculated using normalised corrected SPH as proposed in [43].The expression used are given below for completeness.Zero order consistent Shepard functions W (|x I − x J | , h) are used for approximation of field variables (as kernels).For instance, displacements are approximated as And the corrected approximations of differential operators (first order consistent), e.g.displacement gradient is approximated as where B is and finally, ⟨∇u(x)⟩ in discretised form is where Eq. ( 43) applied to particle velocities is used to approximate local velocity gradient.Finally, local stress divergence is approximated as These expressions complete definitions of all discrete operators required to discretise the equations of balance of mass, momentum and energy and at the same time to calculate parameters c and c ρ .
According to d'Alembert principle, we may also substitute f = −ρ ü, where ρ is mass density and ü = ∂ 2 u /∂t 2 , t is time.Thus, Eq. ( 45) where ψ i , Φ j , Ψ i are functions independent of δu j .Note that only the first of the integrals in Eq. ( 46) involve the interior domain Ω ′ , which excludes the boundary layer Ω b .According to the principle of virtual work, Eq. ( 46) must hold for any kinematically admissible variation, δu j (x).Choosing δu j (x) to be zero outside Ω ′ , and nonzero and arbitrary within Ω ′ , it follows from the fundamental lemma of the variational calculus that the expression in parentheses in the first integral must vanish for all x.This yields the continuum equation of motion: which is a partial difference-differential equation.Note that this equation holds only within Ω ′ , i.e., it does not apply within the boundary layer Ω b .
The continuum equations of motion for the surface and the boundary layer Ω b follow, mainly, from the second to fifth integrals in Eq. ( 46).However, their form requires special attention and it is addressed separately.

Determination of the partition of unity parameter c
As already stated, c characterises the fraction of the material behaving in a local manner.For c = 0, we have the nonlocal continuum and for c = 1, we have the classical local continuum as special cases.One possible way of determining c for an elastic continuum is based on the requirement that the local and nonlocal stresses combined should be able to define a constant stress field (i.e. are in equilibrium).We start with the assumption that the effects of the nonlocal stress σ i j ( ε i j ) are equal to the effects of the local stress τ i j ( ε i j ) , i.e. (1 − c) σ i j = c • τ i j where the local strain ε i j is calculated using normalised corrected SPH as proposed in [43].
In Eq. ( 48), representing nonlocal and local stress equilibrium, the fourth order tensors C i jkl and C i jkl define local and nonlocal materials stiffness properties.In general, C i jkl and C i jkl are different [32].However, in this consideration in order to simplify Eq. ( 48) we assume that local C i jkl and nonlocal C i jkl are the same.This allows for c to be defined as: To be consistent in discretisation of the balance laws, we have to use local/nonlocal form of the balance of mass equation where D • v is nonlocal velocity divergence and c ρ characterises the fraction of the material undergoing a local change of density.Assuming that the local and nonlocal density in a homogeneous body have to be the same c ρ can be determined as An alternative way to determine c ρ based on the assumption of volume conservation can be as follows:

SPH stability assessment for elastic continuum
The instability of the Classical SPH method which typically demonstrates itself in tension (consequently named tensile instability) has been one of the main shortcomings and sources of criticism of the SPH method.In addition to tensile instability SPH, similar to the finite element method, exhibits the instability related to zero strain energy modes of deformation.Here we limit our attention to two key publications related to SPH stability analysis without any disrespect to other researchers who worked on this problem.
The first comprehensive stability assessment of the SPH method for a 1D elastic continuum was performed by Swegle [25].This von Neuman analysis resulted in the well-known instability criterion where the positive product of stress and second derivative of the kernel function, i.e.W ′′ σ > 0, implies instability.
Six years later, Belytschko at al. [44] published a unified stability analysis of meshless particle methods which included assessment of SPH, again for a 1D elastic continuum, and confirmed Swegle's findings.It is important to notice that, unlike in Swegle's analysis, the SPH stability analysis in [44] was based on the constitutive equation integrated over the domain (neighbourhood) of the particle I instead of domains (neighbourhoods) of the neighbouring J particles.The effects of this assumption on the stability analysis were not considered in the paper.
These two publications considered a version of the SPH method based on the combination of nonlocal derivatives of velocity and local derivatives of stress.Their findings are directly relevant to the Classical SPH form of the balance of momentum equation typically used in contemporary SPH codes.
The SPH stability analysis presented here is performed in two stages.In this section we present only the main findings and provide the detailed analyses in Appendices A and B.
Firstly, we considered stability of the continuous nonlocal and local-nonlocal SPH approximations of the elastic continuum following [32,36] in order to establish stability requirements for the kernel functions used in the SPH approximations.In order for the SPH approximations, to be applicable to an elastic continuum they must satisfy the following two requirements [36]: Requirement I -If the stresses, σ(x), are everywhere zero, the strains in a stable material must be also zero, i.e., no unresisted deformation (zero -strain energy deformation mode) may be permitted by the theory.
Requirement II -In a stable material, the wave propagation velocity, v, must be real with an appropriate dispersion relation.
The consideration of the Requirement I showed that the Fourier transform of the SPH kernel function W * (k) must be positive for all k including k = 0.
From the dispersion relationship we see that, in the case of nonlocal continuum (c = 0), velocity v is real and the continuum is stable if and only if the autocorrelation of the SPH smoothing function |W * (k)| 2 is positive for all k except k = 0.The local continuum, (c = 1) is unconditionally stable.Consequently, the partition of unity parameter c can be used to stabilise potentially unstable behaviour of nonlocal continuum.
Details of the analysis which led to these conclusions are given in Appendix A. Second, we assessed stability of the discretised nonlocal SPH form by repeating Swegle's [25] analysis.Swegle's criterion for instability in terms of the stress state and the second derivative of the kernel function in the case of nonlocal SPH is W ′′ σ < 0 which is reverse of the original Swegle's instability criterion (illustrated in Fig. 7(a)).The analysis showed that the method's stability in tension and in compression coincides with the stability of elastic continuum (illustrated in Fig. 7(b)).The instability of the Classical SPH originating from an effective stress with a negative modulus being produced by the interaction between the constitutive relation and the kernel function was not present in the case of nonlocal SPH.
This analysis is presented step by step for a 1D elastic continuum in Appendix B.

Numerical examples
This section provides a number of 1D and 2D examples designed to investigate and illustrate behaviour of the various SPH formulations derived in the previous sections when implemented in a numerical code.This includes the discrete local, nonlocal SPH and mixed local-nonlocal forms when implemented in a numerical code.As already stated, the local form used to approximate velocity gradient is given in Eq. ( 43) and the form used to approximation stress divergence is given in Eq. (44).
The discrete form used to approximate non-local relative velocity gradient is: Note that Eq. ( 53) is identical to the Classical SPH approximation of velocity gradient.The non-local form of the momentum equation, consistent with the approximation of gradient of relative velocity (53) is given in Eq. ( 54).This form is directly related to Eq. ( 37) in the derivation of nonlocal SPH.
⟨ dv I dt The second form of the balance of momentum equation considered in this section is the Classical SPH Eq. ( 55) In addition to nonlocal SPH and Classical SPH we also considered two mixed formulations based on partition of unity.The first mixed formulation is local-nonlocal form and the second mixed formulation is local-Classical SPH form.
In both mixed formulations the balance of mass and the balance of momentum equations include a partition of unity parameter, so the density is updated as And the acceleration is updated as Through the choice of partition of unity parameter c for the momentum equation and c ρ for the balance of mass equation the behaviour of the mixed local forms as well as the special cases of the non-local continuum and local continuum were investigated.The numerical code uses the second-order explicit central difference algorithm for time integration.The algorithm is identical to the one used by the LS-DYNA code [45].This algorithm advances the position of all particles from time t n to time t n+1 using where the particle velocity is updated using the acceleration determined from Eq. (56), A linear hypo-elastic constitutive model is used for the stress update.The local stress, used in the local form of the momentum equation and the non-local stress used in the non-local form of the momentum equation were calculated and stored separately within the code.In the results presented below, where stress is plotted it is the total stress, i.e.
For the uniform tension/compression test cases and the Swegle stability test no artificial viscosity is used in order not to mask any effects of instability.For the subsequent wave propagation and impact test, the standard Monaghan and Gingold artificial viscosity term [5] is used.The smoothing length is kept constant throughout the simulations.

Uniform tension and compression in 1D
In these tests a 1D bar of length 1.0 is given a linearly varying initial velocity (v = ±0.001x) to generate a uniform strain over the length of the bar.In this example the origin of the x-axis was located at the mid-point of the   and 9.These results show good agreement with the analytical solutions in both tension and compression.The one exception was the sum form (Classical SPH) for which the calculation was unstable and did not reach a solution time of t = 10.0.
The values of the partition of unity parameter that are calculated for the combined case remain essentially constant at c ρ = 0.501 and c = 0.502.

Swegle numerical stability test in 2D
Swegle demonstrated the SPH tensile instability by applying a perturbation velocity to a single particle in the centre of a square 2D domain (solid block) prestressed with tensile or compressive stress [25].In this numerical example we used a direct copy of Swegle's numerical test.The initial particle distribution is a two dimensional array of 441 particles arranged in a rectangular lattice with particle coordinates between −1.0 ≤ x ≤ 1.0 and −1.0 ≤ y ≤ 1.0 with a particle spacing of 0.1, as shown in Fig. 10.The loading was generated by defining the initial material density as a fraction of the referential density.If this fraction is greater than 1.0 then the material is exposed to a constant isotropic compressive stress, if the fraction is smaller than 1.0 then the material is exposed to a constant isotropic tensile stress.An initial velocity perturbation of v x = 1.0 × 10 −10 is applied to the centre particle.The outer three layers of particles were fully fixed throughout the calculations to remove any influence of the domain boundary.
The material properties are representative of hydrodynamic (zero yield strength) aluminium with density ρ 0 = 2.8, Bulk Modulus K = 0.686275, Shear Modulus = 0.0.All artificial viscosities were switched off.To consider stability in this simulation the absolute velocity of the perturbed particle was recorded.If the behaviour is stable the velocity amplitude remains below the value of the initial perturbation.If the behaviour is unstable then the amplitude grows and the particles move in a non-physical way.
For the sum form (Classical SPH) the nonphysical particle motion (clumping) is illustrated in Fig. 11 for two response times t = 30 and t = 42.The velocity history for the perturbed (central) particle is given in Fig. 12(a) displaying the behaviour observed by Swegle [25], i.e. stable behaviour in compression and unstable behaviour in tension.Note that the rate of growth of the instability is proportional to the level of stress.
Stable behaviour of the non-local difference form in Swegle's test is illustrated in Fig. 12(b) showing no growth of the perturbed particle velocity magnitude.
The same test was also performed on the mixed local-Classical SPH form to investigate stabilising effect of the local SPH on the Classical SPH form.The test results are given for two levels of background stress, i.e.Fig. 13(a) for ρ 0 = 0.990 and (b) for ρ 0 = 0.950.The partition of unity parameter c calculated was determined according to Eq. ( 49) varied in the range 0.497< c < 0.502 and consequently the curves for c = 0.5 and c-calculated are almost identical.The increase of c had stabilising effect, manifested in a reduced rate of instability growth.

A 3D bar under uniform compression and tension
This numerical example was designed to validate and illustrate behaviour of the mixed local-nonlocal form in 3D.For that purpose, a 3D model of an elastic bar with a constant cross section under the uniform tension and compression test was generated, see Fig. 14.The boundary conditions and number of particles along the X-axis was 100 as in the 1D test above.For the free surfaces Y = ±0.035and Z = ±0.035symmetry plane boundary conditions were used to generate the 1D strain state in the bar.A 7 by 7 particle discretisation was used in the Y and Z directions to ensure that the particles along central region of the bar have no ghost particle neighbours used to enforce the symmetry boundary conditions.The resulting stress distributions for the central row of particles along the x-axis for the response times t = 1.0 and t = 10.0 are given in Fig. 15 for compression and in Fig. 16 for tension.These 3D results agree with the results from the equivalent test in 1D given above.

Compressive wave propagation tests with free surface boundary
This wave propagation test was performed in 1D.This test consisted of a 5 cm long bar impacting a rigid wall at the origin.This rigid wall was modelled with a symmetry plane.The bar initial velocity was 50 m/s.The material was linear elastic with Young's modulus equal to 200 GPa and Poisson's ratio equal to 0.3.The material density was 8000 kg/m 3 .For this wave propagation case, and all subsequent impact test results, Gingold and Monaghan artificial viscosity term [5] was used in this test case, with a coefficient of 0.5, and the bar was discretised with 500 particles.The ratio of smoothing length to particle spacing was 1.3.This test case was intended firstly to verify basic performance of the different SPH formulations described in this paper, and second to illustrate the need for a different treatment of particles in the boundary layer, as described in Section 2.
Upon impact with the rigid wall a compressive stress wave was generated and propagated into the bar.Fig. 17 shows the results before the compressive wave reaches the free end for the Classical SPH (Equations 1.39 and 1.40) and the nonlocal SPH (Equations 1.38 and 1.39), the purely local formulation (Equations 1.31 and 1.32), and mixed     of different values of the partition of unity coefficient c is also illustrated in Fig. 17.As expected the curves lie between the curves for c values equal to 0 and equal to 1.
As pointed out in Section 2, the nonlocal formulation requires special treatment of the boundary layer.This is illustrated in Fig. 18 which shows the same wave propagation problem at a later response time after the wave has reflected from the free end and the release wave has travel some distance into the bar.The analytical solution, the thick solid line, in Fig. 18(a) shows the section of the bar that has unloaded at response time t = 10 µs.It is clear from this figure that the response on the local, nonlocal and mixed local-nonlocal formulations is non-physical.The exception, as expected, is result for Classical SPH, which does reproduce the free end boundary condition.
Based on this property of Classical SPH the following approach to treatment of the free boundary was investigated.The Classical SPH approximation of stress divergence (plus form) for particles in the boundary layer and the nonlocal SPH (difference form) for the remainder of the domain were used.Specifically, in this example the Classical SPH form was applied to the last four particles at the free end.The results for this approach are shown in Fig. 18(b).It is evident that this simple boundary treatment worked effectively in this 1D example as the results for all formulations now closely match the analytical solution.

Tensile wave propagation
This second wave propagation test problem, schematically represented in Fig. 19, was adapted from the wave propagation problem that was used by Baǎzant and Belytschko in [1,33].In these papers the constitutive behaviour of the bar included a strain softening response, but for the purpose of this work a pure linear elastic model was more appropriate to study the propagation of tensile stress waves.The Young's modulus was 70.Therefore, two tensile step waves are generated at the two ends which travel towards the centre of the bar.When the two waves meet at the centre they superpose resulting in two times higher stress magnitude which propagates from the centre towards the bar ends.
The problem was modelled in 1D and 3D.In the 1D case the bar had a cross section area A = 1, and in 3D the models had a 1-by-1 cm cross section with the use of symmetry planes in the +Y/−Y and +Z/−Z directions  to generate a 1D strain state, as shown in Fig. 20.The bar was discretised with 200 particles along the length for both 1D and 3D models, with boundary layers of seven particles at each end where the velocity was prescribed (see Fig. 20).
In both 1D and 3D cases (see Fig. 21(a) and (b)) it can be clearly observed that the Classical SPH (sum form) results in an unstable solution.For this form the instability starts to develop as soon as the tensile wave is generated at the ends of the bar and then propagates into the bar as the wave progresses.As can be seen in Fig. 21 after 20 µs the response has become unphysical with large errors in the stress levels across the bar, and including large unphysical oscillations.In contrast, the non-local difference form and the mixed local-non-local formulation does not show instability, and match the analytical solution closely.The effect of using different values of the partition of unity parameter c is illustrated in Fig. 22.It can be seen that irrespective of the value used the response appears stable.
The SPH solutions are also in close agreement with the equivalent FEA solution.The SPH wave front is not as steep as the FEA results as this is inherent to the SPH interpolation.The FEA results were obtained using the LS-Dyna software.The FEA mesh consisted of a line of 200 under integrated hexahedral elements.The symmetry boundary conditions were applied to the side surfaces of the bar to ensure a 1D strain deformation state.The prescribed velocities were applied at both ends of the bar.

2D and 3D impact tests
The purpose of these tests is to assess the behaviour of the local-nonlocal form on representative solid mechanics problems.This requires treatment of the free-surface boundary condition that is not represented in the nonlocal form  of the momentum equation (54).This is unlike the Classical SPH form of the momentum equation (55), where the error which is the consequence of the neighbour particle deficit at the boundary, results in approximation of the free surface.The wide use of this form is a consequence of this ability to represent free surface.As rigorous treatment of the free surface boundary condition for the local-nonlocal form is future work, in these tests the free surface behaviour was approximated through the use of a boundary layer of particles that use the Classical SPH rather than the nonlocal form of the momentum equation.Consequently, this boundary layer is unstable in tension.
The first test example is the 2D rubber shell or 'tennis ball' impact originally used by Swegle and co-workers [25] to illustrate SPH behaviour in tension.This problem consists of two 4 cm outer diameter, 3 cm inner diameter, shells each with an initial velocity of 50 m/s.The material properties are density 1.01 g/cm 3 , elastic modulus 439.9 MPa and Poisson's ratio 0.4.In the results presented here we use 40 particles through the shell thickness, Fig. 23, this is significantly higher than the 10 through the thickness used in the original Swegle results, the reason is that the significant tensile strains that occur on the inner shell surface close to the initial contact point result in instability in the boundary layer particles that prevents the calculation form continuing.The use of 40 particles ensures that the boundary layer thickness is small compared to the total shell thickness, allowing the calculation to continue successfully.The boundary layer is two particles thick, so representing the particles that have a reduced number of neighbour particles due to the surface.The results for the rubber shell are shown in Fig. 24, demonstrating that the local-nonlocal form does allow the shells to contact and separate without failure.
The second example is the 3D symmetric impact of two cylinders, Fig. 25.The purpose of this test was to further illustrate the behaviour of the local-nonlocal form in 3D.Each cylinder has length 4 cm and diameter 1.9 cm with an inter-particle distance of approximately 0.1 cm, resulting in 11,360 particles per cylinder.The cylinders use a linear elastic material model with density 2.8 g/cm 3 , elastic modulus 1.0 GPa and Poisson's ratio 0.3.Each cylinder has an initial velocity of 50 m/s.Results for this case are shown in Fig. 26, showing the impact and separation of the cylinders.The level of deformation in this analysis is small as increasing the impact velocity or reducing the   elastic modulus resulted in instability in the boundary layer across the contact interface, the region of greatest tensile stress.The interaction of this instability with the contact algorithm resulted in error termination of the analyses.
In order to investigate the 3D behaviour under larger deformation, the model was modified to consist of a single cylinder and a symmetry plane, removing the need for boundary particles across the impacting surface.To increase deformation the initial velocity of the cylinder was increased to 100 m/s and the elastic modulus decreased to 100 MPa.Results for this case are shown in Fig. 27, showing successful treatment of the larger elastic deformations in this case.Instability did begin in the tensile region of boundary layer shortly after the final time shown.
The results from these three cases verify the capability of the 3D local-nonlocal form for solid mechanics applications.Robust treatment of the free-surface boundary condition remains a future challenge.

Conclusions
In this paper we have revisited the derivation of the Eulerian version of the SPH method.In addition to reviewing the Classical SPH and the normalised corrected/local forms of SPH, we provide a consistent derivation of a nonlocal These developments are applicable to Eulerian SPH characterised by variable neighbourhood approximations which allow for modelling of large deformations.
Our stability analysis of the continuous nonlocal SPH approximation for an elastic continuum defined a stability requirement that the auto correlation of the kernel function has to be positive.We then demonstrate that the discrete form of nonlocal SPH, unlike Classical SPH, does not suffer from tensile instability.
However, the nonlocal SPH requires additional treatment of free surface boundary conditions.In the numerical examples for the nonlocal SPH, free surface boundary conditions were imposed by applying the Classical SPH approximation for stress divergence over a boundary layer of particles (particles with incomplete kernel support), i.e. these examples were treated as volume constrained problems.
The mixed SPH forms offer a way of modelling local and nonlocal effects in continuum, e.g.progressive localised damage through use of local and nonlocal constitutive models.LNL-SPH does not suffer from tensile instability.In CLL-SPH the local SPH improves stability of Classical SPH.
Further work is required in derivation of rigorous approach to treatment of free surface boundary conditions in nonlocal SPH and LNL-SPH, as well as development of an efficient numerical implementation of the LNL form to minimise the additional computational cost of this form.In addition, the proposed nonlocal and mixed formulations can be extended to SPH in the Total Lagrangian framework.
It is important to observe that the approximation (A.60), of the relative displacement and relative velocity fields in the continuum, results in nonlocal continuum description.
where ⟨⟩ indicates kernel approximation/smoothing.And the related nonlocal stress can be labelled the nonlocal elastic moduli C i jkm where the SPH kernel W ( ⏐ ⏐ x − x ′ ⏐ ⏐ , h) corresponds to the weighting function in nonlocal theory of elasticity [28].In the SPH approximation of the equation of motion the smoothing integral (A.60) is applied to the right-hand side of the equation of motion, i.e. the term Note, when the kernel support size h → 0, the kernel function tends to Dirac's delta function lim h→0 W ( ) ⇒ 0 and the nonlocal strain ε i j (x) ⇒ 0 leading to the classical, local continuum as a special case, i.e.
In order for this type of approximation, to be applicable to an elastic continuum it must satisfy the following two requirements [36]: Requirement I -If the stresses, σ(x), are everywhere zero, the strains in a stable material must be also zero, i.e., no unresisted deformation (zero -strain energy deformation mode) may be permitted by the theory.
Requirement II -In a stable material, the wave propagation velocity, v, must be real with an appropriate dispersion relation.It is shown below that from these two requirements it follows that the Fourier transform of the kernel function and its differential, W ( , must be positive for all real k.To prove this statement and for the sake of clarity, we consider an infinite 1D elastic body subjected to deformation, for which Eq. (A.61) may be simplified as According to Requirement I σ (x) = 0 ⇔ ε(x) = 0 and using Eq.(A.62), The nonlocal strain should be zero, i.e.Eq. (A.63) should be satisfied, only when relative displacements u R ( are equal to zero (no deformation).A general displacement may be approximated as u (x) = a exp (ikx) which yields where a is a real amplitude, and k > 0 is a real constant, should satisfy Eq. (A.63).
Substitution of Eq. (A.64) into Eq.(A.63) yields the condition that the equation must not have any solution, which after dividing by a exp (ikx) becomes where W * (k, h) and ) respectively as given in Eq. (A.66).
An alternative approach, presented below, based on Fourier transform of W (s, h) leads to the same requirement.
So, we can conclude that the kernel function Fourier transform W * (k, h) must not be zero for any k.Hence, W * (k, h) must be continuous and it must be either positive or negative everywhere on Ω .The requirement that W * (k, h) must be positive cannot be proven without Requirement II.
Consider the Requirement II.We restrict attention to small deformations, such that in 1D local strain is ε (x) = ∂u(x) ∂ x , where u (x) is displacement.The equation of motion for a nonlocal continuum is , where t is time and ρ is mass density.Substitution of Eq. (A.62) into the equation of motion then yields where nonlocal terms are where |W * (k)| 2 is autocorrelation of the SPH smoothing function From the dispersion relationship Eq. (A.72) we see that, in the case of nonlocal continuum (c = 0), v is real and the continuum is stable if and only if the autocorrelation of the SPH smoothing function |W * (k)| 2 is positive for all k except k = 0.In addition, from the Requirement I the Fourier transform W * (k) must be positive also for k = 0, as already proven.The local continuum, (c = 1) is unconditionally stable.Consequently, the partition of unity parameter c can be used to stabilise potentially unstable behaviour of nonlocal continuum.
It is important to note the following two points.
• Requirement I is not applicable to inelastic materials with residual stresses.However, Requirement I could then be replaced, with equal results, by the requirement that, if the stress rate, σ , is zero, the strain rate ε must also be zero, provided that the tangent modulus, E t = ∂σ ∂ε , is positive.• The above stability assessment and the derived stability requirements apply to the continuous form of the SPH method.Discretised form of SPH is characterised with additional stability requirements for details see Swegle at al. [25] and Belytschko at al. [44].

Appendix B. Stability assessment of the discretised SPH method for elastic continuum
The main assumption used in the analysis are as follows.The smoothing length is equal to the initial interparticle distance, consequently only the nearest neighbours contribute to the particle sums.The next-nearest neighbours located at distance 2h do not, since both the kernel function and its derivative are zero for |x J − x I | ≥ 2h.Density is assumed to be constant and unaffected by small perturbations.The dependence of the kernel function on the smoothing length is not considered.Further, instead of using symmetrised form Eq. (B.73).
the balance of linear momentum is approximated by Eq. (B.74) in order to simplifies the stability analysis.
where Π I J is artificial viscosity required to smooth shock waves which has to be consistent with numerical resolution [46,47], µ I J is relative velocity of particles I and J, α and β linear and quadratic term coefficients.Based on the above assumptions Eq. (B.74) can be restated as , the particles have been numbered in order of increasing position in the x direction, so that I − 1 is the index of the nearest neighbour in the negative direction, while I + 1 is the index of the nearest neighbour in the positive direction.In Eq. (B.75) anti-symmetry of the differential of the kernel function was used.In 1D, W has dimension of length −1 , while the mass m should be interpreted as mass per unit length, with the cross-sectional area numerically equal to one.The dimension of W ′ is consequently length −2 .
It is important to notice that Eq. (B.75) would be the same if instead of Eq. (B.74) the simplification is performed on Eq. (B.73) characterised by the stress plus form (Classical SPH), i.e.
This leads to the conclusion that the results of Swegle's stability analysis apply to the sum form of the momentum equation.
If the same simplification is applied to Eq. ( 31) we obtain the form of the momentum equation compatible with and convenient for the stability analysis In order to simplify and shorten the writing of Eq. (B.76) we introduce the following approximations/substitutions Which, after ime derivatives are approximated by central difference expressions and the second substitution T I +1 = σ I +1 + Q I,I +1 , the equation assumes the form ẋI,n+ Following Swegle's approach, stress is defined as a function of density as where K is bulk modulus, c is speed of sound and η is the volume strain.Eq. (B.77) is used to update positions of particles according to the following expression Perturbation Propagation Equations -The next step is to obtain the linearised equations of first variation, which describe the propagation of small perturbations in the original equations.This is done by applying perturbations of the form x → x + δx to each of the variables, subtracting the unperturbed equations and keeping only the terms which are first order in δx.All unperturbed quantities in the equations are frozen; the coefficients of the perturbations are considered constants.The resulting equations are known as the perturbation propagation equations.
Each position and velocity, regardless of its spatial and temporal index, is replaced by a perturbed value given in Eq. (B.80).The equations describing the evolution of the perturbations are obtained by substituting the perturbed quantities into the discretised equation of balance of linear momentum (B.77), retaining only those terms linear in the perturbations, and subtracting the original equations.
And the perturbation of the kernel derivative is approximated using a first-order Taylor series expansion of W ′ and is independent of the form of the kernel function.W ′′ is the second derivative of W with respect to its argument, and is thus the slope of W ′ .describes the spatial variation of the perturbed quantities.In this expression X is the initial, or Lagrangian coordinate, so that ∆X is the initial uniform spacing between particles, and E I represents E being raised to the power I, rather than an index.The position of particle I at time zero is given by The perturbation is assumed to be periodic with wavenumber, k, which is the circular inverse of the perturbation wavelength, k, so that  The first term in brackets multiplying δx n involves kernel variations at constant stress, while the second term in brackets multiplying δx n involves stress variations at a constant value of the kernel.The term multiplying it δ ẋn− 1 2 involves the artificial viscosity.The shortest wavelength perturbation which can be resolved by the discrete system is λ min = 2∆X .At this wavelength the term involving stress variations goes to zero, so that the variation of stress has no effect on the propagation of perturbations at the shortest wavelength.For the resulting system of equations, we find the amplification matrix A, defined by U n+1 = AU n where U n+1 is the vector of values at the new time step, and U is the vector of values at the old time step.The eigenvalues of A, which depend on the wavenumber of the perturbation, determine the stability of the system of equations.If the largest eigenvalue exceeds unity, the amplitude of the value vector is exponential in time.
The perturbation propagation equations can be rewritten in the form δ ẋn+ 1 Therefore, a sufficient condition for unstable growth of the shortest wavelength (twice the particle spacing) is where T is negative in compression and positive in tension.On the other hand, if the system is conditionally stable, which means that the time step must be limited in order to achieve stability.

Fig. 2 .
Fig. 2. Central particle I (in red) and its neighbours (in blue) used in the SPH approximation.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 3 .Fig. 4 .
Fig. 3. SPH discretisation of an elastic bar of length L, 4h -kernel support size, ∆p -inter-particle distance with illustration of overlapping kernel domains, Particle I shown in red with its neighbours shown in black.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 5 .
Fig. 5. Various domains of integration, Ω kernel support, Ω ′ domain discretised with complete kernel supports Ω , Ω b boundary layer domain which cannot be discretised with complete kernel supports Ω , S 0 surface of the domain Ω b ∪ Ω .

Fig. 8 .
Fig. 8. Stress profiles for 1D uniform compression problem at the response time t = 1.0(a) and the response time t = 10.0 (b).

Fig. 9 .
Fig. 9. Stress profiles for 1D uniform tension problem at the response time t = 1.0(a) and the response time t = 10.0 (b).Results for the Classical SPH form are not shown due to instability.

Fig. 10 .
Fig. 10.Initial rectangular lattice used in the stability tests.

Fig. 11 .
Fig. 11.Plot showing particle positions and x (horizontal) velocity at the response timet = 30 (a) and the response timet = 42 (b) for the sum form with initial density ρ0 = 0.990.

Fig. 12 .
Fig. 12. Velocity amplitude for the perturbed (central) particle (a) for Classical SPH (sum form) showing presence of tensile instability, (b) for the non-local (difference) form showing no tensile instability.

Fig. 13 .
Fig. 13.Velocity amplitude for the perturbed (central) particle for the unstable (mixed) Classical-local SPH form (a) for ρ 0 = 0.990 and (b) for ρ 0 = 0.950 illustrating stabilising influence of the partition of unity parameter c.

Fig. 14 .
Fig. 14.Initial particle distribution for 3D bar model of uniform compression and tension models.

Fig. 15 .
Fig. 15.Stress distribution for the uniform bar compression for 1D and 3D at the response time t = 1.0(a) and the response time t = 10.0 (b).

Fig. 16 .
Fig. 16.Stress distribution for the uniform bar tension for 1D and 3D at the response time t = 1.0(a) and the response time t = 10.0 (b).

Fig. 17 .
Fig. 17.Compression wave travelling from the rigid wall at the origin into the bar (towards the right) for local-no-local formulation with different values of the c coefficient: 0, 0.25, 0.5, 0.75, 1.0 and calculated using density-based criterion at response time t = 10 µs.

Fig. 18 .
Fig. 18.Release wave travelling from the bar free end towards the rigid wall (from right to left) for Classical SPH, local SPH, nonlocal SPH and local-nonlocal SPH: (a) without boundary treatment, (b) with boundary treatment, at response time t = 10 µs.

8
GPa, Poisson's ratio 0.108 and the density 1550 kg/m 3 .A bar 0.2 m long bar was loaded at both ends with a constant velocity v = 70 m/s (x = −0.1 m, v = −70 m/s; x = 0.1 m, v = 70 m/s) as shown in Fig. 20.

Fig. 22 .
Fig. 22.Effect of different values of c parameter on results at response time t = 20 µs.

Fig. 23 .
Fig. 23.Detail view of initial particle distribution for elastic shell problem.Boundary particles are highlighted in red.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 26 .
Fig. 26.Results for 3D symmetric elastic cylinder impact, showing impact and separation.Particles coloured by von Mises stress (Mbar).Half the cylinder is blanked to show the internal stress distribution.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 27 .
Fig. 27.Results for 3D elastic cylinder impact.Particles coloured by von Mises stress (Mbar).Half the cylinder is blanked to show the internal stress distribution.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

2 = ( 1 − 2 + s∆tδx n −∆tδ ẋn+ 1 2 +( 2 −
r ∆t)δ ẋn− 1 δx n+1 = δx n Amplification Matrix Eigenvalues -Eq.(B.99) can be rewritten in matrix form asLU n+1 = RU n (B.101)where the vector of new velocities and positions is this set of equations is determined by the eigenvalues of the amplification matrix A, in the equation U n+1 = AU n .From Eq. (B.101) it follows thatA = L −1 R (B.105)Determination of the eigenvalues of A can be simplified by noting that|A − λI| = 0 ⇔ |L| |A − λI| = |R − λL| = 0 (B.106)Thus, an equivalent procedure is to find the eigenvalues ofR − λL = [ (1 − r ∆t − λ) s∆t ∆tλ 1 − λ ] (B.107)The resulting eigenvalue equation isλ 2 + ( r ∆t − s∆t 2 − 2 ) λ + 1 + r ∆t = 0 (B.108)The system is unstable if the larger root of Eq. (B.108) exceeds unity, i.e. max {λ 1 , λ 2 } > 1.This leads to exponential growth of the perturbations.Determination of stability criterion can be simplified by writing the eigenvalue equation in the following formλ 2Bλ + C = 0 (B.109)whereB = 1 + θ, θ = r ∆t + s∆t 2 2 C = 1 + r ∆t (B.110)The value of the maximum eigenvalue depends on the value of the discriminant D, whereD = B 2 − C (B.111)so that the eigenvalues are given byλ 1,2 = B ± √ D (B.112)There are several cases to consider based on the sign of D and the magnitude of B, but in the current analysis all cases reduce to the statement that D > 0 is a sufficient condition for instability.Using Eq. (B.110) D can be expressed asD = s∆t 2 + θ 2 (B.113)so a sufficient condition for instability iss > 0 (B.114)Having in mind that the stability limit is controlled by λ min the parameter s becomes