<figure><img src="image_1.png" /></figure>

HAL
open science

# Taylor-SPH meshless method

Hajar Idder, Mokhtar Mabssouta

► To cite this version:

Hajar Idder, Mokhtar Mabssouta. Taylor-SPH meshless method. 12e Colloque national en calcul des structures, CSMA, May 2015, Giens, France. hal-01517304

HAL Id: hal-01517304

https://hal.science/hal-01517304

Submitted on 3 May 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Taylor-SPH meshless method.

H. Idderª,*, M. Mabssoutª

ª Laboratoire de Mécanique et Génie Civil, Faculté des Sciences et Techniques, UAE, Tanger, Maroc
Corresponding author : *idderhajar@gmail.com

**Abstract** — In this work, we present the Taylor-SPH meshfree method for dynamic problem. The PDEs are written in terms of stress and velocity. A corrected Lagrangian kernel is employed for spatial discretization and two different sets of particles are used for time discretization. This method is used to solve the propagation of shock waves in elastic-viscoplastic media. The Taylor-SPH method is shown to be efficient and it provides solutions of good accuracy.

**Keywords**— Taylor-SPH, stability, shock wave, viscoplastic, localized deformation

## 1. Introduction

Finite Element Method (FEM) and Finite Volume Method have been widely used to solve many engineering problems in Solid and Fluid Mechanics. Despite of their great success, these 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. One of the oldest meshfree methods is the Smoothed Particle Hydrodynamics (SPH) [10,2]. The classical formulation of SPH method presents several difficulties when dealing with Dynamics and shock wave propagation. We have published in previous works some different alternatives to solve the propagation of shock waves and localized strain in viscoplastic solid using FEM [6,7,8]. In this paper, a new time discretization algorithm using SPH for solving the propagation of shock waves in viscoplastic media is presented. This method is based on two points; discretization in time based on the Taylor expansion in two steps and discretization in space based on the corrected SPH method. The Taylor-SPH method is applied to solve the problem of propagation of shock waves in an elastic-viscoplastic media.

The paper is organized as follows. First the governing equations for dynamic problems in viscoplastic media are given in terms of stress and velocity. In section 3, numerical discretization is presented using the Taylor-SPH method. Some numerical examples in 1D and 2D are described in Section 4.

## 2. Governing equations

Neglecting the body forces, and for small strains, the balance of momentum and constitutive equations can be written in terms of stress s and velocity v as follows:

- Equilibrium equation

$$ \rho \frac{\partial v}{\partial t} = \text{div} \, \sigma \qquad (1) $$

- Constitutive equation

It will be assumed that the material behaviour can be described by Perzyna's viscoplastic law [9]:

$$ \frac{\partial \sigma}{\partial t} = D^e : \left( \frac{\partial \varepsilon}{\partial t} - \frac{\partial \varepsilon^{vp}}{\partial t} \right) \qquad (2) $$

Where $\rho$ is the density, $D^e$ is the elastic constitutive tensor, $\varepsilon$ the strain tensor and $\varepsilon^{vp}$ is the viscoplastic strain given by the Perzyna law [9]:

$$ \dot{\varepsilon}_{vp} = \bar{\gamma} \left\langle \left( \frac{F - \bar{\sigma}}{\bar{\sigma}} \right)^N \right\rangle \frac{\partial F}{\partial \sigma} \qquad (3) $$

In (3), <.> is the Macaulay brackets and $\gamma$ and $N$ are model parameters. $F$ is a function describing a convex surface in the stress space. Von Mises yield criterion has been used. Equations (1) and (2) can be written as:

$$ \frac{\partial U}{\partial t} + \text{div}F = S \qquad (4) $$
Being $U^T = (\sigma_{11} \ \sigma_{22} \ \sigma_{12} \ v_1 \ v_2)$ the unknown vector, and $F$ and $S$ the flux and source terms respectively.

$$ 
\begin{aligned}
F_x^T &= (-D_{11}v_1 - D_{12}v_2 - D_{33}v_2 - \frac{\sigma_{11}}{\rho} - \frac{\sigma_{12}}{\rho}) \\
F_y^T &= (-D_{12}v_{21} - D_{22}v_2 - D_{33}v_1 - \frac{\sigma_{12}}{\rho} - \frac{\sigma_{22}}{\rho}) \\
S^T &= (-D_{11}\dot{\epsilon}_1^{VP} - D_{12}\dot{\epsilon}_2^{VP} - D_{12}\dot{\epsilon}_1^{VP} - D_{22}\dot{\epsilon}_2^{VP} - D_{33}\dot{\epsilon}_3^{VP} \quad 0 \quad 0)
\end{aligned}
\quad (5) $$

