The Taylor-SPH Meshfree Method: Basis and Validation
H. Idder1, M. Mabssout1, *, M. I. Herreros2
1Faculty of Science and Technology, Laboratory of Mechanics and Civil Engineering, Tangier, Morocco
2CEDEX, Madrid, Spain
To cite this article:
H. Idder, M. Mabssout, M. I. Herreros. The Taylor-SPH Meshfree Method: Basis and Validation. Applied and Computational Mathematics. Vol. 4, No. 4, 2015, pp. 286-295. doi: 10.11648/j.acm.20150404.17
Abstract: This paper presents the basis and validation of the Taylor-SPH meshless method formulated in terms of stresses and velocities which can be applied to Solid Dynamic problems. The proposed method consists of applying first the time discretization by means of a Taylor series expansion in two steps and a corrected SPH method for the space discretization. In order to avoid numerical instabilities, two different sets of particles are used in the time discretization. To validate the Taylor-SPH method, it has been applied to solve the propagation of shock waves in elastic materials and the results have been compared with those obtained with a corrected SPH discretization combined with a 4th order Runge-Kutta time integration. The Taylor-SPH method is shown to be stable, robust and efficient and it provides more accurate results than those obtained with the standard SPH along with the Runge-Kutta time integration scheme. Numerical dispersion and diffusion are eliminated and only a reduced number of particles is required to obtain accurate results.
Keywords: Taylor-SPH (TSPH), Meshfree, Runge-Kutta, Shock Wave, Stability, Dynamics
Mesh based numerical methods, such as the Finite Element Method (FEM), have been widely used to solve Solid Mechanics problems. Despite of their great success, mesh based numerical methods suffer from some difficulties when dealing with problems where large deformations, discontinuities and crack propagation are involved.
To overcome these difficulties, some meshfree methods were developed in the past. They are based on constructing approximating functions at arbitrary discrete points in the domain, and, as a consequence, the need for a mesh is eliminated. One example of meshfree method is the Smoothed Particle Hydrodynamics (SPH). The SPH method was developed in 1977 by Lucy  and Gingold and Monaghan  to be applied in Astrophysics. In 1991, it was extended to solve problems in Solid Mechanics where traditionally FEM fails .
It is well known that the original SPH method suffers from tensile instabilities  and lack of consistency . Over the past years, different corrections have been introduced to improve the accuracy of the SPH solution. Liu et al  proposed a correction function that restores the first order completeness of the kernel function. Some years later, Randles and Libersky  developed a transformation of the gradient that provides the correct values of the strains for linear fields; this normalization of the gradient is a generalization of the work proposed before by Johnson and Bessel .
To solve the tensile instability, Dyka et al   introduced the stress point method into SPH. This approach was later extended to higher dimensions by Randles et al. . Belytschko et al.  showed that the stress point technique stabilizes SPH by removing the instability that arises due to rank deﬁciency, i.e. spurious singular modes, while tensile instabilities can be avoided by using a Lagrangian kernel. Other recent studies to improve the convergence of the SPH method can be found in  .
Modelling of shock waves propagation in solids is most demanding area since a wide variety of problems are encountered. Even in the linear range, the short wavelengths are subjected to numerical damping and dispersion. The shocks are smoothed and leading or trailing oscillations appear. The situation is more involved in the non-linear regime, where strain can localize along shear bands which may be affected by the mesh if FEM is used in the analysis.
The key advantage of using SPH is that the absence of a mesh enables it to deal with larger local distortion than Finite Element Methods, and therefore phenomena such as fracture, and other material instabilities are more easily modelled. In addition to that, mesh dependence problems such as mesh alignment or mesh size dependence, are solved in a straightforward manner, given the meshfree nature of the method.
Nevertheless, the SPH method still presents several difficulties when dealing with dynamics and shock wave propagation: numerical damping, depending on relative wavelengths, and numerical dispersion, consisting of short waves propagating faster or slower and causing oscillations close to the front.
The authors of the present paper have published in previous works some different alternatives to solve the propagation of shock waves in solids using FEM . In Solid Dynamics, the classical approach of numerical analysis is based on a displacement formulation using the Newmark scheme. It is well established that the Newmark family of methods is not suitable to solve discontinuous phenomena, and many efforts have been devoted in the past to solve this problem. The shortcomings of the classical approach can be classified in: (a) Numerical damping and numerical dispersion; (b) In classical displacement formulations low-order elements cannot be used as they are not accurate and the results depend on mesh alignment and mesh size. To solve the problems mentioned in (b) the authors have proposed in previous works a numerical model formulated in terms of stress and velocity based on a Taylor–Galerkin scheme. This method has been proved to perform well for localized failure of viscoplastic materials of Von Mises type. It was shown that low-order elements, such as linear triangles, exhibited an excellent behaviour and reasonably numerical diffusion and dispersion.
In this paper, a new meshfree method (Taylor-SPH) for solving the propagation of shock waves in linear or non-linear media is presented. It uses a two-steps time discretization algorithm by means of a Taylor series expansion and a corrected SPH method for the spatial discretization. In order to avoid numerical instabilities, two different sets of particles are considered to perform the time discretization and a Lagrangian kernel is used. Both, Lagrangian kernel and its gradient, are corrected to satisfy the consistency conditions.
In order to validate the proposed method, the same problem of a shock wave propagating in an elastic material is solved using the 4th order Runge Kutta time integration scheme combined with a corrected SPH for the space discretization.
The purpose of this paper is to show how the proposed Taylor-SPH method provides solutions of more accuracy than those obtained with classical SPH methods.
The paper is organized as follows. First, the governing equations for dynamic problems in elastic media are given in terms of stress and velocity. In Section 3.1, equations are discretized using the proposed Taylor-SPH method. In Section 3.2, the partial differential equations are discretized using the 4th order Runge Kutta time integration combined with corrected SPH method. To assess the performance of the proposed method, some numerical applications in 1D and 2D are described in Section 4.
2. Governing Equations
Consider a linear elastic domain for with boundary. The governing equations are written in terms of stress and velocity. Prescribed stress and velocity are imposed on and respectively. Neglecting the body forces, and for small strains, the governing equations can be written as
where σ is the stress tensor, ρ the density and v the velocity vector.
are given functions.
In 2D problems, the differential operator S is defined as
The equation (1) has to be completed by a constitutive equation which gives the relation between the stress and the strain tensors. Here, the material behaviour is assumed to be elastic. In 2D problems, the constitutive equation is written as
In the case of 2D plane stress problems, the elastic matrix is given by
where E is the elastic modulus and n is the Poisson's coefficients.
In small strain analysis, the rate of strain tensor is related to the gradient of velocity as
being S the differential operator defined in (2).
Finally, the governing equations of the problem can be written as
For 2D plane stress problems, equations (6) can be written as
where Dij are the components of the plane stress elastic matrix.
The system of equations (7) can be written in a conservation form as
being U the unknown vector and F the flux vector.
It can be seen that equation (8) represents a system of first order hyperbolic equations.
It is interesting to note that in 1D problem, equations (7) can be written as
The velocity of wave propagation is given by, which correspond to two waves propagating in opposite directions.
3. Numerical Discretization
To solve numerically the system of partial differential equations (8), the traditional SPH method applies first the SPH spatial discretization, obtaining a set of simultaneous ordinary differential equations with respect to time, and this set of equations is then integrated in time using one of the standard techniques such as Euler, predictor-corrector or one of the Runge-Kutta schemes. Nevertheless, in the case of discontinuous functions, such as shock waves, these standard methods still present some numerical problems, such as numerical dispersion and diffusion close to the discontinuities.
To overcome the problem of shock waves propagation, the authors propose in this paper the use of an alternative method: the Taylor-SPH method (TSPH)   , which consists of applying first the time discretization by means of a Taylor series expansion in two steps and thereafter the spatial discretization using a corrected SPH.
3.1. Proposed Numerical Method: Taylor-SPH
3.1.1. Taylor-SPH Time Discretization
Time discretization of equation (8) is carried out by means of a Taylor series expansion in time of U up to second order accuracy:
The first order time derivative of the unknowns can be calculated using equation (8)
The second order derivative with respect to time is given by
First step: In order to obtain the time derivatives of fluxes at time tⁿ, the values of the unknowns at an intermediate time tn+1/2 will be obtained first
Using the computed value of Un+1/2, fluxes can be evaluated as
which can now be substituted in equation (13), resulting on
Second step: Substituting now the expressions obtained for the first and second order time derivatives, (12) and (16), in the Taylor series expansion (11), the following expression is obtained for the values of the unknowns at time tn+1
3.1.2. Spatial Discretization by the Corrected SPH Method
Once the equations are discretized in time, the corrected SPH method will be used to discretize the equations in space. A brief summary of the basic SPH method and its corrected form is presented.
(i). Integral Approximations in SPH Method
SPH is based on the concept of the integral approximation of a given function f and its derivatives. In SPH the integral approximation of a function f(x) is defined by
where the brackets denote the integral approximation, x and x’ are vectors, W(x-x’,h) is the kernel function and h is a measure of the size of the kernel support. It is clear that the kernel approximation converges to the exact function f(x) if the kernel function W(x-x’,h) tends to the Dirac function d(x-x’). The accuracy of the SPH method depends on the properties of the kernel. It must fulfil the following conditions:
- The kernel function must be positive
- The kernel function must be normalized over its support domain
- The kernel function must have a compact support
being k a constant positive parameter
- The kernel function must verify the Dirac delta function as h approaches to zero
- The kernel function must be a symmetric function of (x-x').
The SPH approximation can be formulated in terms of the Eulerian kernel or the Lagrangian kernel. It has been shown by Belytschko et al.  that the Lagrangian kernel eliminates the tensile instability, and therefore it will be used here. Using the Lagrangian kernel, the size of the kernel support remains constant h=ho during the simulation.
In this work the B-spline function based on the cubic spline function, previously used by Monaghan and Lattanzio , has been chosen as the kernel function:
being; the scaling factor C is given by , and in 1D, 2D and 3D respectively.
Concerning the SPH integral representation of the spatial derivative of a function f(x), it is given by
where is the derivative of W with respect to x´.
(ii). Corrected form of the SPH Discrete Approximation of Function
In the SPH method, a continuum is represented by a set of particles, thus the SPH approximation of a function f at point I is given by
where the summation subscript J denotes a particle label and runs over all particles N inside the domain, such that .
In equation (22), denotes the value of the kernel centred at node I at position J; fI = f(xI) and mJ and rJ are the mass and density associated to particle J.
It is well known that for boundary particles, the consistency conditions are not satisfied when using approximation (22). For a particle J near the boundary, it support extends out of the problem domain W and the kernel function is truncated by the boundary. Therefore the normality condition of the kernel is no longer valid. To overcome this problem, it has been used here the corrective kernel approximation of a function
where is the volume associated to particle J and is given by
It is important to note here that the denominator of above equation is unity for those particles whose support domain does not intersect the boundary. This correction of the approximating function satisfies the zeroth order completeness.
(iii). Corrected form of the SPH Discrete Approximation of Derivatives
Taking into account that , the discrete form of (21) can be written as
where ; ; ; and
Since, the completeness is given by the order of the polynomial which can be represented exactly; equation (25) satisfies only the zeroth order derivative completeness condition.
In order to fulfil the first order completeness, the corrected form of the derivative of the approximating function must be used:
This modification of the SPH approximation for the gradient of a function f, allows the method to fulfil the first-order completeness condition, enabling the derivatives of constant or linear fields to be reproduced exactly.
3.1.3. Discretized Equations Using the Proposed Method: Taylor-SPH
The time discretization is carried out in two steps by means of a Taylor series expansion. To perform the time discretization, it will be necessary the use of an auxiliary set of particles that will be called "virtual" particles. These "virtual" particles will be interspersed among the "real" particles in a similar manner it was done in stress-point integration methods. Fig. 1 shows the arrangement of "real" and "virtual" particles in one and two dimensions.
Thus, time discretization of model equations is carried out in two steps:
- In the first step, the values of the field variables at time tn+1/2 are computed by means of equation (14) at the positions of the Nv "virtual" particles.
- In the second step, the values of the field variables at time tn+1 are computed, using equation (17), at the positions of the Nr "real" particles.
First step: Applying the corrected SPH spatial discretization to the first step of time discretization, (14), we obtain
The subscript VP indicates that the values of U and are computed at the positions of the "virtual" particles.
Using the corrected form for the approximation of derivatives given by equation (26), we obtain the values of the variable U at tn+1/2
being J the "real" particles, such that .
Second step: Applying the corrected SPH spatial discretization to equation (17), we obtain
The subscript RP indicates that the values of U and are computed at the positions of the "real" particles.
Using expression (26), the values of the variable U at tn+1 are given by
where J are the "virtual" particles, such that .
In the 2D procedure, a structured quadrilateral particle arrangement is used. As it is shown in Fig. 2; 4 "real" particles are arranged so they form a square and a "virtual" particle for the calculations at tn+1/2 is placed in the centroid of the square.
In this particular case, computation of particle volumes is especially easy. The square is then subdivided into 4 other squares for computation of the volumes, and once the coordinates of the "real" and "virtual" particles are known, quadrilateral volumes can be easily computed as it is shown in Fig. 2. It is important to note here that in the first step, (29), only volumes of "real" particles are considered, while in the second step, (31), only volumes of "virtual" particles are taken into account. Note that particle masses and volumes need to be computed only once for a Lagrangian kernel.
It is also important to mention that the method presented here does not require any special treatment of the boundary conditions.
3.2. Runge-Kutta Time Discretization for SPH
To solve the system of partial differential equations (8), the traditional SPH method applies first the SPH spatial discretization, obtaining a set of simultaneous Ordinary Differential Equations, and this set of equations is then integrated in time using one of the standard techniques, such as the Runge-Kutta schemes.
3.2.1. SPH Discretization
Applying the corrected SPH discretization presented in Section 3.1.2 to equation (8), it results in the following semi-discrete equation
where being J the neighbouring particles around particle I, such that .
3.2.2. Time Integration: Runge-Kutta Schemes
The basic idea of Runge–Kutta (RK) methods is to evaluate the right-hand side of the differential system of equations given above (32) at several values of the unknowns UI in the interval between tn and tn+1, and to combine them in order to obtain a high order approximation of the field variables at tn+1. The number of intermediate values is referred to as the "Runge–Kutta stages".
Applying the general explicit M-stage Runge-Kutta method to solve (32) it results in
Where Dt = tn+1- tnand = UI (tn), and the coefficients aij, bj and cj () define the accuracy and stability of a given Runge-Kutta scheme. These coefficients are often written in tabular format known as the Butcher table 
Second order Runge-Kutta scheme: To integrate in time the equation (32) by using the second order Runge-Kutta scheme, it is necessary to perform two stages (M = 2) with b1 = b2 = 1/2; c2 = a21 = 1:
Fourth order Runge-Kutta scheme: To obtain the value of UI at time tn+1 from equation (32) by using the 4th order Runge-Kutta scheme, it is necessary to perform four stages (M = 4) with the coefficients given in the Butcher table 
4. Numerical Examples
The purpose of this section is to present some examples which will show the performance of the Taylor SPH method. The examples which will be considered next have been chosen in order to illustrate the following points:
(i) The TSPH method eliminates numerical instabilities
(ii) The TSPH method presents good wave propagation properties
(iii) The TSPH method provides solutions of more accuracy than those obtained using the Runge-Kutta time integration schemes for the corrected SPH.
(iv) The TSPH method presents good performance in bending dominated situations using a reduced number of particles
4.1. Propagation of a Shock Wave on a 1D Elastic Bar
In this section it will be solved the problem of a shock wave propagating in 1D elastic bar, using the proposed method (TSPH). The results will be compared with those obtained using the Runge–Kutta time integration schemes with the corrected SPH.
The problem consists of a bar of length L = 1 m with a unit section. The initial and boundary conditions are the following:
·Dirichlet boundary conditions:
The material properties are the density and the elastic modulus . With these parameters, the wave speed is
The bar has been spatially discretized by 51 "real" particles and 50 "virtual" particles for TSPH. The distance between two consecutive "real" particles is .
There is an analytical solution available for this elastic problem: the incoming wave will propagate towards the right boundary without any distortion and keeping its initial amplitude of and . It will reflect at the fixed end (L = 1 m) and the amplitude of s at this point will be doubled to a value of , while the velocity of the wave after reflexion will propagate along the bar with .
Figs. 3-4 show the stress and the velocity after reflection of the wave using the TSPH method. The time-step used for the analysis has been chosen to be which corresponds to a Courant number. It can be observed that no instability appears, neither numerical dispersion nor numerical diffusion appear in the front of the wave and the results are in complete agreement with the analytical solution.
In order to show the better performance of the TSPH in comparison with the Runge-Kutta time integration schemes, the same example of the propagation of shock wave in an elastic bar is solved using fourth order Runge-Kutta (RK4) and TSPH schemes under the following initial and boundary conditions.
·Dirichlet boundary conditions: and
With the TSPH method, the bar is discretized with 101 "real" particles and 100 "virtual" particles. The distance between two consecutive "real" particles is . The time step is . However, using the 4th order Runge-Kutta time integration with a corrected SPH, the bar is discretized with 101 particles and the time step used is
The results of the comparison are presented in Figs. 5-6. Fig. 5 shows the comparison for stress at point x=1 m. The comparison for velocity at point x=0.5m is given in Fig. 6. It can be observed that when using the 4th order Runge-Kutta time integration schemes with a corrected SPH, numerical dispersion and diffusion are present. On the contrary, the numerical solution using the TSPH method is free of oscillations and diffusion and the result is in complete agreement with the analytical solution.
These results show the good performance and efficiency of the Taylor-SPH (TSPH) method in comparison with the corrected SPH with a Runge-Kutta time integration (RK-SPH).
4.2. Numerical Stability of Taylor-SPH Method
The aim of this section is to study the numerical stability of the Taylor-SPH method with respect to the smoothing length ho and the time step Dt.
4.2.1. Stability Analysis with respect to the Smoothing Length
It is well known that the smoothing length is very important in the SPH method since it has direct influence on the efficiency of the computation and the accuracy of the solution. In the present work a Lagrangian kernel is used and the smoothing length ho has a constant value.
In order to accomplish a sensitivity analysis of the proposed method (TSPH) with respect to the smoothing length ho, the problem of propagation shock wave in an elastic 1D bar of length L = 1 m has been solved considering a fixed distribution of 51 "real" particles and 50 "virtual" particles. The time step used for the analysis is Dt=10-4 s. The value of parameter ho has been gradually increased and stability analysis has been carried out using the L2 norm to estimate the error.
The error estimation using the L2 norm is given by
where vexact is the exact solution and vhis the numerical solution given by Taylor-SPH.
Fig. 7 gives the error as a function of ho/Dx using the TSPH method. The values for h0/Dx considered in this study are within the range of 0.6-2.5. For lower values than ho/Dx=0.6 there are not enough particles within the support domain to accomplish the approximation. It can be observed that when ho/Dx is within the range of 0.6-1.5 the error is about 10-3%, the solution preserves its accuracy and it is in good agreement with the analytical solution. If the ratio ho/Dx is increased the solution loses its accuracy until it becomes highly distorted.
4.2.2. Stability Analysis with respect to the Time Step
In order to perform the stability analysis of the Taylor-SPH method with respect to the time step or the Courant number, a similar example as above has been treated. The parameter ho/Dx has been kept constant equal to 1.5 and the time step has been gradually increased from 10-5s to 10-4s. The stability analysis of the numerical solution as a function of the Courant number is given in Fig. 8. It can be observed that when the time step is Dt=10-4s which corresponds to Courant number the accuracy of the numerical solution gets its maximum with an error of about 10-3%. As the value of the time step is decreased, the solution loses its accuracy becoming oscillatory. For C > 1 the numerical solution becomes unstable.
4.3. Propagation of a Shock Wave on a 2D Elastic Bar
To investigate the performance of the TSPH method, a simple 2D example is studied. It consists of two opposite velocity shock waves which propagate into a 2D elastic bar. The bar is 1m long and 0.1m width. The material properties are density and elastic modulus . A rectangular impulse with a velocity of 1m/s is applied to both ends of the bar. The bar is spatially discretized using a structured particle arrangement of 306 "real" particles. The distance between two consecutive "real" particles is 0.02m. Fig. 9 shows the horizontal stress at time t=2ms for both methods: Taylor-SPH and the 4th order Runge-Kutta scheme (RK-SPH). It can be observed that the results are very dispersive when using RK-SPH while TSPH provides a more accurate solution and minimum dispersion.
4.4. Bending Problem
The purpose of this example is to show the behaviour of the proposed method in bending dominated situations. Consider a cantilever beam subjected to a vertical load is Fo = 8 102 N at the free end. Material is assumed to be elastic, with elastic modulus E = 8 109 Pa , Poisson's ratio ν = 0.0, density r=2 10³ Kg/m³ and the beam dimensions are L = 1 m and b = 0.02 m. The vertical displacement at x=L calculated by the TSPH method is compared to the exact solution ymax = 0.05 m. A structured arrangement of 63 "real" particles and 40 "virtual" particles is used as it is shown in Fig. 10.
Fig. 11 depicts the vertical displacement of the beam obtained with the TSPH method. The computed numerical value of the maximum deflection obtained at the right end is 4.9 10-2 m, which is in very good agreement with the analytical solution.
As it can be observed, the result obtained with the proposed algorithm reaches an acceptable accuracy with a small number of degrees of freedom (63 "real" particles) as well as it seems to be free of instabilities.
In order to study the improvement in the numerical solution for an increasing number of particles a convergence analysis has been accomplished. Table 1 shows the vertical displacement at the right end of the beam as a function of the increasing number of particles.
In this analysis it is shown that the TSPH method converges quickly to the analytical solution and only a small number of particles is required to obtain very accurate results.
|"Real" particles||"Virtual" particles||ymax exact||ymax TSPH||Error (%)|
A new meshfree method Taylor-SPH for dynamic problems has been presented. The method consists of a two-steps time discretization scheme based on a Taylor series expansion of the stress and velocity fields using a corrected SPH. Two different sets of particles have been used for the computations at each time step and a Lagrangian kernel has been used in order to avoid numerical instabilities.
The problem of the propagation of a shock wave in an elastic bar has been analyzed using the proposed method. The results have been compared with those obtained using the 4th order RK-SPH, and it has been demonstrated how the proposed TSPH method provides more accurate solutions, since numerical problems such as diffusion and dispersion are absent.
It has been proved that the proposed method:
- eliminates numerical instabilities
- presents good shock wave propagation properties
- provides solutions of more accuracy than those obtained using the Runge-Kutta time integration schemes for the corrected SPH.
- presents good performance in bending dominated situations using a reduced number of particles
In addition to that, the method has been shown to be stable, robust and only a reduced number of particles is required to obtain accurate results.