Taylor-SPH Method for Viscoplastic Damage Material
Hajar Idder, Mokhtar Mabssout*
Laboratory of Mechanics and Civil Engineering, Faculty of Science and Technology, Abdelmalek Essaâdi University, Tangier, Morocco
To cite this article:
Hajar Idder, Mokhtar Mabssout. Taylor-SPH Method for Viscoplastic Damage Material. Applied and Computational Mathematics. Vol. 4, No. 3, 2015, pp. 162-173. doi: 10.11648/j.acm.20150403.19
Abstract: In this paper, we apply the meshless method Taylor-SPH to solve the propagation of shock wave in viscoplastic material coupled to damage. The equations are written in terms of stress and velocity. Taylor-SPH method is based on the Taylor series expansion of stress and velocity and on the corrected SPH approximation. Numerical stability of the method as a function of the smoothing length and the Courant number is analysed in the elastic case. The Taylor-SPH method is used to simulate localization in a one dimensional viscoplastic damage problem. The numerical results show that the Taylor-SPH method is able to model localization phenomena in viscoplastic damage material without lose of hyperbolicity of partial differential equations.
Keywords: Taylor-SPH, Meshless, Viscoplastic, Damage, Shock Wave, Stability
Finite Element Method (FEM) has been widely used to solve the Partial Differential Equations (PDEs) in Mechanic and Engineering problems. Despite of it great success, FEM suffers from some difficulties when dealing with problems where extremely large deformations, moving boundaries, discontinuities or crack propagation are involved. To overcome these difficulties, various meshless methods were developed in the 90s. The first meshless method is the Smoothed Particle Hydrodynamics (SPH). It was originally introduced by Lucy  and Gingold and Monaghan  to model astrophysical problems. The continuum is represented by a discrete set of particles. Each particle carries field variables such as velocity, stress, mass, density, etc. However original SPH suffers from some numerical difficulties: tensile instabilities  and lack of consistency . Over the last decades, many corrections have been introduced to improve the consistency and the accuracy of the SPH solution [5,6,7]. To eliminate the tensile instability, Dyka et al.  have proposed to insert additional points called stress points into SPH formulation. Nevertheless, the SPH method still presents some difficulties: numerical damping and dispersion when dealing with dynamics and shock wave propagation leading to poor accuracy of the solution. In the previous work, we have used FEM to solve the problem of propagation of shock wave in an elasto-viscoplastic material; the results present small oscillations [9,10].
In this paper, we present a new meshless method Taylor-SPH [11,12,13] for solving the propagation of shock waves in solids. This meshfree method 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 the tensile instability, two different sets of particles are used. Both Lagrangian kernel and its gradient are corrected to satisfy the consistency conditions.
The purpose of this article is to show that Taylor-SPH avoids the tensile instability, minimizes the numerical damping and dispersion and it performance to predict the damage of materials in dynamic conditions. To avoid the loss of hyperbolicity of PDEs during the damage evolution, viscoplastic law coupled with damage model is adopted.
The paper is organized as follows. In the next section, the stress-velocity mixed formulation for dynamic problems in elasto-viscoplastic damage material is presented. In Section 3, governing equations are discretized using the proposed Taylor-SPH method. To assess the performance of the proposed algorithm, some numerical examples in 1D are described in Section 4 using elasto-viscoplastic and damage materials.
2. Governing Equations
The governing equations used in this work consist of the balance of momentum equations and a constitutive law describing material behaviour. For the sake of clarity of the method, this study is limited to one dimensional case. Neglecting the body forces and for small strains, the governing equations can be written as follows.
2.1. Balance of Momentum Equations
For one dimensional bar of length L, the balance of momentum equation is given by
where is the stress, v is the velocity and ρ is the density.
2.2. Constitutive Law
In this work, an elasto-viscoplastic law coupling with damage has been chosen to describe the material behaviour.
2.2.1. Viscoplastic Law
For an elasto-viscoplastic law, the relation between the stress and strain rates is given by
where E is the Young’s modulus, e is the total strain andis the viscoplastic strain which is given by the Perzyna law :
are the Macaulay brackets such that: and g is a viscoplastic parameter.
j( f ) is a flow function chosen here as
N is a model parameter. In the case of Von Mises yield criterion, the functiondepends only on the second invariant of the deviatoric stress vector
The size of the yield stress fo will vary according to a suitable hardening or softening law. Here it will be assumed a linear dependence on the equivalent deviatoric plastic strain.
where H is the hardening modulus.
Taking into account that; the balance of momentum and constitutive equations can be written in 1D form as
The velocity of wave propagation is given by the eigenvalues of the matrix of the system (7): , which remain real.
2.2.2. Elasto-viscoplastic Damage Model
where D represents the isotropic damage variable . For the initial undamaged material D=0 and for completely failed material D=1. The damage evolution law used here is chosen as
is a model parameter, is the initial damage threshold and is the viscoplastic strain.
Differentiating eq. (8) with respect to time yields to
As damage increases in the material, the viscoplastic part of the model must incorporate the damage parameter. By writing the Von Mises yield criterion (5) based on the effective stress, the flow function (4) becomes
The system of first order PDEs for 1D elasto-viscoplastic damage model is given by
In compact form
It is important to notice that the eigenvalues of the matrix A are always real:
Therefore the system of eqs. (14) remains hyperbolic and the problem still well-posed during the damage evolution. The use of the viscoplastic law regularizes the damage model and avoids the PDEs to become ill-posed.
3. Numerical Discretization: Taylor-SPH Meshless Method
To solve the PDEs (14), the conventional SPH method applies first the SPH space discretization, obtaining a set of 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 Runge–Kutta schemes. In the case of discontinuous functions, such as shock waves, these standard methods present numerical problems, as numerical dispersion and diffusion, close to the discontinuities. To solve the problem of the shock waves propagation, we have used an alternative meshless method, Taylor-SPH [11,12,13] which consists of applying first the time discretization by means of a Taylor series expansion in two steps and thereafter the space discretization using a corrected SPH.
3.1. Time Discretization in Two Steps
Time discretization of Eq. (14) is carried out by means of a Taylor series expansion in time of U up to second order of accuracy:
The first order time derivative of the unknown can be calculated using Eq. (14)
The second order derivative with respect to time is given by
First step: In order to obtain the time derivatives at the RHS of expression (18) at time tn, the values of the unknowns at an intermediate time tn+1/2 will be obtained first
Using the computed value of Un+1/2, we have
Substituting (20) and (21) in equation (18), we obtain
Second step: Substituting now the expressions obtained for the first and second order time derivatives, (17) and (22), in the Taylor series expansion (16), we obtain the values of the unknowns at time tn+1:
3.2. Space Discretization by the Corrected SPH Method
In this section, a brief summary of the basic SPH method and its corrected form is presented.
3.2.1. Discrete Approximation of Function
In the SPH approximations, each function f(x) is represented by its integral approximation defined by:
where the brackets denote a kernel approximation; W(x-x’, h) is the kernel function; h is the smoothing length that defines the size of the kernel support and (x;x') are the coordinates vectors.
The Lagrangian kernel  is used here where the neighbours of influence do not change during the calculation, remaining h=ho as a constant value.
In this work the B-spline function introduced by Monaghan and Lattanzio  has been used as the kernel function:
Being ; the scaling factor C is given by , and in 1D, 2D and 3D respectively.
In the SPH method, a continuum is represented by a set of particles, thus it is necessary to approximate the integral (24) in a discrete manner:
where the summation subscript J denotes a particle label and runs over all particles N inside the domain, such that . WIJ = denotes the value of the kernel, centred at node I at position J, and are the mass and density associated to particle J and .
It is clear that for boundary particles, the consistency conditions are not satisfied when using approximation (26). Many works have been proposed in the past to improve the particle consistency of the SPH method [5,6,7]. The corrective kernel approximation of a function used here is :
The approximating function (27) satisfies the zeroth order completeness.
3.2.2. SPH Discrete Approximation of Derivatives
The SPH integral representation of the derivative of a function f(x) is given by:
Taking into account that, the discrete form of eq. (29) is:
where ; ; ; ;
The approximation (30) satisfies only the zeroth order derivative completeness condition. The corrected form of the kernel approximation for the derivatives is :
The approximation (31) fulfills the first-order completeness condition.
3.3. Taylor-SPH Discretization
To perform the time discretization presented in Section 3.1, it is necessary the use of an auxiliary set of particles called virtual particles. These virtual particles will be interspersed among the real particles as it is shown in Figure 1.
Therefore, the time discretization of model equations is carried out in two steps:
- In the first step, Un+1/2 is computed using (19) at the positions of the Nv virtual particles.
- In the second step, Un+1 is computed using (23) at the positions of the Nr real particles.
First Step: Applying the corrected SPH spatial discretization to equation (19), and using the corrected form for the SPH approximation given by (27) and (31), we obtain the vector UI at tn+1/2:
where I = virtual particle; J = real particles such that .
Second Step: Applying the corrected SPH spatial discretization to equation (23), and using expressions (27) and (31), we obtain the vector UI at time tn+1:
where I = real particles ; J = virtual particles, such that .
It is important to mention that with Taylor-SPH method, the PDEs are written in terms of stress and velocity and thus only essential boundary conditions must be considered.
4. Numerical Examples
The examples which will be considered in the following have been chosen in order to verify the following points:
• The Taylor-SPH method eliminates the SPH tensile instability
• The stability condition of the Taylor-SPH method
• The performance of the Taylor-SPH method concerning numerical damping and dispersion properties when dealing with shock wave propagation
• The Taylor-SPH method is useful to model the damaged material in dynamic
4.1. 1D Elastic Bar: SPH Tensile Instability
As for any numerical method, the first step is to ensure that the method is able to provide accurate results for problems having an analytical solution. Here, the Taylor-SPH method is applied to a shockwave propagation problem into a one-dimensional elastic bar. The problem has been sketched in Fig. 2 and consists of a bar of length L=0.1333 m with a unit section.
The applied boundary and initial conditions are the following:
- Boundary conditions:
- Initial conditions:
Therefore, the bar is initially under tension.
The material properties of the bar are: Young’s modulus and density This example has been treated by Dyka et al.  to address the tensile instability. Material properties, boundary conditions and discretization parameters are similar to those used by Dyka et al. . Solving this problem with a traditional SPH algorithm is impossible due to the tensile instability .
There is an analytical solution for this problem: the incoming wave will propagate towards the right boundary without any distortion and keeping its initial amplitude of . The wave speed is . It will reflect at the fixed end (x=L) and the amplitude of the stress at this point will be doubled, while the velocity of the wave after reflection will propagate along the bar with .
To apply the Taylor-SPH method, the bar has been discretized using 41 real particles and 40 virtual particles. The space between two consecutive real particles is. Dyka et al.  used the time step, therefore the Courant number is. Dyka el al. has chosen this time step as to not exceed the Courant stability limit. Fig. 3 depicts the stress at point obtained using the Taylor-SPH method. As it can be observed, with this Courant number the result presents numerical diffusion and oscillations which remain bounded during the simulation. The result shows that the Taylor-SPH method is stable.
The calculation has been repeated with Courant number C=1 and therefore the time step is Figs. 4 and 5 show the stress and the velocity at point respectively. The numerical solution using Taylor-SPH is free of oscillations and diffusion and the results are in complete agreement with the analytical solution. These results for the shock wave propagation problem in an elastic bar show the good performance of the Taylor-SPH method and demonstrate that the proposed method eliminates the SPH tensile instability and have a good shockwave propagation property.
4.2. Numerical Stability of Taylor-SPH
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 courant number C. The case study is another problem of propagation of a velocity shock wave in a 1D elastic bar. The bar is 1m length (L=1m) and unit section.
The boundary conditions are:
The initial conditions are:
Material properties are: density and Young's modulus. With these parameters, the wave speed of the elastic material 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 so=4x105Pa and vo=1m/s. It will reflect at the fixed x=1m and the amplitude of stress at this point will be doubled to a value of 8x105 Pa, while the velocity of the wave after reflection will propagate along the bar with v = -1m/s.
4.2.1. Stability Analysis with Respect to the Smoothing Length
To apply the Taylor-SPH method, the bar has been discretized using 51 real particles and 50 virtual particles. The space between two consecutives real particles is. The space between real-virtual particles is. Here the Courant number, C, is defined as a function of Dx. It is clear that if the Courant number is defined as a function of Dx', the value of C would be doubled. The time step used for the analysis has been chosen to be which corresponds to a Courant number C=1.
The error estimation is computed using the L2 norm as well as the energy norm.
• The L2 norm of velocity error is obtained by:
where vexact is the exact solution and vhis the numerical solution given by Taylor-SPH.
• The error in the energy norm is obtained by
where sexact is the exact solution and sh is the numerical solution.
Table 1 gives the L2 norm of the velocity error and the error in the energy norm as function of ho/Dx. Fig. 6 shows the error in the L2 norm for different values of ho/Dx. It can be observed that for values of ho/Dx between 0.6 and 1.5, the results are in good agreement with the analytic solution. The error is the order of 10-3 % in the L2 norm and 10-8 % in energy norm. In this case, the number of real particles inside the compact domain is between 3 and 7 particles.
By increasing the value of ho/Dx≥1.6, the error in both norms increases. The numerical solution loses its accuracy and it becomes highly distorted. For higher values of ho/Dx, the solution diverges.
|ho/Dx||r = 2h0||Number of real particles inside the compact domain||Error L2 norm (%)||Error Energy norm (%)|
|0.6 - 1.5||0.024 - 0.06||3 - 7||10-3||10-8|
|1.6 - 2.5||0.072 - 0.1||7 - 11||13 - 27||5 - 8|
|Courant number||Dt(s)||Error L2 norm (%)||Error Energy norm (%)|
4.2.2. Stability Analysis with Respect to the Courant Number
In order to perform the stability analysis of the Taylor-SPH method with respect to the Courant number C, a similar example as above has been treated. The parameter ho/Dx has been kept constant equal to 1.5 and the courant number has been gradually increased from 0.1 to 1. The L2 norm of the error and the error in the energy norm as function of the Courant number are given in Table 2. Fig. 7 shows the error in the L2 norm for different values of C.
It can be observed that for Courant number C=1 which corresponds to time step, the accuracy of the solution is maximal with an error of about 10-3 % in L2 norm norm. As the value of the Courant number is decreased, the solution loses its accuracy and becomes oscillatory which results in an increasing error. For , the numerical solution diverge.
4.3. Shock Wave Propagation in Viscoplastic Damage Bar
In this Section, the shock wave propagation in elasto-viscoplastic damage bar with a length of L=2m has been studied. As mentioned above, the viscoplastic law is introduced to regularize the damage model. The parameters used in this example are: density r = 2400 kg/m3; Young’s modulus E=36000 MPa; viscoplastic model parameters (g,N)=(5,1); softening parameter H= -E/10; yield stress so=10MPa; damage parameters a =5000; ko=10-4.
The boundary conditions are the following: the bar is fixed at the right end x=L, and at the left end x=0 the velocity is imposed: where H is the Heaviside function and to=5x10-4s. The bar is disretized with 200 real particles and 199 virtual particles. The space between two consecutive real particles is Dx=0.01 m. The wave speed in the elastic material is c=3873m/s and the corresponding stress amplitude is. The Courant number chosen in the computation is C=1 which corresponds to the time step Dt=2.58x10-6 s. Fig. 8 shows the stress s(x=2m,t) at the right end of the bar after the reflection of the shock wave occurs for the following cases: (a) elastic bar (b) viscoplastic bar (c) viscoplastic damage bar. Before the reflection, the behaviour of the bar is elastic because the stress along the bar does not exceed the yield stress. After reflection the amplitude of the stress doubles and becomes . The bar behaves non-linearly and irreversible viscoplastic strains occur followed by damage of the material. After reflection of the wave, the stress exceeds the yield stress and therefore viscoplastic strain and damage are accumulated which results in the decrease of the stresses. The decrease of the curves depends on the parameters of viscoplastic and damage model.
In Fig. 9 is represented the evolution of damage variable D along the bar at different times. The maximum value of damage Dmax = 0.94 is reached at time t = 1.016 ms. The bar starts to be damaged at the fixed end x = 2m just after reflection of the wave. Then damage propagates to locate on a reduced distance of the bar. Fig. 10 shows the evolution of the damage at x = 2m versus time. After reflection of the wave, the stress value is doubled and the damage accumulates and converges to the limit value 0.94. Fig. 11 represents the stress-strain curves for two different refinements: 100 and 200 real particles. It can be observed that the results are similar and do not depend on the number of particles.
This example illustrates the performance of the Taylor-SPH method in numerical simulation of localization phenomena for viscoplastic damage material.
Modelling of shock wave propagating in linear or non linear materials presents several problems such as numerical diffusion and dispersion, volumetric locking, mesh dependency, accurate determination of failure surfaces, etc. Here, we have introduced a new meshless method Taylor-SPH that is based on a first-order system of PDEs involving velocities and stresses. It consists of a two-steps time discretization algorithm by means of a Taylor series expansion and a corrected SPH method for the spatial discretization. Two different sets of particles have been used for the numerical computation and a Lagrangian kernel has been used in order to avoid numerical instabilities. It is important to mention that using the two sets of particles ensures stability and not being necessary to introduce artificial viscosity in the model to solve the SPH tensile instability. This meshless method has been applied to shock waves propagating in 1D elastic bar and viscoplastic damage bar. From the obtained results, we can conclude that the main advantages of Taylor-SPH method are:
• It avoids the tensile instability
• It minimizes the numerical damping and dispersion when dealing with shock waves
• It is able to model localization phenomena for viscoplastic damage material without lose of hyperbolicity of PDEs.
• It requires a small number of particles to obtain accurate results.
We would like to thank Prof. M. I. Herreros for helpful discussions.