### 3. Taylor-SPH method (TSPH)

The Taylor-SPH method [3,4,5] is used to solve the PDEs (4), and it 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. Taylor-SPH: Time discretization

Time discretization of equation (4) is carried out by means of a Taylor series expansion in time of $U$ up to second order of accuracy in two steps [7, 8]:

$$ \textit{First step: } U^{n+1/2} = U^n + \frac{\Delta t}{2}(S - \Delta F)^n \quad (6a) $$

$$ \textit{Second step: } U^{n+1} = U^n + \Delta t(S - \Delta F)^{n+1/2} \quad (6b) $$

#### 3.2. Taylor-SPH: Spatial discretization

The Taylor-SPH spatial discretization is carried out using two steps and two sets of particles named “real and virtual” particles [3,4,5].

*First Step:* Applying the corrected SPH spatial discretization to (6a), we obtain:

$$ U_{vp}^{n+1/2} = U_{vp}^{n} + \frac{\Delta t}{2} \left( \sum_{j=1}^{Nr} \frac{m_j}{\rho_j} S_j^n \tilde{W}_{ij} - \sum_{j=1}^{Nr} \frac{m_j}{\rho_j} F_j^n \tilde{\nabla} w_{ij} \right) \quad (7a) $$

$$ \text{Where } \tilde{W}_{ij} = \frac{W_{ij}}{\sum_{j=1}^{N} W_{ij} \Omega_j}; \quad \tilde{\nabla} W_{ij} = \left(\sum_{j=1}^{N} X_{ij} \otimes \nabla W_{ij} \Omega_j\right)^{-1} \nabla W_{ij} $$

The subscript VP refers to the virtual particles. $W$ is the kernel function and $ho$ is the smoothing length that defines the size of the kernel support. $J$ is the “real” particles, such that $|X_J - X_{vp}| \le 2h_0$. $m_J$ and $\rho_J$ are the mass and density associated to particle $J$ and $\Omega_J$ is the volume associated to particle $J$.

*Second Step:* Applying the corrected SPH spatial discretization to equation (6b), we obtain

$$ U_{Rp}^{n+1} = U_{Rp}^{n} + \Delta t \left( \sum_{j=1}^{Nr} \frac{m_j}{\rho_j} S_j^{n+1/2} \tilde{W}_{ij} - \sum_{j=1}^{Nr} \frac{m_j}{\rho_j} F_j^{n+1/2} \tilde{\nabla} w_{ij} \right) \quad (7b) $$

where RP= Real particles and J are the “virtual” particles such that $|X_J - X_{Rp}| \le 2h_0$.

### 4. Numerical examples

#### 4.1. 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 problem of propagation of a velocity
shock wave in a 1D elastic bar. The bar is 1m length (L=1m) and unit section. The B-spline function [1], with a support domain of radius 2h₀, has been considered.

### 4.1.1. Stability analysis with respect to the smoothing length $h_0$

It is well known that the smoothing length $h$ is very important in the SPH method since it has direct influence on the efficiency of the computation and the accuracy of the solution. If $h$ is too small, there may be not enough particles in the support domain of dimension $\kappa h$ to exert forces on a given particle, which results in very low accuracy. On the contrary, if the smoothing length is too large, local properties may be smoothed out, and the accuracy decreases too. Therefore the particles approximation used by the SPH method depends on having a sufficient and necessary number of particles within the support domain of $\kappa h$.

To apply the Taylor-SPH method, the bar has been discretized using 50 “real” particles separated from each other a distance of $\Delta x=0.02$ m. The time-step used for the analysis has been chosen to be $\Delta t=2.10^{-4}$ s, which corresponds to a Courant number $C = \frac{c\Delta t}{\Delta x} = 2$, where c is the wave speed. Fig.1 (a)

shows the error in the L2 norm of the velocity error for different values of $h_0/\Delta x$. It can be observed that for values of $h_0/\Delta x$ between 0.6 - 1.5, the results are in good agreement with the analytic solution. And the error is the order of $10^{-8}$ %. By increasing the value of $h_0/\Delta x \ge 1.8$, the error increases. The numerical solution loses its accuracy and it becomes highly distorted. For higher values of $h_0/\Delta x$, the solution diverges.

### 4.1.2. Stability analysis with respect to the Courant number

In order to accomplish a sensitivity analysis of Taylor-SPH (TSPH) with respect to the time-step, $\Delta t$, a similar analysis as above has been carried out. Therefore the same problem of the propagation of a shock wave in an elastic bar of length L = 1 m has been solved considering a fixed distribution of 50 “real” particles separated from each other a distance of $\Delta x = 0.02$ m. The parameter $h_0/\Delta x$ has been chosen to be equal to 1.5. The value of parameter $\Delta t$ has been gradually increased and the accuracy and stability of the solution for different values of the Courant number have been studied.

