Updated Lagrangian Taylor-SPH method for elastic dynamic problems

This paper presents a discussion on the properties of the collocation meshfree method, the Updated Lagrangian Taylor-SPH (UL-TSPH), for dynamic problems in solid mechanics. The PDEs are written in mixed form in terms of stress and velocity for the elastodynamics problems. Two sets of particles are used to discretize the partial differential equations, resulting on avoiding the tensile instability inherent to classical SPH formulations. Numerical examples ranging from propagation of a shock wave in an elastic bar to a stationary Mode-I semi-Inﬁnite cracked plate subjected to uniaxial tension are used to assess the performance of the proposed method.


Introduction
The Finite Element Method (FEM) is the most popular mesh-based method used in computation solid mechanics, and many references can be found in the literature. Although the FEM has proved successful in modelling a range of problems, FEM has still a few drawbacks for modelling large deformations problems, presence the discontinuities and dynamic fracture.
To overcome the difficulties encountered by FEM in large deformation and problems involving discontinuities, several meshfree methods have been developed in the past being used in computational solid mechanics with a great success. In contrast to FEM, meshfree methods do not require a mesh connectivity procedure and they can perform numerical solutions without using a computational mesh-grid since they are based only on particles where field variables and their derivatives are estimated. Therefore, meshfree methods can deal in a simple manner with failure problems and large deformation without the problems encountered in mesh-based methods where mesh refinement is usually required.
Generally, meshfree methods are formulated using either the Galerkin meshfree methods based on the weak formulation of the PDEs or the collocation meshfree methods based on the strong formulation of the PDEs.
Galerkin meshfree methods require a special technique to implement the essential boundary conditions and the selection of a suitable numerical integration scheme, while the collocation meshfree methods are based on the strong form of the PDEs which leads to efficient formulations without numerical integration and essential boundary condition procedures.
A wide variety of meshfree methods have been proposed over the past years including the Element Free Galerkin Method [3], Reproducing Kernel particle Method (RKPM) [14], Smoothed Particle Hydrodynamics (SPH) developed by Lucy [17] and Gingold and Monaghan [8] and the Taylor-SPH (TSPH) method developed by Herreros and Mabssout [9,11,18,19]. The SPH method has been extended to problems involving fluid-structure interaction [1,12]. In this paper, we will use the TSPH method, which has been successfully applied to solve linear and nonlinear problems [9,11,18,19].
Generally, the SPH method suffers from two principal problems: (i) the lack of consistency on the boundary that can lead to poor accuracy. Various corrections have been developed to improve the precision of the SPH method and to restore the zeroth and first order completeness of the kernel function [2,4,10,13,16,21] and (ii) tensile instability which appears in dynamic problems with material strength as reported by Swegle et al. [23]. To eliminate SPH tensile instabilities, Dyka et al. [5,6] proposed a stress point approach where the unknowns such as density, stress and energy are estimated at some prescribed stress points placed among the SPH particles. This approach was later extended to two dimensions by Randles et al. [22].
Taylor-SPH (TSPH) meshfree method developed by the authors over the past years [9,11,18,19] has been demonstrated to exhibit a good performance when applied to dynamic problems. The TSPH is a collocation meshfree method where the system of equations is written as a set of first-order hyperbolic partial differential equations. The method firstly consists a time discretization by mean of a Taylor series expansion in two steps and thereafter the space discretization using a corrected form of the SPH method. To solve the tensile instability, two different set of particles are used during simulation. It is worth mentioning that the classical SPH methods are written in simple terms of displacements and therefore Dirichlet boundary conditions for displacements and natural boundary conditions for traction must be considered. However, when using the TSPH meshfree method, the PDEs are formulated in mixed terms of stress and velocity, and thus only Dirichlet boundary conditions must be taken into account and are enforced in a simple manner at each particle on the border.
In previous works, the authors have extended the TSPH to large deformation problems [11]. The governing equation are written using the updated Lagrangian formulation. It was shown that the updated Lagrangian Taylor-SPH (UL-TSPH) produces very accurate results for solid mechanics involving large deformations and dynamic fracture. In this paper, the UL-TSPH is presented and it is applied to the propagation of the shock wave, to the vibration of a bar and to the elastodynamics crack problems, showing that the proposed numerical method is able to compute the Dynamic Intensity Factors and therefore it opens a new way to simulate dynamic fracture of a cracked material.
The paper is organized as follows. The partial differential equations using the Updated Lagrangian description are presented in Section 2. In Section 3, the PDEs are discretized using the UL-TSPH method. To show the performance of the proposed method, some numerical examples are presented in Section 4.

Mathematical model
The governing equations used in this work consists of (i) the momentum balance equations and (ii) a constitutive law defining the solid behavior. Here, we have restricted our analysis to an elastic material. In order to get an objective rate of the stress tensor, the Jaumann rate of the Cauchy stress is used. In terms of stress and velocity and neglecting the body forces, the governing equations are given by where ρ is the material density, σ is the Cauchy stress tensor, v is the velocity vector, d is the rate of deformation tensor, ω is the spin tensor and D e is the elastic tensor. The term S in (2) come from the Jaumann rate of the Cauchy stress, where ω is given by the antisymmetric part of the velocity gradient or and d is given by the symmetric part of the velocity gradient or Because the updated Lagrangian formulation is chosen here, the above equations are expressed in the current configuration which will be taken as the referential configuration.
In this paper, 2D plane stress problems are considered and therefore the stress and strain rate can be arranged in a matrix form as

Numerical model
The UL-TSPH method [11] is used to discretize (1) and (2). UL-TSPH consist of applying first the time discretization using a Taylor series expansion in two steps and thereafter the space discretization using a corrected form of the SPH method.

Time discretization
Time discretization of (1) and (2) are carried out using Taylor series expansion up to the second order of the unknowns v and stress σ To obtain the first order time derivatives of S, σ and d in (8)-(9), we compute first the unknowns σ and v at an intermediate time t Therefore, the first-order time derivatives of σ, S and d at time t n are computed as follows: Inserting (11) into (8) and (9), the unknowns (σ, v) at time t n+1 are computed as follows: The equations (10) and (12) represent the first and the second step, respectively, using UL-Taylor-SPH meshfree method.

Spatial discretization: Corrected SPH approximation
To perform the time discretization, the spatial discretization is required to obtain the numerical solution of the problem. A corrected form of SPH method is used to discretize (10) and (12) in space.
Here, a summary of the basic SPH method and its corrected form is presented in this section. SPH is based on the integral approximation of a given function f (x) and its derivatives. The kernel approximation of function f (x) by SPH is given by where W (x − x, h) is the kernel function and h is the smoothing length defining the size of the kernel support. The integral approximation of the gradient of a function f (x) is given by where The functions used as kernels in SPH approximations are required to fulfill some specific conditions [15]. They must be even function and must fulfill the following specific conditions: In (15b), k is a positive parameter defining the smoothing support size of the kernel function.
In this work, the cubic spline function (16), proposed by Monaghan and Lattanzio [20], has been chosen as the kernel function the scaling factor C is given by 1 h , 15 7πh 2 and 3 2πh 3 in 1D, 2D and 3D, respectively.
In the SPH method, a continuum is represented by a set of particles, thus the kernel integral (13) is approximated by a summation over a set of discretized particles where N is number of the neighboring particles of I such that represents the value of the kernel function at position J and Ω J = m J ρ J is the volume related to particle J; m J and ρ J are the mass and density of particle J, respectively.
The particle approximation for the spatial derivative (14) is given by where SPH approximations (17) and (18) do not fulfill the consistency conditions for a given particle near to the boundary of the domain. The kernel function is truncated at the boundary. Therefore, the normality condition of the kernel function and its gradient is no longer valid. To restore the zeroth order and first order completeness conditions, a corrected form of the SPH meshfree method based on the normalization of the kernel function and its derivatives is used in this work. Corrected form of the kernel approximation: In order to fulfill the C 0 consistency, the corrective kernel approximation for function f (x I ) is given by [2,4] The kernel approximation (19) represent constant functions and then it fulfills the C 0 consistency. Corrected form of the first derivatives approximations: In order to fulfill the C 1 consistency, the corrected SPH approximation of the gradient of function f is given by [2,4] where∇ and SPH approximation (20) satisfies the zeroth and the first order derivative consistency; that is, the approximation of the zeroth and first order polynomials is exactly approximated.

Updated Lagrangian Taylor-SPH: Spatial discretization
As noted above, the time discretization is carried out in two steps using a Taylor series expansion.
To perform this space discretization, the computational domainΩ = Ω∪∂Ω is discretized in two set of particles: {x k } NR 1 real particles and {x k } NV 1 virtual particles where N R and N V denote real and virtual particles. The virtual particles will be interspersed among the real particles in a similar manner it was done in stress-point integration [6,22] to stabilize the SPH method. Fig. 1 shows the distribution of real and virtual particles in 1D and 2D. Using the UL-TSPH algorithm, the spatial discretization is carried out in two steps: First step: Applying the corrected SPH spatial discretization (19) and (20) to the first-step time discretization (10) Second step: Applying the corrected SPH spatial discretization (19) and (20) to the second step time discretization (12) Finally, at the end of a certain time step, the positions of the real particles are updated where particle I refers to real particle.
The new position of the virtual particles is updated as the centroids of the Delaunay triangulation of the new real particles position.

Propagation of a shock wave on a 1D elastic bar: Numerical stability
In this section, we analyze the influence of the smoothing length and the Courant number on the convergence of the UL-TSPH.