Fig.1 (b) shows the error in the L2 norm for different values of C. It can be observed that for Courant number $C=2$, the accuracy of the solution is maximal with an error of about $10^{-7}$ %. As the value of the Courant number is decreased, the solution loses its accuracy and becomes oscillatory which results in an increasing error. For $C < 2$, the numerical solution diverge.

* The L2 norm of velocity error is obtained by:

$$ \|Err\|_{L_2} = \frac{\|v^{exact} - v^h\|_{L_2}}{\|v^{exact}\|_{L_2}} \quad \text{with} \quad \|v^h\|_{L_2} = \sqrt{\int_0^L v^2 dx} \quad (8) $$

Where $v^{exact}$ is the exact solution and $v^h$ is the numerical solution given by Taylor-SPH.

## 4.2. Propagation of a shock wave on a 2D viscoplastic bar

This problem consists of two opposite velocity shockwaves propagate into a bi-dimensional bar. The bar is 1m long and its cross section has a diameter of 0.1m. 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.

Fig. 2(a) shows the viscoplastic strain at the middle of the bar where the stress reaches the yield surface. The stretching of the bar causes the narrowing at the middle. This example shows that the Taylor-SPH algorithm is capable to describe and reproduce localized modes of failures. Fig. 2(b) illustrates the viscoplastic strain along the bar using the Taylor-SPH algorithm for two different refinements: 306 and 1111 real particles. It can be observed that the viscoplastic strain is localized at the middle of the bar and the results do not depend on the number of particles used in the simulation.
<figure><img src="image_2.png" /><figcaption>Figure 1: (a) Error in L2 norm for different values of Courant number
(b) Error in L2 norm for different values of $h_0/\Delta x$</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2: (a) Localized viscoplastic strain at the middle of the bar with $\gamma = 50$ and $N = 1$
(b) Influence of the number of particles: viscoplastic strain along the bar</figcaption></figure>

**5. Conclusion**

The Taylor-SPH method applied to shock waves propagating in viscoplastic media has been presented. The results show that the Taylor-SPH algorithm is capable to describe and reproduce localized viscoplastic strain. The Taylor-SPH method has been proved to be stable, robust and efficient, requiring only a reduced number of particles to obtain accurate results.

**6. References bibliographiques**

[1] J.J. Monaghan, J.C. Lattanzio. A refined particle method for astrophysical problems, Astronomy and Astrophysics, 149: 135-143, 1985.

[2] L.B. Lucy. A numerical approach to the testing of fusion process. Astronomical J, 82:1013-1024, 1977.

[3] M. Mabssout & M.I. Herreros. Runge-Kutta vs Taylor-SPH: Two time integration schemes for SPH with application to Soil Dynamics. Applied Mathematical Modelling, 37:3541-3563, 2013.

[4] M.I. Herreros & M. Mabssout. A two-steps time discretization scheme using the SPH method for shock wave propagation. Comput. Methods Appl. Mech. Engrg, 200: 1833-1845, 2011.

[5] M. Mabssout & M.I. Herreros. Taylor-SPH vs Taylor-Galerkin for shock waves in viscoplastic continua. European Journal of Computational Mechanics, 20/5-6: 281-308, 2011.

[6] M. Mabssout & M. Pastor M., M.I. Herreros, M. Quecedo. A Runge-Kutta, Taylor-Galerkin scheme for hyperbolic systems with source terms. Application to shock wave propagation in viscoplastic geomaterials. Int. J. Numer. Analytical Meth. Geomechanics, 30(13) 1337-1355, 2006.

[7] M. Mabssout & M. Pastor. A two step Taylor-Galerkin algorithm for shock wave propagation in soils. Int. J. Numer. Analytical Meth. Geomechanics, 27: 685-704, 2003.

[8] M. Mabssout & M. Pastor. A Taylor-Galerkin algorithm for shock wave propagation and strain localization failure of viscoplastic continua. Comput. Methods. Appl. Mech. Engrg, 192:955-971, 2003.

[9] P. Perzyna. Fundamental problems in viscoplasticity. Recent Advances in Applied Mechanics, Academic press, New York, 9: 243-377, 1966.

[10] R.A. Gingold & J.J.Monaghan. Smoothed particles hydrodynamics: Theory and application to nonspherical Stars. Monthly Notices of the Royal Astronomical Society, 181:375-389, 1977.