Sensitivity of the numerical results with respect to the smoothing length
It is well known that in the SPH method, the smoothing length h plays an important role since it has direct influence on the accuracy of the solution [18].
If h is too small, there may be not enough particles number in the smoothing support of dimension kh to interact with a given particle, which results in very low precision. By contrast, if h is too large, local properties perhaps smoothed out, and the precision decreases too. In order to perform a sensitivity testing of the UL-TSPH method with respect to the smoothing length h, the problem of the propagation of a shock wave in an elastic bar has been analyzed.
The initial and Dirichlet boundary conditions are as follows: Factor α defines the radius of the kernel function which relates the interparticle spacing to the smoothing length. For the cubic spline kernel, the smoothing domain has radius 2h. The value of factor α has been gradually increased and the sensitivity of the numerical solution for different h has been studied. Factor α was given values of α = 1.2, 1.8 and 2.2, corresponding to 5, 7 and 9 neighbor real particles inside the support domain. Fig. 2 illustrates the stress σ(L, t) for the three values of α. It can be observed that for α = 1.2, the numerical solution converges well to the exact solution with an error about 2.5 × 10 −3 %. When α is increased (α = 1.8; 2.2) the numerical oscillations appear; the error increases to 25% and the solution diverges. Using the L 2 norm, the error of the numerical method as function of factor α is given in Table 1.  This example was calculated with a non-uniform particle distribution along the bar. Fig. 3 presents the stress history at point x = 0.44 m. As expected, the solution is stable but it exhibits small oscillations at discontinuities of the wave. These oscillations are present where the distance between two consecutive particles is not constant and the local Courant number is not equal to 1. For complex geometries where a non-uniform particle distribution is required, a local maximum time step satisfying the CFL condition can be used. That allows to avoid the use of the smallest time step imposed by the smallest distance between two consecutive particles for the whole domain.

Sensitivity of the numerical results with respect to the Courant number
In order to accomplish a sensitivity analysis of the UL-TSPH with respect to the Courant number, a similar study as above has been carried using different values of C. The factor α is kept constant equal to 1.2 (α = 1.2). The value of time step Δt has been gradually increased and the accuracy and stability of the numerical solution for different Courant number C have been tested. Fig. 4 shows the stress σ(L, t) for different Courant number after reflection of the wave using the proposed method. It can be observed that when the C= 1, the accuracy of the solution Fig. 4. Stress at the fixed end of the bar for different Courant number using UL-TSPH gets its maximum, the error is about 2.5 × 10 −3 %. As can we see, the numerical results loses its accuracy when the Courant number is decreased and the solution becoming oscillatory. For C > 1, the numerical solution becomes unstable. The error using the L 2 norm is given in Table 2 for different Courant number.   Fig. 6 gives time histories of kinetic energy, elastic strain energy and total energy obtained using the UL-TSPH method with C = 1. The proposed algorithm is shown to possess excellent energy conservation property.

Vibration of a bar with fixed ends
To show the performance of the proposed method, we consider a vibrating elastic bar fixed at both ends as shown in Fig. 7 where c = E/ρ = √ 2 m/s represents the elastic wave speed in the bar. The Fig. 8 presents the stress and displacement along the x axis using the analytical solution and the UL-TSPH method at time t = 0.353 8 s. The results are free of diffusion and dispersion and show an excellent agreement between the numerical and the exact solutions. Fig. 9 shows the contour of stress along the bar at different time. The results are stable and free of oscillations.
In order to check the convergence of the method, this example has been calculated for three particles spacing: 0.02 m; 0.01 m; 5 × 10 −3 m. The time step is kept constant equal to 3.53 × 10 −3 s. Table 3 presents the L 2 error for each case, showing clearly the convergence of the proposed method with a small number of particles.

Stationary Mode-I semi-infinite crack
In this example, the dynamic stress intensity factors (DSIFs) for an infinite plate with a semiinfinite lateral crack is evaluated using UL-TSPH meshfree method. The plate loaded by uniform tensile stress perpendicular to the crack-face with the amplitude σ 0 = 500 MPa as shown in Fig. 10. The applied boundary conditions are as follows: (ii) on Γ 2 ∪ Γ 3 , σ xy = v 1 = 0, (iii) finally, stress at Γ 4 is σ yy = σ(t).
The plate dimensions are taken as follows: length L = 10 m, vertical half-height of the plate H = 2 m and crack length a = 5 m. The material properties are taken as: Young's modulus E = 210 GPa, Poisson ratio ν = 0.3 and density ρ = 8 000 kg. The analytical solution is given by [7] K dy (0, t) = 0 if t < tc, where t c is the time needed by the stress wave to reach the crack-tip and given by where c d is the dilatational wave speed. The normalized DSIFs is given bỹ where r represents the distance from the crack tip along the line θ = 0 • . Fig. 11 presents the normalized DSIFs as a function of a normalized time obtained by using the UL-TSPH method. As can be seen, the numerical results converge well to the exact solution with 101 × 41 SPH particles.

Conclusion
A collocation meshfree method Updated Lagrangian-Taylor-SPH for large deformation in dynamic problems has been presented. The governing equations are written in mixed terms of stress and velocity. The time discretization is carried out using Taylor series expansion and the space discretization based on corrected form of SPH method.
From the examples presented in Section 4, the following can be deduced: • The UL-TSPH presents good shock wave propagation properties minimizing numerical diffusion and dispersion.
• It avoids numerical instabilities in a natural manner.
• It is an efficient tool to assess the Dynamic Stress Intensity Factors. Therefore, it opens a new way to simulate dynamic fracture of a cracked material.
• The method is efficient, robust, easy to implement and only few particles are needed to get accurate results.